Hi,
I have a working example that solves the simple linear PDE \partial_t T - \nabla^2T - r = 0 using the method of lines, where petsc ts is used (see my mwe below):
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, log
from dolfinx.io import XDMFFile
from dolfinx.fem import (functionspace, Function, locate_dofs_topological, dirichletbc, form)
import dolfinx.fem.petsc as petsc
from basix.ufl import element
from ufl import TestFunction, dx, grad, inner, TrialFunction
from petsc4py import PETSc
# 3d geometry parameters
box_length = 1
cells_number_per_dim = 4
# time stepping
t = 0
t_final = 10
dt_init = 1.
msh = mesh.create_box(MPI.COMM_WORLD,[[-box_length/2, -box_length/2, -box_length/2], [box_length/2,box_length/2, box_length/2]],
[cells_number_per_dim, cells_number_per_dim, cells_number_per_dim])
CG1_elem = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, CG1_elem)
T_test = TestFunction(V)
T = TrialFunction(V)
# Dirichlet boundary conditions on cube surface
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_conditions = [
(lambda x: np.isclose(x[0], -box_length/2), 'dof_bdry_1'),
(lambda x: np.isclose(x[0], box_length/2), 'dof_bdry_2'),
(lambda x: np.isclose(x[1], -box_length/2), 'dof_bdry_3'),
(lambda x: np.isclose(x[1], box_length/2), 'dof_bdry_4'),
(lambda x: np.isclose(x[2], -box_length/2), 'dof_bdry_5'),
(lambda x: np.isclose(x[2], box_length/2), 'dof_bdry_6')
]
dof_boundaries = {}
for marker, dof_name in boundary_conditions:
facets_bdry = mesh.locate_entities_boundary(msh, dim=fdim, marker=marker)
dof_boundaries[dof_name] = locate_dofs_topological(V=V, entity_dim=fdim, entities=facets_bdry)
class T_exact_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = np.exp(-self.t) * (x[0]**2 + 2*x[1]**2)
return values
class r_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = -np.exp(-self.t) * (x[0]**2 + 2*x[1]**2 + 6)
return values
T_exact = Function(V)
T_exact.interpolate(T_exact_fun(t))
r = Function(V)
r.interpolate(r_fun(t))
bc1_fun = Function(V)
bc1_fun.interpolate(T_exact_fun(t))
bc1 = dirichletbc(bc1_fun, dof_boundaries['dof_bdry_1'])
bc2_fun = Function(V)
bc2_fun.interpolate(T_exact_fun(t))
bc2 = dirichletbc(bc2_fun, dof_boundaries['dof_bdry_2'])
bc3_fun = Function(V)
bc3_fun.interpolate(T_exact_fun(t))
bc3 = dirichletbc(bc3_fun, dof_boundaries['dof_bdry_3'])
bc4_fun = Function(V)
bc4_fun.interpolate(T_exact_fun(t))
bc4 = dirichletbc(bc4_fun, dof_boundaries['dof_bdry_4'])
bc5_fun = Function(V)
bc5_fun.interpolate(T_exact_fun(t))
bc5 = dirichletbc(bc5_fun, dof_boundaries['dof_bdry_5'])
bc6_fun = Function(V)
bc6_fun.interpolate(T_exact_fun(t))
bc6 = dirichletbc(bc6_fun, dof_boundaries['dof_bdry_6'])
Dirichlet_bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
# Mass matrix
M_ufl = T_test * T * dx
M = petsc.assemble_matrix(form(M_ufl), Dirichlet_bcs)
M.assemble()
diffusion_ufl = inner(grad(T_test), grad(T)) * dx
diffusion = petsc.assemble_matrix(form(diffusion_ufl), Dirichlet_bcs)
diffusion.assemble()
reaction_ufl = - T_test * r * dx
reaction = petsc.assemble_vector(form(reaction_ufl))
reaction.assemble()
def calc_function(ts, t, T, T_dot, F):
F.set(0)
F += M * T_dot + diffusion * T - reaction
ts = PETSc.TS().create(comm=MPI.COMM_WORLD)
ts.setType('beuler')
ts.setIFunction(calc_function)
ts.setTime(t)
ts.setTimeStep(dt_init)
ts.setMaxTime(t_final)
ts.setExactFinalTime(PETSc.TS.ExactFinalTime.MATCHSTEP)
opts = PETSc.Options()
ts.setFromOptions()
T_vec = T_exact.x.petsc_vec.copy()
ts.solve(T_vec)
Now, I wonder how it is possible to extend this example, such that I can solve the nonlinear PDE \partial_t T - \nabla \cdot (T^2 +1) \nabla T - \hat{r} = 0.
Especially I am interested how I need to modify the diffusion_ufl expression, since TrialFunctions can only be used for linear PDEs. Also it seems that I would rather need a vector than a matrix in the nonlinear case. To be consistent, I also need to find a corresponding expression when defining the diffusion part in
def calc_function(ts, t, T, T_dot, F):
F.set(0)
F += M * T_dot + diffusion * T - reaction
Thanks in advance for your help.