Solving a system of coupled PDEs in fenicsx

[dokken]

@Mamadou

a_m \nu \Delta T term for the first PDE nonlinear and the $R \epsilon \frac{\partial m}{t} term for the second PDE.

What makes them non-linear?

@Mamadou

dt*am * inner(grad(m), grad(q)) * dx

This term is linear.

Same holds for the second term, as R and epsilon are constants.

[dokken]

If they were non-linear, you would have to have either: a power of m and T, or a combination of the two, in a single expression. What do you mean when you say they are non-linear?

[Mamadou]

OK, I get it, you are right.
These terms make both PDEs inhomogeneous.

That means I need to change my solver.

[dokken]

Since the problem is linear a nonlinear solver should converge in one iteration. This indicates that there is a bug in your code.

[Mamadou]

Thank you, dear Dokken, I’ll have a look!

[Mamadou]

Yes, the bug could be :

problem = NonlinearProblem(F, u, bcs = bcs)   

which I have to replace by

problem = LinearProblem()  

[dokken]

@Mamadou

right_m_dofs = locate_dofs_geometrical(Q, right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), right_m_dofs, Q)



top_m_dofs = locate_dofs_geometrical(Q, top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.12), top_m_dofs, Q)

right_T_dofs = locate_dofs_geometrical(V, right)
bc_T_right = dirichletbc(PETSc.ScalarType(80.), right_T_dofs, V)

top_T_dofs = locate_dofs_geometrical(V, top)
bc_T_top = dirichletbc(PETSc.ScalarType(80.), top_T_dofs, V)

#Creation of dirichlet conditions
#bcm = [bc_m_right, bc_m_top]     #bc for unknown m
#bcT = [bc_T_right, bc_T_top]     #bc for unknown T

bcs = [bc_m_right, bc_m_top, bc_T_right, bc_T_top]  #dirichlet boundary conditions

None of the boundary conditions are created correctly to work for a problem with Mixed elements. See:
https://jsdokken.com/dolfinx-tutorial/chapter3/component_bc.html
Or Solving Stokes equation with NewtonSolver - #4 by dokken

[raboynics]

Is this universally true? Because I always use the NonlinearProblem solver in FEniCSx, and it usually converges within 3 iterations, on linear problems with mixed elements. I really hope this does not indicate my problems are bugged… (I have posted several such codes in the past here).

[dokken]

Not sure what your question relates to.
If you solve problems in a mixed space, the Dirichlet condition has to relate to either the full space, or a sub space of the mixed space.

[raboynics]

I am talking about the convergence of the Newton solver. Does it really has to converge in 1 iteration if the problem is well posed and linear? Or is 3 iterations acceptable, as I thought?
For example there: Krylov solver's option max_it and a few questions. Newton method converges in 4 steps, I am dealing with a linear problem (albeit with “tensors” in a 2D problem, with mixed elements).

[dokken]

@raboynics

am talking about the convergence of the Newton solver. Does it really has to converge in 1 iteration if the problem is well posed and linear? Or is 3 iterations acceptable, as I thought?

If the problem is truely linear, it should converge in 1 iteration (one eval + one gradient correction).

I guess with the incremental Newton solver it might take 2 iterations, as it has to build the residual.

If you chance the convergence criterion to residual, and use direct solvers for the linearized problem, it should converge in a single iteration.

If you find this confusing, I would read: Custom Newton solvers — FEniCSx tutorial
and think about what this means for a single variable problem:
Content - Newton's method.

[Mamadou]

Following this tutorial Component-wise Dirichlet BC — FEniCSx tutorial, I implement the right boundary condition m(x=1, y, t) = 0.12 :

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (lx, Ly)), n=(nx, ny),
                            cell_type=mesh.CellType.triangle)


m1 = ufl.FiniteElement("CG", msh.ufl_cell(), 2)
T1 = ufl.FiniteElement("CG", msh.ufl_cell(), 2)
ME = FunctionSpace(msh, m1 * T1)

Q = FunctionSpace(msh, m1)
V = FunctionSpace(msh, T1)

def right(x):
    return np.isclose(x[0], lx)  #right (x = 1)

boundary_facets = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_x = locate_dofs_topological(Q, msh.topology.dim - 1, boundary_facets)
bcx = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_x, Q)

Is this correct now?

[dokken]

@Mamadou

boundary_facets = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_x = locate_dofs_topological(Q, msh.topology.dim - 1, boundary_facets)
bcx = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_x, Q)

No, you need to apply the Dirichlet BC to the subspace (ME.sub(0))

[Mamadou]

Thank you, dear Dokken, for your quick reply.
That means in the implementation, I will have:

boundary_facets = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_x = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets)
bcx = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_x, ME.sub(0))

and

