# Solving a system of PDE's in fenicsx

Dear all,

I’m trying to solve the porous media equation by means of a sort of L-scheme.

Practically this means I’m trying to solve this system of coupled equations in the weak form:

(u_n^i -u_{n-1}, \phi) + \tau(\nabla w_n^i, \nabla \phi) = 0 \\ L(u_n^i - u_n^{i-1}, \xi) = (w_n^i - \Phi(u_n^{i-1}), \xi)

In this case, \Phi(u) = u^m for some m>1 and L>0 big enough, and this system will solve for \partial_{t}u = \Delta u^m. \tau is the time-step discretization, and for enough iterates in i, this will converge to the time-discretized solution. For our finite element solution, we need \phi \in CG1 and \xi \in DG1.

Here I copied pasted the relevant parts of my try so far (for one iteration of this scheme):

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.mesh import locate_entities, meshtags
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from dolfinx import io

# Creating an interval as mesh
h = 0.1
x_min, x_max = -2,2
nx = int((x_max - x_min)/h)
domain = mesh.create_interval(MPI.COMM_WORLD, nx, points = (x_min, x_max))
x = ufl.SpatialCoordinate(domain) # Coordinates of domain used for defining general functions

el1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
el2 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Z = FunctionSpace(domain, mel)

V = FunctionSpace(domain, ("CG", 1))
U = FunctionSpace(domain, ("DG", 1))

(w_trial, u_trial) = ufl.TrialFunctions(Z)
(w_phi, u_phi) = ufl.TestFunctions(Z)
Z_sol = fem.Function(Z)

# Constants
m  = 2
alpha_par = 1/3
beta_par = 1/3
kappa_par = 1/12
t_start  = 0.3

def Phi(u):
return u**m

def exact_sol(x,t):
return ufl.operators.max_value(0,1/(t**alpha_par) * (C - kappa_par * (ufl.operators.sqrt(x[0]**2)/(t**beta_par))**2 )**(1/(m-1)))

initial_u = exact_sol(x,t_start)
expr_u = fem.Expression(initial_u, U.element.interpolation_points())

# Defining the function spaces for all the functions
u_n = fem.Function(U)
u_n.name = "u_n"

u_n_i = fem.Function(U)
u_n_i.name = "u_n_i"

w_n = fem.Function(V)
w_n.name = "w_n"

w_n_i = fem.Function(V)
w_n_i.name = "w_n_i"

# Defining the solution variables
uh = fem.Function(U)
uh.name = "uh"

wh = fem.Function(V)
wh.name = "wh"

initial_u = exact_sol(x,t_start)
expr_u = fem.Expression(initial_u, U.element.interpolation_points())
expr_w = fem.Expression(Phi(initial_u), V.element.interpolation_points())

u_n.interpolate(expr_u)
u_n_i.interpolate(expr_u)
w_n.interpolate(expr_w)
w_n_i.interpolate(expr_w)

L_par = 1
L = fem.Constant(domain, ScalarType(L_par))

a = (u_trial * w_phi + dt*ufl.inner(ufl.grad(w_trial), ufl.grad(w_phi)) + L* u_trial * u_phi - w_trial * u_trial) * ufl.dx
L = (u_n * w_phi + L * u_n_i * u_phi - Phi(u_n_i) * u_phi) * ufl.dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)

b = fem.petsc.assemble_vector(linear_form)
b.assemble()
A = fem.petsc.assemble_matrix(bilinear_form)
A.assemble()

solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY) # PREONLY means you only use the preconditioner
solver.setTolerances(rtol=1e-12, atol = 1e-60)
solver.getPC().setType(PETSc.PC.Type.LU) # The LU preconditioner just solves the system using LU factorization.

solver.setOperators(A)

solver.solve(b, Z_sol.vector)
(wh, uh) = Z_sol.split()


The only thing I’m not sure on, is if I have set up my bilinear form up correctly, and if this is indeed the way to solve these kinds of systems. I could not find any tutorials or information on how to do it in dolfinx. If I run this I get an error at bilinear_form = fem.form(a), which says TypeError: cannot unpack non-iterable bool object.

Does anyone know how to fix it, or alternatively know of a different way of solving two PDE’s as a simultaneous system?

Many thanks,
Robin

