Hi all,
I am having an issue with implementing time dependent boundary conditions such that I don’t need to redefine the NonlinearProblem at each time step. Currently, as seen below in the first MWE, I am effectively creating a new problem with boundary conditions given by functions input1 & input2 (at a single point).
However, 1) this feels wildly inefficient and 2) it means that I cannot for example set the krylov solver to not precondition at each step (as a new solver is defined at every iteration). I think that this can be remedied by making a class and interpolating the function at each step to update BCs (from following other threads).
My attempt to introduce what I have mentioned above is in the second MWE but the output is just zeros. Any help would be greatly appreciated! 
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh, default_scalar_type, log
from dolfinx.io import XDMFFile
import ufl
import dolfinx.fem.petsc
import dolfinx.nls.petsc
import time
start_time = time.time()
## Create mesh
L = 1.0 # Length of the rod
nx, ny, nz = 10, 10, 100 # Number of elements in x, y, and z directions
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array(
[0.1, 0.1, L])], [nx, ny, nz], cell_type=mesh.CellType.hexahedron)
## Create function spaces
V = fem.VectorFunctionSpace(domain, ("CG", 1))
fdim = domain.topology.dim - 1
# Physical constants
E = 200e9
nu = 0.27
mu = E / (2.0 * (1 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
rho = 7990
# Time integration parameters
am = 0.2
af = 0.4
g = 0.5 + af - am
alpha_m = fem.Constant(domain, am)
alpha_f = fem.Constant(domain, af)
gamma = fem.Constant(domain, g)
beta = fem.Constant(domain, (g + 0.5)**2 / 4.0)
T = 5.0e-4
Timp = 5.0e-4
Nsteps = 100
dt = fem.Constant(domain, T / Nsteps)
## Time dep. boundary conditions
Amp=1e-8 #amplitude of perturbation
def point1(x):
return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0)), np.isclose(x[0], 0.05))
def point2(x):
return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0.1)), np.isclose(x[0], 0.05))
def input1(t):
return np.array([Amp*t, 0.0,0.0], dtype=default_scalar_type)
def input2(t):
return np.array([-1*Amp*t, 0.0,0.0], dtype=default_scalar_type)
ti=0
bc_fun1=input1(ti)
bc_fun2=input2(ti)
point1_dofs = fem.locate_dofs_geometrical(V, point1)
bc1 = fem.dirichletbc(bc_fun1, point1_dofs, V)
point2_dofs = fem.locate_dofs_geometrical(V, point2)
bc2 = fem.dirichletbc(bc_fun2, point2_dofs, V)
# Functions and fields
u_ = ufl.TestFunction(V)
u = fem.Function(V, name='Displacement')
u_old = fem.Function(V)
v_old = fem.Function(V)
a_old = fem.Function(V)
# Measures
dx = ufl.Measure("dx", domain=domain)
## GA method functions
def sigma(r):
return lmbda*ufl.nabla_div(r)*ufl.Identity(len(r)) + 2*mu*ufl.sym(ufl.grad(r))
def m(u, u_):
return rho*ufl.inner(u, u_)*dx
def k(u, u_):
return ufl.inner(sigma(u), ufl.grad(u_))*dx
## Update variables
def update_a(u, u_old, v_old, a_old, ufl_=True):
if ufl_:
dt_ = dt
beta_ = beta
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
def update_v(a, u_old, v_old, a_old, ufl_=True):
if ufl_:
return v_old + dt*((1-gamma)*a_old + gamma*a)
else:
return v_old + float(dt) * ((1 - float(gamma)) * a_old + float(gamma)*a)
def update_fields(u, u_old, v_old, a_old):
u_vec, u0_vec = u.x.array[:], u_old.x.array[:]
v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl_=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl_=False)
v_old.x.array[:] = v_vec #numpy.ndarray stored inside dolfinx.fem.function.Function
a_old.x.array[:] = a_vec
u_old.x.array[:] = u_vec
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
a_new = update_a(u, u_old, v_old, a_old, ufl_=True) #ufl.algebra.Sum
v_new = update_v(a_new, u_old, v_old, a_old, ufl_=True)
## Residual to solve
res = m(avg(a_old, a_new, alpha_m), u_) + k(avg(u_old, u, alpha_f), u_)
## Create file
xdmf_file = XDMFFile(domain.comm, "output_xdmfs/rod_straight.xdmf", "w")
xdmf_file.write_mesh(domain)
xdmf_file.write_function(u, 0)
# Write function at each time step
tt=0
for i in range(Nsteps):
tt += float(dt)
bc_fun1=input1(tt)
bc_fun2=input2(tt)
bc1 = fem.dirichletbc(bc_fun1, point1_dofs,V)
bc2 = fem.dirichletbc(bc_fun2, point2_dofs,V)
allbcs=[bc1,bc2] #set BCs
problem = fem.petsc.NonlinearProblem(res, u, bcs=allbcs) #problem setup
## Solver setup
solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
ksp.setFromOptions()
num_its, converged = solver.solve(u)
assert converged
u.x.scatter_forward()
# Update old fields with new quantities
update_fields(u, u_old, v_old, a_old)
# Write the displacement to the XDMF file
xdmf_file.write_function(u, tt)
xdmf_file.close()
if MPI.COMM_WORLD.rank==0:
print("--- %s seconds runtime ---" % (time.time() - start_time))
## rest of code ##
## Time dep. boundary conditions
Amp=1e-8 #amplitude of perturbation
def point1(x):
return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0)), np.isclose(x[0], 0.05))
def point2(x):
return np.logical_and(np.logical_and(np.isclose(x[2], 0.0), np.isclose(x[1], 0.1)), np.isclose(x[0], 0.05))
class input1class:
def __init__(self):
self.t=0.0
def eval(self,x):
return np.stack((np.full(1, Amp*self.t))), np.zeros(1),np.zeros(1)))
class input2class:
def __init__(self):
self.t=0.0
def eval(self,x):
return np.stack((np.full(1, -1*Amp*self.t))), np.zeros(1),np.zeros(1)))
input1_expr=input1class()
input1_expr.t=0.
bc_fun1=fem.Function(V)
bc_fun1.interpolate(input1_expr.eval)
input2_expr=input2class()
input2_expr.t=0.
bc_fun2=fem.Function(V)
bc_fun2.interpolate(input2_expr.eval)
point1_dofs = fem.locate_dofs_geometrical(V, point1)
point2_dofs = fem.locate_dofs_geometrical(V, point2)
bc1 = fem.dirichletbc(bc_fun1, point1_dofs)
bc2 = fem.dirichletbc(bc_fun2, point2_dofs)
allbcs=[bc1,bc2] #set BCs
## rest of code ##
for i in range(Nsteps):
tt += float(dt)
input1_expr.t=tt
bc_fun1.interpolate(input1_expr.eval)
input2_expr.t=tt
bc_fun2.interpolate(input2_expr.eval)
bc1 = fem.dirichletbc(bc_fun1, point1_dofs)
bc2 = fem.dirichletbc(bc_fun2, point2_dofs)
allbcs=[bc1,bc2]
fem.set_bc(u.vector, allbcs)
## rest of code ##