boundary_facets = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_x = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets)
bcx = dirichletbc(PETSc.ScalarType(80.), boundary_dofs_x, ME.sub(1))

for m(x=1, y, t) = 0.12 and T(x=1, y, t) = 80.

[dokken]

Yes, then you are restricting the correct dofs

[Mamadou]

I use the NonlinearProblem solver in Fenicsx to solve my linear problem with mixed elements as Raboynics has done here: Krylov solver's option max_it and a few questions . I also change the convergence criterion.

With solver.convergence_criterion = "residual" I got:

Step 1: num iterations: 1
Step 2: num iterations: 0
Step 3: num iterations: 0
Step 4: num iterations: 0
Step 5: num iterations: 0
Step 6: num iterations: 0
Step 7: num iterations: 0
Step 8: num iterations: 0
Step 9: num iterations: 0
Step 10: num iterations: 0
Step 11: num iterations: 0
Step 12: num iterations: 0
Step 13: num iterations: 0
Step 14: num iterations: 0
Step 15: num iterations: 0
Step 16: num iterations: 0
Step 17: num iterations: 0
Step 18: num iterations: 0
Step 19: num iterations: 0
Step 20: num iterations: 0
Step 21: num iterations: 0
Step 22: num iterations: 0
Step 23: num iterations: 0
Step 24: num iterations: 0
Step 25: num iterations: 0
Step 26: num iterations: 0
Step 27: num iterations: 0
Step 28: num iterations: 0
Step 29: num iterations: 0
Step 30: num iterations: 0
Step 31: num iterations: 0
Step 32: num iterations: 0
Step 33: num iterations: 0
Step 34: num iterations: 0
Step 35: num iterations: 0
Step 36: num iterations: 0
Step 37: num iterations: 0
Step 38: num iterations: 0
Step 39: num iterations: 0
Step 40: num iterations: 0
Step 41: num iterations: 0
Step 42: num iterations: 0
Step 43: num iterations: 0
Step 44: num iterations: 0
Step 45: num iterations: 0
Step 46: num iterations: 0
Step 47: num iterations: 0
Step 48: num iterations: 0
Step 49: num iterations: 0
Step 50: num iterations: 0
2023-12-17 13:02:16.338 (   7.638s) [          BA8740]   vtkTextureObject.cxx:1025   ERR| vtkTextureObject (0x55e487efdac0): Attempt to use a texture buffer exceeding your hardware's limits. This can happen when trying to color by cell data with a large dataset. Hardware limit is 65536 values while 73728 was requested.
ERROR:root:Attempt to use a texture buffer exceeding your hardware's limits. This can happen when trying to color by cell data with a large dataset. Hardware limit is 65536 values while 73728 was requested. 

and for solver.convergence_criterion = "incremental'' I got:

DOLFINx version: 0.6.0 based on GIT commit: 30920d400ccb64d190de13509c43b1755d2606cc of https://github.com/FEniCS/dolfinx/
Step 1: num iterations: 2
Step 2: num iterations: 2
Step 3: num iterations: 2
Step 4: num iterations: 2
Step 5: num iterations: 2
Step 6: num iterations: 2
Step 7: num iterations: 2
Step 8: num iterations: 2
Step 9: num iterations: 2
Step 10: num iterations: 2
Step 11: num iterations: 2
Step 12: num iterations: 2
Step 13: num iterations: 2
Step 14: num iterations: 2
Step 15: num iterations: 2
Step 16: num iterations: 2
Step 17: num iterations: 2
Step 18: num iterations: 2
Step 19: num iterations: 2
Step 20: num iterations: 2
Step 21: num iterations: 2
Step 22: num iterations: 2
Step 23: num iterations: 2
Step 24: num iterations: 2
Step 25: num iterations: 2
Step 26: num iterations: 2
Step 27: num iterations: 2
Step 28: num iterations: 2
Step 29: num iterations: 2
Step 30: num iterations: 2
Step 31: num iterations: 2
Step 32: num iterations: 2
Step 33: num iterations: 2
Step 34: num iterations: 2
Step 35: num iterations: 2
Step 36: num iterations: 2
Step 37: num iterations: 2
Step 38: num iterations: 2
Step 39: num iterations: 2
Step 40: num iterations: 2
Step 41: num iterations: 2
Step 42: num iterations: 2
Step 43: num iterations: 2
Step 44: num iterations: 2
Step 45: num iterations: 2
Step 46: num iterations: 2
Step 47: num iterations: 2
Step 48: num iterations: 2
Step 49: num iterations: 2
Step 50: num iterations: 2
2023-12-17 13:19:05.367 (  46.251s) [        DA8FA740]   vtkTextureObject.cxx:1025   ERR| vtkTextureObject (0x555ad4b60cf0): Attempt to use a texture buffer exceeding your hardware's limits. This can happen when trying to color by cell data with a large dataset. Hardware limit is 65536 values while 73728 was requested.
ERROR:root:Attempt to use a texture buffer exceeding your hardware's limits. This can happen when trying to color by cell data with a large dataset. Hardware limit is 65536 values while 73728 was requested.

