Thank you Dr. Dokken for your reply.
I am trying to approximate the following system of PDEs (with m and T are the unknowns functions).
\dfrac{\partial {m}}{\partial t} \ = \ a_m \Delta m + a_m \nu \Delta T \ , \hspace{1cm} in \ \Omega \times [0, tf] \hspace{4cm} \textbf{(1)}
c \rho_{0} \dfrac{\partial {T}}{\partial t} \ = k \Delta T + R \epsilon \dfrac{\partial {m}}{\partial t} \ , \hspace{1cm} in \ \Omega \times [0, tf] \hspace{4.2cm} \textbf{(2)}
with initial conditions m(x, y, 0) = 0.30 and T(x, y, 0) = 10.
The boundary boundary conditions are:
m(x=1, y, t) = 0.12, m(x, y=1, t), T(x=1, y, t) = 60 and T(x, y=1, t) = 60.
To do so, I am following this tutorial Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial, using a mixed space element as in Solving Stokes equation with NewtonSolver - #4 by dokken and https://fenicsproject.org/olddocs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html.
Here is the implementation of the mesh, mixed element, boundary conditions and initial conditions, variational formulation. Using a NonProblem solver, it converges in two iterations.
import os
import matplotlib.pyplot as plt
import numpy as np
import ufl
from dolfinx.io import VTXWriter
import math
from math import sqrt
# For plotting
import matplotlib.pyplot as plt
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, dot
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType
import math
from math import sqrt
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
r_b = 1 # the coordinate of right boundary
# Define mesh #nums = [8, 16, 32,64,128,255] number of cells
nx, ny = 15, 15
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([r_0, r_0]), np.array([r_b, r_b])],
[nx, ny], mesh.CellType.quadrilateral)
#
###################################################
#### Parameters : time and load and material ######
###################################################
#
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.))
alpha = k/(c*rho_0)
beta = R/(c*rho_0)
epsilon = fem.Constant(msh, ScalarType(0.1))
nu = fem.Constant(msh, ScalarType(1.e-02))
## Time parametrization
t = 0 # Start time
Tf = 10. # End time
num_steps = 100 # Number of time steps
dt = (Tf-t)/num_steps # Time step size
#
###################################################
########## Mixed Space of resolution ##############
###################################################
#
# Define ME function space
# Lagrange" is type of the finite element, and 1 is the degree
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)
# Create the initial and solution functions of space
u = Function(ME) # current solution
u_n = Function(ME) # solution from previous converged step (i.e., initial solution)
#du = TrialFunction(ME)
#We create four python functions for determining the facets to apply boundary conditions to
def right(x):
return np.isclose(x[0], r_b) #right (x = 0.1)
def top(x):
return np.isclose(x[1], r_b) #top (y = 0.1)
#
###################################################
######## Dirichlet boundary conditions #############
###################################################
#
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(60.), 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(60.), boundary_dofs_top, ME.sub(1))
bcs = [bc_m_right, bc_m_top, bc_T_right, bc_T_top] #dirichlet boundary conditions
#
#
###################################################
############## Initial Values #####################
###################################################
#
Mi = 0.30 #Initial moisture content
#
M_n_, M_n_to_ME = ME.sub(0).collapse()
FM_n_ = Function(M_n_)
with FM_n_.vector.localForm() as initial_local:
initial_local.set(ScalarType(Mi))
# # Re assign in u_n
u_n.x.array[M_n_to_ME] = FM_n_.x.array
u_n.x.scatter_forward()
Ti = 10. #Initial temperature content
T_n_, T_n_to_ME = ME.sub(1).collapse()
FT_n_ = Function(T_n_)
with FT_n_.vector.localForm() as initial_local:
initial_local.set(ScalarType(Ti))
# # Re assign in u_n
u_n.x.array[T_n_to_ME] = FT_n_.x.array
u_n.x.scatter_forward()
#
# Identify the unknowns from the function
# Split mixed functions
m, T = ufl.split(u)
m_n, T_n = ufl.split(u_n)
# Set up the test functions (q, v)
q, v = ufl.TestFunctions(ME)
#
###################################################
############## Variationnal form ##################
###################################################
#
#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
F1 = inner(T, v) * dx - inner(T_n, v) * dx + dt*alpha * inner(grad(T), grad(v)) * dx - beta*epsilon * inner(m, v) * dx + beta*epsilon * inner(m_n, v) * dx
F = F0 + F1
# Non linear problem definition
du = TrialFunction(ME)
#Compute the jacobian
Jac = derivative(F, u, du)
problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)
#Solve variational problem
# Create nonlinear problem and Newton solver
#problem = NonlinearProblem(F, u, bcs = bcs) #Replace
#problem = NonlinearProblem(F, u, bcs = bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental" #residual #incremental
# Absolute tolerance
solver.atol = 5e-10
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 = 2
info = solver.solve(u)
print(info)
The next step is to create a loop over time. I built three xdmf files with xdmf1 and xdmf2 containing the solutions m and T.
#
###################################################
################ Pre-processing ###################
###################################################
#
# Create an output xdmf file to store the values
# xdmf = XDMFFile(mesh.comm, "Column.xdmf", "w")
# xdmf.write_mesh(mesh)
# We prepare output files for the velocity and pressure data, and write the mesh and initial conditions to file
# Create an output xdmf file to store the values
file = XDMFFile(msh.comm, "problem_rectangle_mesh_4th_march2024_alhaatim.xdmf", "w")
file.write_mesh(msh)
#with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", [u_n.sub(0).collapse()], "w") as xdmf1:
with io.XDMFFile(MPI.COMM_WORLD, "moisture_4th_March2024.xdmf", "w") as xdmf1:
xdmf1.write_mesh(msh)
with io.XDMFFile(MPI.COMM_WORLD, "temperature_4th_March2024.xdmf", "w") as xdmf2:
xdmf2.write_mesh(msh)
#
###################################################
################ Processing #######################
###################################################
# time steps to evaluate the pressure in space:
#n0, n1, n2 = 200,400,800
#
t = 0
#L2_p = np.zeros(num_steps, dtype=PETSc.ScalarType)
for n in range(num_steps):
t += dt
try:
num_its, converged = solver.solve(u)
except:
if MPI.COMM_WORLD.rank == 0:
print("*************")
print("Solver failed")
print("*************")
pass
u.x.scatter_forward()
# Update Value
u_n.x.array[:] = u.x.array
u_n.x.scatter_forward()
__m, __T = u.split()
#
# Export the results
__m.name = "Moisture"
__T.name = "Temperature"
xdmf1.write_function(__m,t)
xdmf2.write_function(__T,t)
#
#
file.close()
#
xdmf1.close()
xdmf2.close()
The above doesn’t return an error, but when I try to use Paraview to visualise my results, it doesn’t work (i.e. empty file with no data).
Thank you for your time and your patience.