How to implent a time dependent poisson eqution with a source term, that is described by a differential equation?

[Cou]

The source term f is described by the following equation:

\frac{d}{dt}f(t) = \sqrt{f(t)} * \left ( constant - u_{average\;of\;a\;facet}(x_1,t) \right ),

where u_{average\;of\;a\;facet}(x_1,t) is the average value of u of one of the facets of the domain. Because I currently don‘t know how the average value of a facet is calculated in Fenics, I used the average of u of the whole domain as a dummy value to realize the minimum working code example (MWCE).

The minimum working code example is based on the ft03_heat.py tutorial, which explains a time dependent poisson equation and is governed by the following equations:

{u}'=\Delta u + f \;\;\; in the unit square
u = u_D \;\;\;\;\;\;\;\;\;\;\; on the boundary (Dirichlet boundary conditions)
u = u_0 \;\;\;\;\;\;\;\;\;\;\;\; at t = 0

So, how do I implement the equation of the source term f in Fenics, without using additional software?
Or, in other words, how do I replace the pseudo code in the MWCE underneath with real Fenics code?

from __future__ import print_function
from fenics import *

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

constant = 5       # known Constant to calculate f


# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = 3 # known initial condition of f: f(t=0) = 3

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    
    # Calculate average u of the domain instead of calculate the average of a facet
    vertex_values_u = u.compute_vertex_values(mesh)
    avg_u = np.average(vertex_values_u)
    
    # Pseudo code to calculate f for the current time step
    equation = df/dt - sqrt(f) * (constant - avg_u)  == 0    
    f = solve(equation, bc = f_0)
    
    # Compute solution
    solve(a == L, u, bc)

    # Plot solution
    plot(u)

    # Update previous solution
    u_n.assign(u)

# Hold plot
interactive()

Thanks,
Jan

[A.Sourais]

Hello,
This approach from @kamensky might be helpful :

Combine a set of PDE with a fenics mechanics package

One approach that is straightforward, albeit perhaps difficult to optimize for large problems, is to just add the ODE system unknowns as “Real”-type components of a mixed space, i.e., global constants. Consider a contrived example (which I just made up for a code demonstration, rather than to be a reasonable model any physical system). The PDE problem is: Find u:\Omega\to\mathbb{R} such that \begin{array}{lr} \partial_t u - c_1\Delta u = 0 & \text{ on }\Omega\text{ ,}\\ c_1\partial u/\parti…

[Cou]

Hello,

thanks for introducing me to the concept of a mixed function space .

In order to realize a faster implementation, I changed my initial set of PDEs, , which are more similar to the ones @kamensy used, to the following:

Poisson equation:
u′=Δu+f in the unit square
\frac{d}{dt} f=−\frac{1}{|∂Ω|} \int_{∂Ω}^{}u
f_D=\sin(\pi x) * \sin(\pi y)
u=u_D on the boundary (Dirichlet boundary conditions)
u_D=1+x² + 3y² + 1.2t
u=u_0 at t = 0

Unfortunaly, due to the mixed functionspace, I´m no longer able to interpolate u_D:

u_D = Expression(‘1 + x[0] * x[0] + 3 * x[1] * x[1] + 1.2 * t’,
degree=1, alpha=alpha, beta=beta, t=0)

u_n = interpolate(u_D, V)

because

*** Error: Unable to interpolate function into function space.
*** Reason: Rank of function (0) does not match rank of function space (1).

How is it possible to fix this?

Here is my updated MWCE:

from fenics import *

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
cell = mesh.ufl_cell()

uE = FiniteElement('CG',cell,1)
fE = FiniteElement('R',cell,0)

# Define variational problem
V = FunctionSpace(mesh, MixedElement([uE,fE]))

uf = TrialFunction(V)
uf_old = Function(V)

vr = TestFunction(V)
u,f = split(uf)
u_old,f_old = split(uf_old)
v,r = split(vr)

# Define general paramter
T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Define boundary condition
x = SpatialCoordinate(mesh)
u_D = Expression('1+x[0]*x[0] + alpha*x[1]*x[1] + beta*t', 
                 degree=2, alpha=alpha, beta=beta, t=0)

def boundary(x, on_boundary):
    return on_boundary

# Define initial value
u_n = interpolate(u_D, V)

bc = DirichletBC(V, u_n, boundary)

# ODE for source term f driven by average value of PDE unknown u
f_t = (f-f_old)/dt
f = (f_t - u)*r*ds

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - u_n*v*dx + dt*f
a, L = lhs(F), rhs(F)

# Time-stepping
ufh = Function(V)
uf_old.assign(project(as_vector([sin(pi*x[0])*sin(pi*x[1]),Constant(1.0)]),V))
g_out = File("u.pvd")
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    
    # Compute solution
    solve(a == L, ufh, bc)
    
    uf_old.assign(ufh)
    u_out = ufh.split()[0]
    u_out.rename("u","u")
    g_out << u_out

    # Update previous solution
    u_n.assign(u)

Thanks

[A.Sourais]

Dear @Cou,
interpolate expects an expression that has equal number of components with V . You can define your Dirichlet boundary conditions without the need of interpolating like this:

# Define boundary condition
x = SpatialCoordinate(mesh)
u_D = Expression('1+x[0]*x[0] + alpha*x[1]*x[1] + beta*t', 
                 degree=2, alpha=alpha, beta=beta, t=0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V.sub(0), u_D, boundary)

where V.sub(0) is the first sub space reffering to u in your notation.

Additionally you can use expressions directly in ufl forms:

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - u_D*v*dx + dt*f
a, L = lhs(F), rhs(F)

I hope this was helpful .

[Cou]

Thanks @A.Sourais, that helped a lot.

Here is my running code, in case somebody else faces a similar problem:

from fenics import *

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
cell = mesh.ufl_cell()

uE = FiniteElement('CG',cell,1)
fE = FiniteElement('R',cell,0)

# Define variational problem
V = FunctionSpace(mesh, MixedElement([uE,fE]))

uf = TrialFunction(V)
uf_old = Function(V)

vr = TestFunction(V)
u,f = split(uf)
u_old,f_old = split(uf_old)
v,r = split(vr)

# Define boundary condition
x = SpatialCoordinate(mesh)
u_D = Expression('1+x[0]*x[0] + alpha*x[1]*x[1] + beta*t', 
                 degree=1, alpha=alpha, beta=beta, t=0)

    
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V.sub(0), u_D, boundary)

# ODE for source term f driven by average value of PDE unknown u
f_t = (f-f_old)/dt
f = (f_t - u)*r*ds

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - u_D*v*dx + dt*f
a, L = lhs(F), rhs(F)

# Time-stepping
ufh = Function(V)
uf_old.assign(project(as_vector([sin(pi*x[0])*sin(pi*x[1]),Constant(1.0)]),V))
g_out = File("u.pvd")
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    
    # Compute solution
    solve(a == L, ufh, bc)
    
    uf_old.assign(ufh)
    u_out = ufh.split()[0]
    u_out.rename("u","u")
    g_out << u_out