Note also that the convergence is faster with ‘‘residual’’ than with ‘‘incremental’’.

[dokken]

@Mamadou

Step 1: num iterations: 1
Step 2: num iterations: 0
Step 3: num iterations: 0

So the first step indicates that the problem is linear. The subsequent steps indicates that you are not updating anything within your loop, so the initial guess is the solution to your problem.

@Mamadou

Note also that the convergence is faster with ‘‘residual’’ than with ‘‘incremental’’.

Incremental requires as least 2 steps to converge, as it is part of the convergence criteria.

[Mamadou]

Thanks for the useful information, dear Dokken.
I got the same result when modifying the solver.convergence_criterion = "residual" in the Cahn-Hilliard equation — DOLFINx 0.5.1 documentation :

Step 1: num iterations: 1
Step 2: num iterations: 0
Step 3: num iterations: 0
Step 4: num iterations: 0
Step 5: num iterations: 0
Step 6: num iterations: 0
Step 7: num iterations: 0
Step 8: num iterations: 0
Step 9: num iterations: 0
Step 10: num iterations: 0
Step 11: num iterations: 0
Step 12: num iterations: 0
Step 13: num iterations: 0
Step 14: num iterations: 0
Step 15: num iterations: 0
Step 16: num iterations: 0
Step 17: num iterations: 0
Step 18: num iterations: 0
Step 19: num iterations: 0
Step 20: num iterations: 0
Step 21: num iterations: 0
Step 22: num iterations: 0
Step 23: num iterations: 0
Step 24: num iterations: 0
Step 25: num iterations: 0
Step 26: num iterations: 0
Step 27: num iterations: 0
Step 28: num iterations: 0
Step 29: num iterations: 0
Step 30: num iterations: 0
Step 31: num iterations: 0
Step 32: num iterations: 0
Step 33: num iterations: 0
Step 34: num iterations: 0
Step 35: num iterations: 0
Step 36: num iterations: 0
Step 37: num iterations: 0
Step 38: num iterations: 0
Step 39: num iterations: 0
Step 40: num iterations: 0
Step 41: num iterations: 0
Step 42: num iterations: 0
Step 43: num iterations: 0
Step 44: num iterations: 0
Step 45: num iterations: 0
Step 46: num iterations: 0
Step 47: num iterations: 0
Step 48: num iterations: 0
Step 49: num iterations: 0
Step 50: num iterations: 0

This means that the Cahn-Hilliard equations considered in this implementation are two coupled second order equations and linear.

[Mamadou]

Can you help me visualise the plots of my problem using Paraview?

import os

import numpy as np

import ufl

from dolfinx.io import VTXWriter

#from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
                 #TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
                 #inv, split, derivative)



from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, set_bc)

#from dolfinx.fem import (Constant, dirichletbc,
                         #extract_function_spaces, form,
                        # locate_dofs_topological)
from dolfinx import log, plot
#from dolfinx.plot import vtk_mesh
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, TrialFunction, derivative
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType

from mpi4py import MPI
from petsc4py import PETSc

try:
    import pyvista as pv
    import pyvistaqt as pvqt
    have_pyvista = True
    if pv.OFF_SCREEN:
        pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
    print("pyvista and pyvistaqt are required to visualise the solution")
    have_pyvista = False

# Save all logging to file
log.set_output_file("log.txt")

import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

lx = 1.0  #length of x-axis
Ly = 1.0  #length of y-axis
#nx, ny = 8, 8, 50, 50, 96, 96 
nx = 96
ny = 96

#t = 0
#T = 10

# A rectangle mesh with 8 cells edges in each direction is created,
# and on this mesh a
# {py:class}`FunctionSpace<dolfinx.fem.FunctionSpace>` `ME` is built
# using a pair of linear Lagrange elements.

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (lx, Ly)), n=(nx, ny),
                            cell_type=mesh.CellType.triangle)

# Define ME function space
#el = ufl.FiniteElement("CG", mesh.ufl_cell(), 2)
#mel = MixedElement([el, el])
#ME = FunctionSpace(mesh, mel)

#msh = create_rectangle(MPI.COMM_WORLD, 96, 96, CellType.triangle)
#m1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
#T1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
m1 = ufl.FiniteElement("CG", msh.ufl_cell(), 2)
T1 = ufl.FiniteElement("CG", msh.ufl_cell(), 2)
#el = ufl.FiniteElement("CG", mesh.ufl_cell(), 2)
#m1 = ufl.FiniteElement("DG", msh.ufl_cell(), 4)
#T1 = ufl.FiniteElement("DG", msh.ufl_cell(), 4)
ME = FunctionSpace(msh, m1 * T1)