I don’t have a specific answer to the TypeError but you might find the Cahn-Hilliard demo of value.

Also, doesn’t the function \Phi(u) make your system non-linear, but you are solving it as a linear problem?

1 Like

Thank you. I had already looked at the Cahn-Hilliard demo but I want to solve the system using petsc so I want to construct these fem.forms and am not quite sure if this is the way to do it, in this case. I’m also not quite sure what the difference would be between (if there is any)

el1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
el2 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Z = FunctionSpace(domain, mel)


and I guess equivalently:

P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
DG1 = ufl.FiniteElement("DG", msh.ufl_cell(), 1)
ME = FunctionSpace(msh, P1 * DG2)


You are correct that \Phi(u) makes the PDE non-linear, but the idea of this scheme is that we linearize the system through this second equation, as the non-linearity is now only depenend on a previous iteration.

For those encountering a similar problem, I think to have solved it. I just made a small error in the last term of the bilinear form where I wrote w_trial * u_trial where I meant to write w_trial * u_phi

This code should work

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.mesh import locate_entities, meshtags
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from dolfinx import io

# Creating an interval as mesh
h = 0.1
x_min, x_max = -2,2
nx = int((x_max - x_min)/h)
domain = mesh.create_interval(MPI.COMM_WORLD, nx, points = (x_min, x_max))
x = ufl.SpatialCoordinate(domain) # Coordinates of domain used for defining general functions

el1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
el2 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Z = FunctionSpace(domain, mel)

V = FunctionSpace(domain, ("CG", 1))
U = FunctionSpace(domain, ("DG", 1))

(w_trial, u_trial) = ufl.TrialFunctions(Z)
(w_test, u_test) = ufl.TestFunctions(Z)
Z_sol = fem.Function(Z)

# Constants
m  = 2
alpha_par = 1/3
beta_par = 1/3
kappa_par = 1/12
C = 1/32
t_start  = 0.3
dt = 0.1

def Phi(u):
return u**m

def exact_sol(x,t):
return ufl.operators.max_value(0,1/(t**alpha_par) * (C - kappa_par * (ufl.operators.sqrt(x[0]**2)/(t**beta_par))**2 )**(1/(m-1)))

initial_u = exact_sol(x,t_start)
expr_u = fem.Expression(initial_u, U.element.interpolation_points())

# Defining the function spaces for all the functions
u_n = fem.Function(U)
u_n.name = "u_n"

u_n_i = fem.Function(U)
u_n_i.name = "u_n_i"

w_n = fem.Function(V)
w_n.name = "w_n"

w_n_i = fem.Function(V)
w_n_i.name = "w_n_i"

# Defining the solution variables
uh = fem.Function(U)
uh.name = "uh"

wh = fem.Function(V)
wh.name = "wh"

initial_u = exact_sol(x,t_start)
expr_u = fem.Expression(initial_u, U.element.interpolation_points())
expr_w = fem.Expression(Phi(initial_u), V.element.interpolation_points())

u_n.interpolate(expr_u)
u_n_i.interpolate(expr_u)
w_n.interpolate(expr_w)
w_n_i.interpolate(expr_w)

L_par = 1
L_const = fem.Constant(domain, ScalarType(L_par))

a = (u_trial * w_test + dt * ufl.inner(ufl.grad(w_trial), ufl.grad(w_test)) + L_const * u_trial * u_test - w_trial * u_test) * ufl.dx
L = (u_n * w_test + L_const * u_n_i * u_test - Phi(u_n_i) * u_test) * ufl.dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)

b = fem.petsc.assemble_vector(linear_form)
b.assemble()
A = fem.petsc.assemble_matrix(bilinear_form)
A.assemble()

solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY) # PREONLY means you only use the preconditioner
#solver.setTolerances(rtol=1e-12, atol = 1e-60)
solver.getPC().setType(PETSc.PC.Type.LU) # The LU preconditioner just solves the system using LU factorization.

solver.setOperators(A)

solver.solve(b, Z_sol.vector)
(wh, uh) = Z_sol.split()


Also for those looking for more info on these types of things, Stokes equations with Taylor-Hood elements — DOLFINx 0.6.0 documentation is a very nice demo dealing with similar issues.

@nav thanks for the help !

1 Like