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:
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
!