How to handle Nonlinear PDE with Division and Exponent?

1

Hello, I’m fairly new to Fenics and am trying to simulate the following (modified) reaction-diffusion

\frac{\partial{T(x, t)}}{\partial{t}} = k \frac{\partial^2{T(x, t)}}{\partial{x^2}} - \rho \frac{\partial{T(x, t)}}{\partial{x}} + \alpha S(x,t) \exp(-\frac{\beta}{T(x,t)}) -\alpha \gamma T(x,t) \\ \frac{\partial{S(x, t)}}{\partial{t}} = -\gamma_s S(x,t) \exp(\frac{-\beta}{T(x,t)})

I’ve written both equations in weak form, but am having trouble implementing on Fenics. The nonlinear terms, specifically, have given me a lot of trouble since division and exponentiation is not allowed for function types. Are there any specific ways to go about this? I couldn’t seem to find anything in the tutorials.

For reference, here’s my attempt at writing up the weak form as a function.

# Create the Mesh
mesh = dolfin.IntervalMesh(n, xgrid[0], xgrid[-1])
T = dolfin.FunctionSpace(mesh, 'CG', 2)
S = dolfin.FunctionSpace(mesh, 'CG', 2)


# FEM Discretization
v = dolfin.TestFunction(T)
phi = dolfin.TrialFunction(T)

w = dolfin.TestFunction(S)
theta = dolfin.TrialFunction(S)

clsdform_T = alpha * dolfin.inner(S * **expterm**, v)*dolfin.dx

clsdform_S = -gamma_s * dolfin.inner(S * **expterm**, w)*dolfin.dx

where the expTerm is

\exp(\frac{-\beta}{T(x, t)})

and T and S are the unknowns I’m attempting to solve for, while \alpha, \gamma_s, \beta are constants, and v and w are my test functions.

Thanks for the help!

EDIT Included only relevant parts of code relating to original question and added the relevant PDE.

2

There are lots of things which don’t make sense in your code. I recommend you examine the demos. Particularly for your problem the nonlinear Poisson demo.

3

Thanks for the suggestion @nate - I’m looking through it now. Perhaps I should have left out the code since it’s just a sketch. I guess the heart of my question lies at constructing exponent exp and division operators / when defining the weak form. Is that even possible or is there some other clever way of getting around it?

4

Hi @jarvin,

A few helpful (hopefully) points to consider: first, it’s recommended to post minimal working code examples: Read before posting: How do I get my question answered?. Running this code

gives

  File "/home/connor/Documents/fenics/code.py", line 14
    clsdform_T = alpha * dolfin.inner(S * **expterm**, v)*dolfin.dx
                                          ^
SyntaxError: invalid syntax

It’s easiest for people to assist if they don’t have to debug your code. :slight_smile: It’s also recommended to post the full error message that you receive.

Second, don’t confuse the functions T and S with the function spaces on which they live. In your code, you have defined S and T as dolfin.FunctionSpaces. I’m not an expert, but I wouldn’t expect UFL operations to be supported on function spaces. Instead, you should use the test and trial functions in your weak form.

Exponentiation (using the UFL operator dolfin.exp) and division seem to work for me when I replace the FunctionSpaces S and T with the trial and test functions v, w, phi, and theta defined on those function spaces:

import dolfin

alpha = 1
beta = 1
gamma_s = 1

# Create the Mesh
mesh = dolfin.IntervalMesh(10, 0, 1)
T = dolfin.FunctionSpace(mesh, 'CG', 2)
S = dolfin.FunctionSpace(mesh, 'CG', 2)

# FEM Discretization
v = dolfin.TestFunction(T)
phi = dolfin.TrialFunction(T)

w = dolfin.TestFunction(S)
theta = dolfin.TrialFunction(S)

clsdform_T = alpha * dolfin.inner(theta*dolfin.exp(-beta/phi), v)*dolfin.dx
clsdform_S = -gamma_s * dolfin.inner(phi*dolfin.exp(-beta/theta), w)*dolfin.dx

(Please note that I have not attempted to determine the correct weak form for this problem.)

A final note: since you are solving a nonlinear problem, you should use Functions (not TrialFunctions) in your weak form, as the nonlinear Poisson demo explains.

2 Likes
5

Hi @conpierce8, thanks for the response. I apologize for the messy initial post - I’ll keep your suggestions in mind next time. FWIW, I actually have a Tfun = dolfin.Function(T) and Sfun = dolfin.Function(S) and then define the closed form as shown above but replacing T and S as Tfun and Sfun, respectively.

Additionally, I believe I’ve figured out how to write the nonlinear term (code ran without runtime error). I think my initial problem was that I was defining T as a VectorFunctionSpace (it shows FunctionSpace in my code above, but that was because I was in the middle of debugging) which is an error since T is supposed to be a scalar at all points in my mesh. Once I changed it to FunctionSpace, it ran just fine with the following…

clsdform_T = alpha * dolfin.inner(Sfun * dolfin.exp(-beta / Tfun), v)*dolfin.dx

Here’s to hoping the solution is correct! Cheers and thanks again for the suggestions!

1 Like