Hi dear fenicsx users.
I am trying to solve the following system of PDEs (in 2D) using fenicsx (dolfinx version 0.6.0).
\dfrac{\partial {m}}{\partial t} \ = \ a_m \ \Delta m + a_m \nu \ \Delta T
c \rho_{0} \dfrac{\partial {T}}{\partial t} \ = k \ \Delta T + R \epsilon \ \dfrac{\partial {m}}{\partial t}
with initial conditions
m(x, 0) = 0.30 and T(x, 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.
I had a close look on the Cahn-Hilliard and Navier-Stockes equations and I have tried to implement this problem in fenicsx but I am encountering some issues with my code.
My code is as follow:
import os
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)
#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
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
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 = 50, 50
nx = 8
ny = 8
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)
#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)
ME = FunctionSpace(msh, m1 * T1)
Q = FunctionSpace(msh, m1)
V = FunctionSpace(msh, T1)
m = ufl.TrialFunction(Q)
q = ufl.TestFunction(Q)
T = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# 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))
m_ini = fem.Constant(msh, ScalarType(0.5))
m_lbou = fem.Constant(msh, ScalarType(0.12))
m_rbou = fem.Constant(msh, ScalarType(0.12))
T_ini = fem.Constant(msh, ScalarType(10.))
T_lbou = fem.Constant(msh, ScalarType(80.))
T_rbou = fem.Constant(msh, ScalarType(80.))
dt = 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) #right (y = 1)
right_m_dofs = locate_dofs_geometrical(m1, right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), right_m_dofs, m1)
top_m_dofs = locate_dofs_geometrical(m1, top)
bc_m_top = dirichletbc(PETSc.ScalarType(0.12), top_m_dofs, m1)
right_T_dofs = locate_dofs_geometrical(T1, right)
bc_T_right = dirichletbc(PETSc.ScalarType(60.), right_T_dofs, T1)
top_T_dofs = locate_dofs_geometrical(T1, top)
bc_T_top = dirichletbc(PETSc.ScalarType(60.), top_T_dofs, T1)
bcm = [bc_m_right, bc_m_top]
bcT = [bc_T_right, bc_T_top]
# +
u = Function(ME) # current solution
u_n = Function(ME) # solution from previous converged step
# 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"
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
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs = [bcm, bcT])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
# 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"
ksp.setFromOptions()
# -
# -
# The setting of `convergence_criterion` to `"incremental"` specifies
# that the Newton solver should compute a norm of the solution increment
# to check for convergence (the other possibility is to use
# `"residual"`, or to provide a user-defined check). The tolerance for
# convergence is specified by `rtol`.
#
# To run the solver and save the output to a VTK file for later
# visualization, the solver is advanced in time from $t_{n}$ to
# $t_{n+1}$ until a terminal time $T$ is reached:
# +
# Output file
file = XDMFFile(MPI.COMM_WORLD, "demo_ch/output.xdmf", "w")
file.write_mesh(msh)
from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "poiseuille_m.bp", m0, engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_T.bp", T0, engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
When running the code, I got the following error:
DOLFINx version: 0.6.0 based on GIT commit: 30920d400ccb64d190de13509c43b1755d2606cc of https://github.com/FEniCS/dolfinx/
Traceback (most recent call last):
File "/home/mamadou/Downloads/fawzoulazim2023_december2023_2.py", line 115, in <module>
right_m_dofs = locate_dofs_geometrical(m1, right)
File "/home/mamadou/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 56, in locate_dofs_geometrical
raise TypeError
TypeError
Thanks for your time and attention.