Solving a system of coupled PDEs in fenicsx

dolfinx
41

Thank you very much dear Dokken. I have updated my dolfinx version (0.7.2). That was the issue.

42

Hi dear Dooken,

Thank you for your help so far.

I just now want to study the problem in 1D and then visualize the solution using pyvista.

Following the discussions in https://fenicsproject.discourse.group/t/pyvista-with-mixed-elements-in-dolfinx/9082, I have done so:

import os

import matplotlib.pyplot as plt

import numpy as np

import ufl

from dolfinx.io import VTXWriter




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

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

# parmeters
r_0 = 0.0
r_b = 1.0  # the coordinate of right boundary

msh = mesh.create_interval(MPI.COMM_WORLD, 10, points=(0, 1))


# Define ME function space

m1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
T1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)

ME = FunctionSpace(msh, m1 * T1)

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



# Defining our problem

am = fem.Constant(msh, ScalarType(3.e-06 ))  
k = fem.Constant(msh, ScalarType(0.12))  
c = fem.Constant(msh, ScalarType(1.))     #c = fem.Constant(msh, ScalarType(1284.)) 
rho_0 = fem.Constant(msh, ScalarType(1.))      #rho_0 = fem.Constant(msh, ScalarType(500.))
R = fem.Constant(msh, ScalarType(1.))    #R = fem.Constant(msh, ScalarType(2400.))  

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


dt = 0.001  #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], r_0)  #left (x = 0)


def right(x):
    return np.isclose(x[0], r_b)  #right (x = 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_left = locate_entities_boundary(msh, msh.topology.dim - 1, left)
boundary_dofs_left = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_left)
bc_m_left = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_left, 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(30.), boundary_dofs_right, ME.sub(1))

boundary_facets_left = locate_entities_boundary(msh, msh.topology.dim - 1, left)
boundary_dofs_left = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_left)
bc_T_left = dirichletbc(PETSc.ScalarType(30.), boundary_dofs_left, ME.sub(1))


#dirichlet boundary conditions
bcs = [bc_m_right, bc_m_left, bc_T_right, bc_T_left]  




#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(lambda x: np.full(x.shape[1], 0.25, 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

#Compute the jacobian
#Jac = derivative(F, u, du)

#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 = 1e-12 #0.0 #1e-14 #0.0 #1e-12
#solver.max_it = 70 
solver.report = True
#solver.max_it = 100        #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 1


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


with io.XDMFFile(MPI.COMM_WORLD, "problem_interval_mesh.xdmf", "w") as xdmf:
     xdmf.write_mesh(msh)

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

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


import pyvista as pv
import matplotlib.pyplot as plt
pv.start_xvfb()

u0_h = u.sub(0).collapse()    #m
u1_h = u.sub(1).collapse()    #T
# Dummy function space for pyvista grid
W_dummy = FunctionSpace(msh, ("CG", 1))
#cells, types, x = plot.create_vtk_mesh(W_dummy) #dolfinx.plot.vtk_mesh
cells, types, x = dolfinx.plot.vtk_mesh(W_dummy) #dolfinx.plot.vtk_mesh
grid = pv.UnstructuredGrid(cells, types, x)

plotter = pv.Plotter()
plotter.open_gif("u0_htime.gif", fps=10)
grid.point_data["u0_h"] = u.sub(0).collapse().x.array.real

warped = grid.warp_by_scalar("u0_h", factor=1)

viridis = plt.colormaps.get_cmap("viridis").resampled(25)
sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
             position_x=0.1, position_y=0.8, width=0.8, height=0.1)

renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
                            cmap=viridis, scalar_bar_args=sargs,
                            clim=[0, max(u.x.array)])
plotter.show_axes()

The code works but I don’t see any plot.

Please let me know what I am missing.

Thank you for your time and your patience.

43

It sweems like you are trying to make a gif, but you are not adding to the gif or closing it,
using
plotter.write_frame() and plotter.close() respectively

44

There is something I am missing.
The PDE can be linear, but its solution needs not be linear. Newton method is applied to a function that should resemble the solution one. Therefore, I don’t see any reason for Newton method to converge within 1 iteration or 2 if the solution to the linear PDE is non linear (i.e. not a straight line).

What am I missing?

45

You need to think of newtons method in a multidimensjonalt setting. If the PDE is linear, then the gradient computed by the Newton method goes in the gradient direction for each degree of freedom.

1 Like
46

Thanks for your reply. I have some more reading to do understand things better. By degree of freedom, you refer to the unknowns at each nodes? I.e. for a heat equation, the degrees of freedom would be the values of the temperature at each nodes of the mesh. This is what I have in mind.

I am reading about the NonlinearSolver, and in my case (which I posted a few times in different posts in this website), I define a Jacobian, which, if I understand well, is the initial guess of the Newton’s method, and isn’t updated at each time step if the PDE is linear. For a better convergence, I think I would like to define a different Jacobian than Jac = derivative(weak_form,TempVolt,dTV). I am dealing with a thermoelectric problem, maybe I should use a simpler Jacobian, for example by solving the decoupled thermal and voltage problems and use these as initial guesses for the dofs.

47

I would suggest reading: Custom Newton solvers — FEniCSx tutorial as it goes through newtons method in detail, both for polynomial approximations and for PDEs with Dirichlet conditions

1 Like
48

Thanks dokken. What is still not very clear to me is why, or when, would one want to use a custom Newton solver as opposed to use the one from dolfinx as used in Implementation — FEniCSx tutorial ?
For example, with the one from dolfinx, I can use Dirichlet b.c., I can specify the Jacobian I want, just like with the custom Newton solvers.

49

Custom Newton solvers are good for:

  1. Debugging: If something doesn’t work or give unexpected results with the DOLFINx solver, it can be hard to pinpoint the issue if you are not a core developer.
  2. Customizability: If you want something that is not a normal Newton solver, but has some custom behavior (like a specific adaptivity step for the descent direction), it is easy to add this to a custom solver.
  3. Understanding core concepts: If one haven’t work with non-linear problems before, it is often beneficial to try to implement it once to understand how it works.
1 Like
50

Thanks a lot dokken, this is very informative. I might make a push request in the near future to mention these points in the tutorial, as I believe they are important.