Hello,
I wonder whether there is a more efficient way to implement Dirichlet boundary conditions using dolfinx.fem.petsc NonlinearProblem. Especially, i think that creating a NonlinearProblem within the time stepping loop is not a good way (see code below). However, since only at the call of NonlinearProblem the parameter bcs=Dirichlet_bcs is handed over, i don’t know how to update it otherwise.
Thank you in advance for your help.
Cheers
Domi
import numpy as np
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.io import XDMFFile
from dolfinx.fem import functionspace, Constant, Function, locate_dofs_topological, dirichletbc, Expression
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from basix.ufl import element
from ufl import TestFunction, dx, grad, inner
from petsc4py import PETSc
# 3d geometry parameters
box_length = 1
cells_number_per_dim = 9
# time stepping
t = 0
t_final = 10
dt_init = 1e-3
def generate_mesh(box_length, cells_number_per_dim):
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])
return msh
msh = generate_mesh(box_length, cells_number_per_dim)
CG1_elem = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, CG1_elem)
dt = Constant(msh, dt_init)
T_test = TestFunction(V)
T = Function(V)
T_ = Function(V)
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
# Initial condition
T.interpolate(T_exact_fun(t))
T_.interpolate(T_exact_fun(t))
T_exact = Function(V)
T_fun_value = T_exact_fun(t)
T_exact.interpolate(T_fun_value)
# Weak form
F = T_test * (T - T_) * dx
F += dt * inner( grad(T_test), grad(T) ) * dx
def r_fun(x):
return x[0]**2 + 2*x[1]**2 + 6 * np.exp(-t)
r_space = functionspace(msh, CG1_elem)
r = Function(r_space)
r.interpolate(r_fun)
F -= dt * T_test * r * dx
def create_Dirichlet_bc(V, T_fun, t, boundary_key):
bc_fun = Function(V)
T_fun_value = T_fun(t)
bc_fun.interpolate(T_fun_value)
bc_T = dirichletbc(bc_fun, dof_boundaries[boundary_key])
return bc_T
while (t < t_final):
bc_T_1 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_1')
bc_T_2 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_2')
bc_T_3 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_3')
bc_T_4 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_4')
bc_T_5 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_5')
bc_T_6 = create_Dirichlet_bc(V, T_exact_fun, t, 'dof_bdry_6')
Dirichlet_bcs = [bc_T_1, bc_T_2, bc_T_3, bc_T_4, bc_T_5, bc_T_6]
problem = NonlinearProblem(F, T, bcs=Dirichlet_bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual" # 'residual' or 'incremental'
t += dt.value
n, converged = solver.solve(T)
dt.value *= 1.05
