Hi dear all,
Thanks in advance.
I am trying to solve the time-dependent non-linear partial differential equations using complex Fenicsx.
The system reads:
\begin{equation}
\begin{cases}
\rho_{\alpha} \Bigl(C_{p \alpha} + C_{pv} {W} \Bigr) \Bigl( \dfrac{\partial {T_\alpha}}{\partial t} + V_x \dfrac{\partial {T_\alpha}}{\partial x} \Bigr)\ = \ \dfrac{\alpha \nu \ ({T_g} - {T_\alpha})}{\varepsilon} \vspace{0.5cm} \
\rho_g \Bigl(C_{pg} + C_{pw} {M} \Bigr) \ \dfrac{\partial {T_g}}{\partial t} \ = \ \dfrac{\nu \alpha}{1 - \varepsilon} ({T_\alpha} - {T_g}) \
\hspace{0.7cm} - \Bigl(H_v + C_{pv} \Bigl({T_\alpha} - {T_g} \Bigr) \Bigr) \ \dfrac{\rho_{\alpha} \varepsilon}{1 - \varepsilon} \ V_x \dfrac{\partial {W}}{\partial x} \vspace{0.5cm} \
\rho_{\alpha} \varepsilon \Bigl( \dfrac{\partial {W}}{\partial t} + V_x \dfrac{\partial {W}}{\partial x} \Bigr) \ = \ \Bigl(\varepsilon - 1 \Bigr) \rho_g \dfrac{\partial {M}}{\partial t} \vspace{0.5cm} \
\dfrac{\partial {M}}{\partial t} \ = \ - K\Bigl({M} - M_e\Bigr) ,\vspace{0.3cm} \
\end{cases}
\end{equation}
Except T_{alpha}, T_{g}, W and M which represent the unknown functions of our system, all others parameters are constant.
I am beginner in Fenics but not in solving PDEs with FEM. I have read the Fenics documentation and had a look on the implementation of Cahn-Hilliard (in Fenicsx) and Advection-Diffusion-Reaction equation (in Fenics).
However I get the following error in the weak form:
Shapes do not match: and .
ERROR:UFL:Shapes do not match: and .
Traceback (most recent call last):
File “/home/mamadou/fenicsx-code/meeting28.py”, line 200, in
F3 = inner(((W_k - W_k1) / dt), v_3) * dx + inner(V_x * grad(W_k), v_3) * ufl.dx
File “/home/mamadou/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/operators.py”, line 159, in inner
return Inner(a, b)
File “/home/mamadou/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/tensoralgebra.py”, line 147, in new
error(“Shapes do not match: %s and %s.” % (ufl_err_str(a), ufl_err_str(b)))
File “/home/mamadou/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/log.py”, line 135, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: and .
(fenicsx) mamadou@it-latitude5431:~/fenicsx-code$
Please find my implementation below:
import os
import numpy as np
import matplotlib as plt
import ufl
import dolfinx
from dolfinx import log, plot
from dolfinx.fem import Function, FunctionSpace, Constant
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_interval
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, SpatialCoordinate
from mpi4py import MPI
from petsc4py import PETSc
import pyvista as pv
import pyvistaqt as pvqt
from pyvistaqt import BackgroundPlotter
#import pyvistaqt as pvqt
#from pyvistaqt import BackgroundPlotter
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
clears the terminal and prints dolfinx version
os.system(‘clear’)
prints dolfinx version
print(f"DOLFINx version: {dolfinx.version}")
Save all logging to file
log.set_output_file(“log.txt”)
Set constants
alpha = 1.0 #grain surface area per volume
C_pa = 1.0 #specific heat of the air
C_pg = 1.0 #specific heat of the grain
C_pv = 1.0 #specific heat of the water vapour
C_pw = 1.0 #specific heat of the water
rho_g = 1.0 #specific mass of the grain
rho_a = 1.0 #specific mass of the air
V_x = 1.0 #air velocity
epsilon = 0.7 #porosity
H_v = 1.0
K = 4.13
M_e = 0.15
nu = 1.0
T = 1.0 # final time
num_steps = 5000 # number of time steps
dt = T / num_steps # time step size
L = 1.0 #total lengh
#NOTE: floats must be converted to dolfin constants on domain below
#Creating an interval mesh
NUM_ELEM = 10
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])
print (“Creating an interval mesh”)
################################################################
ENTER MATERIAL PARAMETERS AND CONSTITUTIVE MODEL
#################################################################
x = SpatialCoordinate(domain) # Coordinates of domain used for defining general functions
alpha = Constant(domain, alpha) #grain surface area
C_pa = Constant(domain, C_pa) #specific heat of the air
C_pg = Constant(domain, C_pg) #specific heat of the grain
C_pv = Constant(domain, C_pv)#specific heat of the water vapour
C_pw = Constant(domain, C_pw) #specific heat of the water
rho_g = Constant(domain, rho_g) #specific mass of the grain
rho_a = Constant(domain, rho_a) #specific mass of the air
V_x = Constant(domain, V_x) ##air velocity
epsilon = Constant(domain, epsilon) #porosity
H_v = Constant(domain, H_v)
K = Constant(domain, K)
M_e = Constant(domain, M_e)
nu = Constant(domain, nu)
dt = Constant(domain, dt)
#finite element function space on domain, with trial and test functions
el1 = ufl.FiniteElement(“CG”, domain.ufl_cell(), 1) #finite element function space on domain
el2 = ufl.FiniteElement(“CG”, domain.ufl_cell(), 1) #finite element function space on domain
el3 = ufl.FiniteElement(“CG”, domain.ufl_cell(), 1) #finite element function space on domain
el4 = ufl.FiniteElement(“CG”, domain.ufl_cell(), 1) #finite element function space on domain
mix = ufl.MixedElement([el1, el2, el3, el4])
Z = FunctionSpace(domain, mix) #mixed function space
print(“Number of DOFs: %d” % Z.dofmap.index_map.size_global)
print(“Number of elements (intervals): %d” % NUM_ELEM)
print(“Number of nodes: %d” % (NUM_ELEM+1))
#Once the space has been created, we need to define
#our test functions and finite element functions. Test functions for a mixed function space can
be created by replacing TestFunction by TestFunctions:
#(T_a, T_g, W, M) = ufl.TrialFunctions(Z) #finite element functions
v_1, v_2, v_3, v_4 = ufl.TestFunctions(Z) ##Test functions
s = ufl.TrialFunction(Z)
#These functions will be used to represent the unknowns T_a_k, T_g_k, W_k, and M_k at the new time
level k. The corresponding values at the previous time level k - 1 are denoted by T_a_k1,
T_g_k1, W_k1, M_k1 in our program.
#(T_a_k, T_g_k, W_k, M_k) = ufl.TrialFunctions(Z) #the unknowns T_a_k, T_g_k, W_k, and M_k at
#the new time level k.
#(T_a_k1, T_g_k1, W_k1, M_k1) = ufl.TrialFunctions(Z) #the unknowns T_a_k1, T_g_k1, W_k1, and M_k1 at
#the new time level k - 1.
+
u = Function(Z) # current solution
u0 = Function(Z) # solution from previous converged step
Split mixed functions
s1, s2, s3, s4 = ufl.split(s)
T_a_k, T_g_k, W_k, M_k = ufl.split(u)
T_a_k1, T_g_k1, W_k1, M_k1 = ufl.split(u0)
#T_a_k = ufl.TrialFunction(el1)
#v_1 = ufl.TestFunction(el1)
#print(grad(T_a_k).ufl_shape)
#print(v_1.ufl_shape)
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 - lmbda * inner(grad(c), grad(v)) * dx
#F = F0 + F1
This is a statement of the time-discrete equations presented as part of
the problem statement, using UFL syntax.
```{index} single: Newton solver; (in Cahn-Hilliard demo)
```
The DOLFINx Newton solver requires a
{py:class}NonlinearProblem<dolfinx.fem.NonlinearProblem> object to
solve a system of nonlinear equations
#For the current example, we will set four different boundary conditions. We will set
T_a, T_g, W, M=0 at the walls of the channel; that is, at x=0 and x=1.
These boundary conditions may be defined as follows:
Define boundaries for wheat drying problem
#T_a_walls = ‘near(x[1], 0) || near(x[1], 1)’ #Air boundary
#T_g_walls = ‘near(x[1], 0) || near(x[1], 1)’ #Grain boundary
#W_walls = ‘near(x[1], 0) || near(x[1], 1)’ #Humidity boundary
#M_walls = ‘near(x[1], 0) || near(x[1], 1)’ #Moisture boundary
Define boundary conditions for wheat drying problem
#bcu_noslip1 = DirichletBC(V, Constant((0, 0)), T_a_walls)
#bcu_noslip2 = DirichletBC(V, Constant((0, 0)), T_g_walls )
#bcu_noslip3 = DirichletBC(V, Constant((0, 0)), W_walls)
#bcu_noslip4 = DirichletBC(V, Constant((0, 0)), M_walls )
#bcu = [bcu_noslip1, bcu_noslip2, bcu_noslip3, bcu_noslip4]
#When now all functions and test functions have been defined, we can express the nonlinear
variational grain drying problem :
Weak statement of the equations
print(grad(T_a_k).ufl_shape)
print(v_1.ufl_shape)
print(grad(v_1).ufl_shape)
print(T_a_k.ufl_shape)
#T= dolfinx.fem.form([V_x*inner(grad(T_a_k), v_1) * dx])
#print(grad(v_1).ufl_shape)
#print(W_k.ufl_shape)
print(grad(W_k).ufl_shape)
print(v_3.ufl_shape)
#T= V_x*inner(grad(W_k), v_3) * dx
print(grad(W_k).ufl_shape)
print(v_2.ufl_shape)
#T= V_x*inner(grad(W_k), v_3) * dx
#print(grad(T_a_k)._ufl_shape)
#print(V_x.ufl_shape)
#print(V_x._ufl_shape)
#print(alpha.ufl_shape)
#print(grad(T_a_k).ufl_shape)
#print(grad(T_a_k).ufl_shape)
#print(alpha._ufl_shape)
#print(alpha.ufl_shape)
F1 = inner(rho_a * (C_pa + C_pv * W_k) * ((T_a_k - T_a_k1) / dt), v_1 ) * dx #+ V_x * inner(grad(T_a_k), v_1) * dx
F2 = inner(rho_g * (C_pg + C_pw * M_k) * ((T_g_k - T_g_k1) / dt) ,v_2 ) * dx
F3 = inner(((W_k - W_k1) / dt), v_3) * dx + inner(V_x * grad(W_k), v_3) * ufl.dx
F4 = inner(((M_k - M_k1) / dt), v_4) * dx
G1 = inner((((alpha * nu) * (T_g_k - T_a_k)) / epsilon), v_1) * dx - inner((((alpha * nu) * (T_a_k - T_g_k)) / (1 - epsilon)), v_2) * dx
G2 = inner(((H_v + C_pv * (T_a_k - T_g_k)) * (((rho_a * epsilon * V_x * grad(W_k))) / (1 - epsilon))), v_2) * dx
G3 =inner(((((epsilon - 1) * (rho_g * ((M_k - M_k1) / dt))) / (rho_a * epsilon))), v_3) * dx - inner(K * (M_e - M_k), v_4) * dx
F = F1 + F2 + F3 + F4 - G1 + G2 - G3
F = F1 + F2 + F3 + F4 - G1 - G3
The DOLFINx Newton solver requires a
{py:class}NonlinearProblem<dolfinx.fem.NonlinearProblem> object to
solve a system of nonlinear equations
+
Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
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(domain)