variation form dependent on non-analytical/complicated expression

Question

I am working on a nonlinear temperature field problem; some of the terms variational/weak form depend on the temperature. I plan on implementing solving the nonlinear problem "manually", using a simple iteration scheme: guess a temperature, solve, and update parameters, and solve again until I reach the converge criteria. My problem is transient so this scheme works nicely as you can use the previous temperature field as a good guess of the temperature field for the next time step.

I've seen examples on how to treat nonlinear terms but those cases have the simple expressions, like k = 1+u**2. For my problem, I am relying on a python package that provides an equation of state which tells me what a certain property is given the temperature. The equation of state does not have some analytical or polynomial form and there's a lot going on in the background. How do I go about implementing this in FEniCS? Do I need to create an instance of Expression that has a method that determines the value of the parameter of interest given the temperature value of the node/vertex?

Answer 1

I'm pretty sure that is what you are looking for -> building your own custom expressions. The following example is for a modified custom 3D expression used on the RHS as source in my weak form. You can use your own definitions like my topo-function for "a lot is going on in the background" within the dolfin expression.

import numpy as np
import dolfin as df

def topo(x, y, h=250., xmin=-1e4):

    return h * np.sin((10.0 / xmin) * y) + 0.1 * x

class Custom_Expression(df.Expression):

    def __init__(self, J=1., a=6., origin=[1., 2., -5.], R=100.,
                        **fenics_kwargs):

    self.J = J
    self.a = a
    self.R = R
    self.O = origin

def eval(self, values, x):

    xx = (x[0] - self.O[0])
    yy = (x[1] - self.O[0])

    # you can use either dolfin math ...

    r = df.sqrt(xx**2 + yy**2)
    values[0] = (-(self.J * (x[1]/r))) * df.exp(-(r-self.R)**2 / self.a**2)
    values[1] = ((self.J * (x[0]/r))) * df.exp(-(r-self.R)**2 / self.a**2) 
    values[2] = 0.0

    # ... or something such as own defenitions or numpy math:

    values[0] = x[0]
    values[1] = x[1] * np.sin(x[0])
    values[2] = x[2] * topo(xx, yy)

def value_shape(self):  #  not necessarily needed, at least in Fenics 2016.2
    return (3,)

FEniCS needs some keyword argument such as the element type, so call your expression in every time step t in the weak form as in the following pseudo-code. The custom expression is evaulated on the fly in every time-step you call it with your updated input parameters.

mesh = df.UnitCubeMesh(6, 6, 6)
FS = df.VectorFunctionSpace(mesh, 'CG', 1)
u = df.TrialFunction(FS)
v = df.TestFunction(FS)
J_t = 1.0
R_t = 234.5
C_E_t = Custom_Expression(J=J_t, R=R_t, element=FS.ufl_element())

lhs = df.inner(df.grad(u), df.grad(v))*dx
rhs = df.inner(C_E_t, v)*dx
# solve(lhs, rhs)