Interpolate not working with external library

1

Hello,

I’m attempting to incorporate the external library CoolProp to calculate fluid properties to be used in a finite element Navier Stokes algorithm. I’m doing this in a pretty simplified manner, see code below. It seems that when I try to interpolate the value to create a function object, the values are no longer the same. Not sure what’s going wrong exactly, if someone could point out my error that would be greatly appreciated.

from fenics import *
from dolfin import *    
from mshr import * 
import numpy as np   
from CoolProp.CoolProp import PropsSI

mesh = UnitSquareMesh(100, 100)

Q = FunctionSpace(mesh, "CG", 1)


class ThermalProperty(UserExpression):
    """
    ThermalProperty(property, a, A, b, B)

    property = string for property sought for
    a = string for input state variable
    A = FEniCS function, input state variable function
    b = string for second input state variable
    B = FEniCS function for second input state variable

    """

    def __init__(self, property, a, A, b, B, **kwargs):
        super().__init__(**kwargs)
        self.property = property
        self.a = a
        self.b = b
        self.A = A
        self.B = B
        self.list = []

    def eval(self, value, x):
        value = PropsSI(self.property, self.a, self.A(x),
            self.b, self.B(x), "CarbonDioxide")
        self.list.append(value)
        return

    def value_shape(self):
         return ()


temp_0 = Expression("320", degree = 0)
T_n = interpolate(temp_0, Q)
p_n = interpolate(Constant(52e6), Q)
rho_n = interpolate(Constant(978.15), Q)


C_V = ThermalProperty("CPMASS", "DMASS", rho_n, "T", T_n)

print(interpolate(C_V, Q).vector().get_local())
print(np.array(C_V.list))

[2.0895309e-316 2.0895309e-316 2.0895309e-316 … 2.0895309e-316
2.0895309e-316 2.0895309e-316]
[1741.70968994 1741.70968994 1741.70968994 … 1741.70968994 1741.70968994
1741.70968994]

2

What you want to do is assign the 0-th element of value in the eval method. See, e.g., the following example:

from dolfin import *

class UEx(UserExpression):
    def eval(self,value,x):
        
        # Incorrect: Creates local variable that is no longer connected
        # to function argument.
        #value = 1.0

        # Correct: Overwrites number stored in the input array.
        value[0] = 1.0
        
    def value_shape(self):
        return ()

mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh,"CG",1)
print(assemble(interpolate(UEx(degree=0),V)*dx))
3

Very simple error, thank you very much for pointing that out.