Q = FunctionSpace(msh, m1) #ME.sub(0)
V = FunctionSpace(msh, T1) #ME.sub(1)



# Defining our problem
#p = fem.Constant(domain, ScalarType(1)) # SIMP penalty exponent
am = fem.Constant(msh, ScalarType(3.e-06 ))  
k = fem.Constant(msh, ScalarType(0.12))  
c = fem.Constant(msh, ScalarType(1284.))  
rho_0 = fem.Constant(msh, ScalarType(500.)) 
R = fem.Constant(msh, ScalarType(2400.))   

epsilon = fem.Constant(msh, ScalarType(0.1))  
nu = fem.Constant(msh, ScalarType(1.e-02))        


dt = 5.0e-06  # 0.001  #0.02 #5.0e-06  # time step


#We create four python functions for determining the facets to apply boundary conditions to

def left(x):
    return np.isclose(x[0], 0)  #left (x = 0)


def right(x):
    return np.isclose(x[0], lx)  #right (x = 1)

def bottom(x):
    return np.isclose(x[1], 0)  #bottom (y = 0)


def top(x):
    return np.isclose(x[1], Ly)  #top (y = 1)



boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_right, ME.sub(0))

boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_top, ME.sub(0))


boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_right)
bc_T_right = dirichletbc(PETSc.ScalarType(80.), boundary_dofs_right, ME.sub(1))

boundary_facets_top = locate_entities_boundary(msh, msh.topology.dim - 1, top)
boundary_dofs_top = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_top)
bc_T_top = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_top, ME.sub(1))

bcs = [bc_m_right, bc_m_top, bc_T_right, bc_T_top]  #dirichlet boundary conditions






#Test Functions (q, v)
q, v = ufl.TestFunctions(ME)

u = Function(ME)  # current solution
u_n = Function(ME)  # solution from previous converged step
#du = TrialFunction(ME)

# Split mixed functions
m, T = ufl.split(u)
m_n, T_n = ufl.split(u_n)

# Weak statement of the equations
m_n = Function(Q)
m_n.name = "m_n"
T_n = Function(V)
T_n.name = "T_n"

# The initial conditions are interpolated into a finite element space:


# Zero u
u.x.array[:] = 0.0

# Interpolate initial condition
#u.sub(0).interpolate(PETSc.ScalarType(0.30))  
u.sub(0).interpolate(lambda x: np.full(x.shape[1], 0.50, dtype=np.float64)) #Initial condition for m
u.sub(1).interpolate(lambda x: np.full(x.shape[1], 10., dtype=np.float64)) #Initial condition for T
u.x.scatter_forward()



#Variational formulation of the system
F0 = inner(m, q) * dx - inner(m_n, q) * dx + dt*am * inner(grad(m), grad(q)) * dx + dt*am*nu * inner(grad(T), grad(q)) * dx
F1 = c*rho_0* inner(T, v) * dx - c*rho_0* inner(T_n, v) * dx + dt*k * inner(grad(T), grad(v)) * dx - R*epsilon * inner(m, v) * dx + R*epsilon * inner(m_n, v) * dx
F = F0 + F1



#Solve variational problem
# Create nonlinear problem and Newton solver
#problem = NonlinearProblem(F, u, bcs = bcs)     #Replace
problem = NonlinearProblem(F, u, bcs = bcs)  
#problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)  
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"      #residual    #incremental
solver.rtol = 0.0 #1e-14 #0.0 #1e-12
solver.report = True
solver.error_on_nonconvergence = False
solver.max_it = 1


# 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_max_it"]= 10000
ksp.setFromOptions()



info = solver.solve(u)
print(info)


with dolfinx.io.VTXWriter(msh.comm, "m.bp", [u.sub(0).collapse()], engine="BP4") as vtx:
    vtx.write(0.0)

with dolfinx.io.VTXWriter(msh.comm, "T.bp", [u.sub(1).collapse()], engine="BP4") as vtx:
    vtx.write(0.0)

I got the following error (I am using paraview 5.9.0):

Traceback (most recent call last):
  File "/home/mamadou/Downloads/fawzoulazim2023_december2023_10.py", line 215, in <module>
    with dolfinx.io.VTXWriter(msh.comm, "m.bp", [u.sub(0).collapse()], engine="BP4") as vtx:
TypeError: VTXWriter.__init__() got an unexpected keyword argument 'engine'

[dokken]

You seem to have an old installation of DOLFINx, please remove the engine keyword.