# Simulation of Coupled Heat and Mass Transfer in Granular Porous Media

Hi dear all,

I am trying to solve the time-dependent non-linear partial differential equations using complex Fenicsx.

\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}

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):
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) Please use proper formatting ($\$ surrounding maths, 3x surrounding code), and make sure the code is a minimal reproducible example. This means removing all code not required to reproduce the results

Thanks for your time and attention dear Dokken.

\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}
\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}
\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}
\dfrac{\partial {M}}{\partial t} \ = \ - K\Bigl({M} - M_e\Bigr)

Please find below the code using Fenicsx with conda:

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)


** Test and Trial Functions**

v_1, v_2, v_3, v_4 = ufl.TestFunctions(Z) ##Test functions
u = Function(Z)  # current solution
u0 = Function(Z)  # solution from previous converged step


Split mixed functions

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)


Weak Formulation

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


When I ran the code, I have the follwing error:

Shapes do not match: <ComponentTensor id=140577986816304> and <Indexed id=140577986813584>.
ERROR:UFL:Shapes do not match: <ComponentTensor id=140577986816304> and <Indexed id=140577986813584>.
Traceback (most recent call last):
File "/home/mamadou/fenicsx-code/meeting28.py", line 200, in <module>
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: <ComponentTensor id=140577986816304> and <Indexed id=140577986813584>.


Best regards,

Hello dear Jokken,
Is my question well formulated?

Kind regards,

I don’t have the bandwidth to parse your code and lengthy formulation; however, your error states that somewhere in

    F3 = inner(((W_k - W_k1) / dt), v_3) * dx  +  inner(V_x * grad(W_k), v_3) * ufl.dx


you have terms whose shapes do not make sense. It’s likely the component inner(V_x * grad(W_k), v_3) since v_3 I believe is a scalar and V_x * grad(W_k) will be a vector or tensor.

Yes, you are right. The errors come from those terms:

V_x * inner(grad(T_a_k), v_1) * dx, inner(V_x * grad(W_k), v_3) * ufl.dx and

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
`