Ok I have a reasonably minimal example, below.
Key facts:
- I can replicate the error only when I include the non-linear advection terms (their lines are noted in comments)
- The original behaviour still persists: I can run the code just fine with or without the non-linear terms for some given number of runs. But running them in a big batch causes the error. Specifically, I have two batches of around 200 each running currently. Adding a batch over ~60 now causes the error.
- It uses the nonlinear_assembly code that is part of dolfinx_mpc
Code:
from mpi4py import MPI
from dolfinx.fem import Function, functionspace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import (create_interval, locate_entities_boundary, meshtags)
from basix.ufl import (element, mixed_element)
import numpy as np
import math
from ufl import (TestFunctions, dx, inner, split, derivative, ln)
import dolfinx_mpc
from nonlinear_assembly import *
from pathlib import Path
import sys
import time
import os
comm = MPI.COMM_WORLD
########################################
########## ARGUMENTS + PATH ############
########################################
id_str = sys.argv[1]
dirStr = "data/" +id_str + "/"
os.makedirs(os.path.dirname(dirStr), exist_ok = True)
########################################
############## OPEN FILES ##############
########################################
if comm.rank == 0:
outfileStr = dirStr + "out.dat"
outfile = open(outfileStr, "w")
########################################
############### CONSTANTS ##############
########################################
comm = MPI.COMM_WORLD
mesh_size = 1000
dt = 0.01
T = 10000
loop_count = 0
output_step = int(1 / dt) # number of timesteps per output interval
# domain
x_left = 0.0
x_right = 100.0
########################################
########## INITIAL CONDITIONS ##########
########################################
def f1_init(x):
values = np.zeros((1, x.shape[1]))
for i,v in enumerate(x[0]):
values[0][i] = 0.5 + 0.1 * np.sin(x[0][i] * 2 * np.pi / 100)
return values
def f2_init(x):
values = np.zeros((1, x.shape[1]))
for i,v in enumerate(x[0]):
values[0][i] = 0.1 * np.sin(x[0][i] * 2 * np.pi / 100)
return values
########################################
########## BOUNDARY FUNCTIONS ##########
########################################
def PeriodicBoundary(x):
return np.isclose(x[0], x_right)
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = x_right - x[0]
return out_x
########################################
######## FUNCTION SPACE OBJECTS ########
########################################
# Create mesh and define function space
msh = create_interval(comm, points=(x_left, x_right), nx=mesh_size, ghost_mode=dolfinx.cpp.mesh.GhostMode.shared_facet)
P1 = element("Lagrange", msh.basix_cell(),1)
V = functionspace(msh, mixed_element([P1, P1]))
V_0, _ = V.sub(0).collapse()
num_points_local = V_0.dofmap.index_map.size_local * V_0.dofmap.index_map_bs
bcs = [] # no other boundary conditions
facets = locate_entities_boundary(msh, msh.topology.dim - 1, PeriodicBoundary)
arg_sort = np.argsort(facets)
mt = meshtags(msh, msh.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype = np.int32))
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, 1) # constraint on f1
mpc.create_periodic_constraint_topological(V.sub(1), mt, 2, periodic_relation, bcs, 1) # constraint on psi
mpc.finalize()
u = Function(V) # current solution
u_prev = Function(V) # solution from previous converged step
f1, f2 = split(u)
f1_prev, f2_prev = split(u_prev)
xvals = sorted(msh.geometry.x[:,0])
dx_mesh = xvals[1] - xvals[0]
# all non-autonomous variables must be declared as constant on the mesh to avoid recompilation
t = dolfinx.fem.Constant(msh, 0.0)
del_t = dolfinx.fem.Constant(msh, dt)
########################################
########## INITIAL CONDITIONS ##########
########################################
u.sub(0).interpolate(f1_init)
u.sub(1).interpolate(f2_init)
u.x.scatter_forward()
########################################
########### WEAK FORMULATION ###########
########################################
theta = 0.5
qr, qa = TestFunctions(V)
f1_mid = (1.0 - theta) * f1_prev + theta * f1
f2_mid = (1.0 - theta) * f2_prev + theta * f2
nonlin_f2_mid = (1.0 - theta) * f1_prev * (1 - f1_prev) + theta * f1 * (1 - f1)
nonlin_f1_mid = (1.0 - theta) * f2_prev * (1 - f1_prev) + theta * f2 * (1 - f1)
F = \
( \
# time derivative
inner (f1 - f1_prev, qr) / dt \
+ inner (f2 - f2_prev, qa) / dt \
# diffusion terms
+ inner ((f1_mid).dx(0), (qr).dx(0)) \
+ inner ((f2_mid).dx(0), (qa).dx(0)) \
# advection terms
+ 10.0 * inner ((nonlin_f1_mid).dx(0), qr) \
+ 10.0 * inner ((nonlin_f2_mid).dx(0), qa) \
#
) * dx
J = derivative (F, u)
jit_options = {"cffi_extra_compile_args": ["-Ofast", "-march=native"], "cffi_libraries": ["m"]}
F_compiled = dolfinx.fem.form(F, jit_options=jit_options)
J_compiled = dolfinx.fem.form(J, jit_options=jit_options)
problem = NonlinearMPCProblem(F_compiled, u, mpc, bcs = bcs, J = J_compiled, jit_options = jit_options)
solver = NewtonSolverMPC(msh.comm, problem, mpc)
solver.rtol = 1e-7
u_prev.x.array[:] = u.x.array
########################################
########## INTEGRATE/TIME LOOP #########
########################################
while (t.value < T):
t.value += del_t.value
loop_count += 1
r = solver.solve(u)
if comm.rank == 0:
if (loop_count % output_step == 0):
print('on time ' + str(t.value), file = outfile, flush = True)
u_prev.x.array[:] = u.x.array
########################################
################ TIDY UP ###############
########################################
if comm.rank == 0:
outfile.close()