Unexpected Error/Deviation from Dolfinx simulation

General
1

Hey Guys. I am trying to simulate Cahn-hilliard type equation(with couple with flow equation). I write the code in dolfinx and I think it is working, however, the simulation would sometimes cause local explosion that leads to huge deviation.

So this is what it looks like(nice phase separation and goes stable dynamic as time goes on):

unnamed
unnamed1044×967 129 KB

However, after a while the tolerances will become larger and larger and local explosion occurs and spread to other place, which is clearly not consistent with the equation.

unnamed
unnamed1044×967 141 KB

unnamed
unnamed1044×967 141 KB

So I am wondering is there something about the solver or setting or the coding that I can adjust to avoid this issue. I tried to change dx, dt and other parameters for my simulation. It somewhat helped but not fix the issue totally.

I know I ask kind of general question and you might need more specific information on how I write the code, which type of solver I used, and the equation. I want to keep the original question concise so I didn’t include them here but if you need any other information plz feel free to left a comment and I will upload it.

Any help and advice would be greatly appreciated!

2

Looks like your boundary conditions in the advection discretisation may be inconsistent with the original problem; or the problem is ilposed. To give more help you’d have to give a minimal working example.

3

Hi Nate. Thanks for your comment. Here is a minimal working example code:

import numpy as np
import ufl
from dolfinx import log, plot
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, CellType
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div
from mpi4py import MPI
from petsc4py import PETSc
import pyvista as pv


log.set_output_file("log.txt")


kappa = 0.02
eta = 1.5
dt = 3.0e-02# time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson




msh = create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (30.0, 30.0)), n=(80, 80),
cell_type=CellType.triangle,)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
Mixed_Element = ufl.MixedElement([P1, P1, P1, P2])
ME = FunctionSpace(msh, Mixed_Element)
q, v, m, n = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step


# Split mixed functions
c, mu, p, a = ufl.split(u)
c0, mu0, p0, a0 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0


# Interpolate initial condition
u.sub(0).interpolate(lambda x: 0.50 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()


c = ufl.variable(c)
y = 0.1
f = y * (-2 * (c-1/2)**2 + 4 * (c-1/2)**4)
dfdc = ufl.diff(f, c)


mu_mid = (1.0 - theta) * mu0 + theta * mu


# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - kappa * inner(grad(c), grad(v)) * dx
F3 = -eta * inner(grad(a), grad(n)) * dx - inner(grad(p), n) * dx - inner(c, inner(grad(mu), n)) * dx
F4 = inner(div(a), m) * dx
F = F0 + F1 + F3 + F4


problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e+3
solver.atol = 1e+3
solver.max_it = 1000

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_monitor"] = None
ksp.setFromOptions()

# Output file
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)

# Step in time
t = 0.0
T = 400 * dt

# vector
V0, dofs = ME.sub(0).collapse()

# Prepare viewer for plotting the solution during the computation
topology, cell_types, x = plot.create_vtk_mesh(V0)
grid = pv.UnstructuredGrid(topology, cell_types, x)

# Set output data
grid.point_data["c"] = u.x.array[dofs].real
grid.set_active_scalars("c")
c = u.sub(0)
u0.x.array[:] = u.x.array
time_data = []


while (t < T):
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    file.write_function(c, t)
    time_data.append(u.x.array[dofs].real.copy())




# Create a plotter object
plotter = pv.Plotter(off_screen=True)

# Add the initial mesh to the plotter
grid.point_data["c"] = time_data[0]
plotter.add_mesh(grid, scalars="c", clim=[0, 1])


# Set the view to 2D
plotter.view_xy()


# Open a mp4 file to store the animation
plotter.open_movie("Model-H.mp4")


# Iterate over each time step
for i in range(len(time_data)):
# Update the scalar data
grid.point_data["c"] = time_data[i]
plotter.add_text(f"Time: {((i + 1)*dt):.4f}", name='time-label')
plotter.write_frame()


plotter.close()