# Porting over a Coupled Linear Elasticity-Reaction Diffusion Problem from Legacy FEnicS to FEnicSx

Hey all,

I’ve just been trying to translate some coupled linear elasticity reaction diffusion code from legacy FEnicS to FEnicSx. Originally wasn’t planning to post this here due to the length of the two code pieces but am at wits’ end on what to do… I was wondering if someone could have a quick look over it to see if the problem stems from how I’m generating functions and spaces.

The legacy code produces the desired result of the mesh moving in the positive x-direction (using Paraviews Warp Vector Field Function). However, the FEnicSx code does not, even though the problem and weak formulations are identical.

The Legacy FEnicS code is as follow and produces the desired result:

``````from fenics import *
from mshr import *
from ufl import nabla_div, transpose, inv
import numpy as np
import os
import matplotlib.tri as tri

def main():

#_________________________________________________________________________
#
# C O N T R O L   P A N E L :   I N P U T S
#_________________________________________________________________________
#

# Simulation Displacement Type
sim_multiple = True    # False - Every simulation outputs to output/main (overwriting)
# True  - Every simulation generates new output directory to output/main_i (i = 1, 2, 3, ...)

# Simulation Runtime Parameters
T = 30                 # final time
num_steps = T*10        # number of time steps
dt = T / num_steps      # time step size
output_freq = 1        # number of time steps per output

# Geometric Parameters
mesh_resolution = 32                # Higher number is higher resolution
vol = np.pi*mesh_radius**2          # Initial Volume

# Wave Pinning Kinetic Parameters
b_0 = 4.0                # Basal Activation
b_1 = 0.1                # Chemical Gradient on Basal Activation
delta_0 = 2.0            # Constant Negative Feedback
delta_1 = 50.0           # Area-dependent Negative Feedback
b2 = 5.0                 # Magnitude of Positive Feedback
n = 6                    # Hill's Coeefficient

# Diffusivity
Da = 0.01               # Active Rac Diffusivity
Di = 10.0               # Inactive Rac Diffusivity

# Elasticity Parameters (sim_disp == 2)
T_scale = 0.3                      # Scaling Factor for Polymerization Force
s_0 = 1.0
s_1 = 10.0
E = 0.3                             # Youngs Modulus
nu = 0.45                           # Poisson's Ratio
beta = 10.0                         # Internal Damping: Viscoelastic Modulus
gamma = 20.0                        # External Damping: Stokes Damping
lambda_ = E*nu/((1+nu)*(1-2*nu))    # Lame's First Parameter
mu = 0.5*E/(1+nu)                   # Lame's Second Parameter (Shear Modulus)

#_________________________________________________________________________
#
# F E M   D O M A I N   S E T U P
#_________________________________________________________________________

#-------------------------------------------------------------------------
# M E S H
#-------------------------------------------------------------------------

# Create mesh and define function space
mesh = generate_mesh(domain, mesh_resolution)

#-------------------------------------------------------------------------
# F U N C T I O N   S P A C E S
#-------------------------------------------------------------------------

# Define function space and elements for displacments and concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)        # Function space for concentrations
W = VectorFunctionSpace(mesh, 'P', 2)   # Function space for displacments

#-------------------------------------------------------------------------
# V A R I A T I O N S
#-------------------------------------------------------------------------

# Define variational problem
u = Function(V)
v_1, v_2 = TestFunctions(V)

w = TestFunction(W)

#_________________________________________________________________________
#
# R E A C T I O N   D I F F U S I O N
#_________________________________________________________________________
#

#-------------------------------------------------------------------------
# I N I T I A L   C O N D I T I O N S
#-------------------------------------------------------------------------

# Define initial conditions
u_0 = Expression(('0.0', '2.0'), degree=1)
u_n = project(u_0, V)

#-------------------------------------------------------------------------
# S O U R C E S
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# S P L I T T I N G
#-------------------------------------------------------------------------

# Split system functions to access components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)

#-------------------------------------------------------------------------
# C O N S T I T U T I V E   E Q N :   R E A C T I O N   K I N E T I C S
#-------------------------------------------------------------------------

# Negative Feedback
delta = Expression('delta_0 + (vol-vol_init)*delta_1', \
delta_0=delta_0, vol=vol, vol_init=vol, delta_1=delta_1, degree=2)

# Fields for Kinetic Parameters
b2 = Constant(b2)
n = Constant(n)
Da = Constant(Da)
Di = Constant(Di)

#_________________________________________________________________________
#
# M E C H A N I C S
#_________________________________________________________________________

#-------------------------------------------------------------------------
# I N I T I A L   C O N D I T I O N S
#-------------------------------------------------------------------------

# Define initial value
disp_D = Expression(('0.0','0.0'),degree=2)

#-------------------------------------------------------------------------
# S O U R C E    T E R M S
#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
# K I N E M A T I C S
#-------------------------------------------------------------------------

# Need the following defined for all three simulations modes (sim_disp)
disp = TrialFunction(W)
disp_t1 = interpolate(disp_D, W)

d = disp_t1.geometric_dimension()           # Define space dimension
I = Identity(d)                             # Identity Tensor (Matrix)

vol_init = assemble(J*dx)
delta.vol_init = vol_init
delta.vol = vol_init

# Kinetmatics at previous time step
disp_t2 = interpolate(disp_D, W)
I_n = Identity(disp_t2.geometric_dimension())

b = Expression('b_0 + x[0]*b_1', b_0=b_0,b_1=b_1,degree=2)

# Wave-Pinning Model
G = (b + b2*(pow(u_1,n)/(1+pow(u_1,n))))*u_2 - delta*u_1

#-------------------------------------------------------------------------
# B O U N D A R Y   C O N D I T I O N S
#-------------------------------------------------------------------------

# Actin Polymerization Magnitude of Force (traction)
T_mag = T_scale /(1+exp(-2*s_1*(u_n[0]-s_0)))

# Unit Normal Vector on Reference Domain
Normal_Ref = FacetNormal(mesh)
# Unit Normal Vector on Deformed Domain
Normal_Def1 = Normal_Def/sqrt(Normal_Def[0]**2+Normal_Def[1]**2)

# Actin Polymerization Force
Trac_disp = T_mag*Normal_Def1

# Small Strain Approximation
def epsilon(disp_func):

#-------------------------------------------------------------------------
# C O N S T I T U T I V E   E Q N :   S T R E S S - S T R A I N
#-------------------------------------------------------------------------

# Linear Elasticity
def sigma_p(disp_func):
return lambda_*nabla_div(disp_func)*Identity(d) + 2*mu*epsilon(disp_func)

# Define expressions used in variational forms
beta = Constant(beta)
gamma = Constant(gamma)

#-------------------------------------------------------------------------
# R D   B O U N D A R Y   C O N D I T I O N S
#-------------------------------------------------------------------------
# _______________________________________________________________________
#
# F U N C T I O N A L S
#_________________________________________________________________________
#

# Time Step
k = Constant(dt)

#Displacement
F_disp_damping = (gamma / k) * inner((disp - disp_t2),w)*dx
F_disp_elasticity = inner(sigma_p(disp), epsilon(w))*dx
F_disp_viscous = (beta / k) * inner(sigma_p(disp - disp_t2), epsilon(w))*dx
F_disp_traction = - dot(Trac_disp, w)*ds

F_disp = F_disp_damping + F_disp_elasticity + F_disp_viscous + F_disp_traction
F_disp_lhs = lhs(F_disp)
F_disp_rhs = rhs(F_disp)

disp = Function(W)

#Active Rac
F_1_temporal  =   J*((u_1 - u_n1) / k)*v_1*dx
F_1_reaction  = - J*G*v_1*dx
F_1_dilution  =   u_1*v_1*((J - J_n) / (k))*dx
F_1 = F_1_temporal + F_1_diffusion + F_1_reaction + F_1_dilution

#Inactive Rac
F_2_temporal  =   J*((u_2 - u_n2) / k)*v_2*dx
F_2_reaction  = + J*G*v_2*dx
F_2_dilution  =   u_2*v_2*((J - J_n) / (k))*dx
F_2 = F_2_temporal + F_2_diffusion + F_2_reaction + F_2_dilution

# Functional for Conservation of Mass
F_rd = F_1 + F_2

#_________________________________________________________________________
#
# S O L V I N G
#_________________________________________________________________________
#

#-------------------------------------------------------------------------
# O U T P U T    S E T U P
#-------------------------------------------------------------------------

dir_count = 1
output_dir_original = 'output/main'
output_dir = output_dir_original

if sim_multiple == True:
while os.path.isdir(output_dir):
output_dir = output_dir_original + "_" + str(dir_count)
dir_count += 1

# Outputs
ActiveGTPase_file = XDMFFile(output_dir + '/ActiveGTPase.xdmf')
InActiveGTPase_file = XDMFFile(output_dir + '/InActiveGTPase.xdmf')
Displacement_file = XDMFFile(output_dir + '/Displacement.xdmf')

#-------------------------------------------------------------------------
# I N I T I A L   C O N D I T I O N S   O U T P U T
#-------------------------------------------------------------------------

# Time-stepping
t = 0
i = 0

#Output Labels
_u_1, _u_2 = u_n.split()
_u_1.rename("ActiveGTPase","label")
_u_2.rename("InActiveGTPase","label")
disp.rename("Displacement","label")

# Calculate Volume and Update Feedback
disp_t1.assign(disp)
vol = assemble(J*dx)    #dv = JdV, Jacobian relates Reference Volume to Deformed Volume
delta.vol = vol         #Updating Delta Feedback
total_rac = assemble((u_1+u_2)*J*dx)

#Writing Outputs
ActiveGTPase_file.write(_u_1,t)
InActiveGTPase_file.write(_u_2,t)
Displacement_file.write(disp,t)

if sim_multiple == True:
os.system('mkdir ' + output_dir + '/frames')

#-------------------------------------------------------------------------
# T I M E    L O O P
#-------------------------------------------------------------------------

for n in range(1,num_steps+1):

# Update current time
t += dt

# Mechanical Deformation Solver
solve(F_disp_lhs == F_disp_rhs, disp)

# Calculate Volume and Update Feedback
disp_t1.assign(disp)
vol = assemble(J*dx)    #dv = JdV, Jacobian relates Reference Volume to Deformed Volume
delta.vol = vol         #Updating Delta Feedback

# Reaction Diffusion Solver
solve(F_rd == 0, u)

# Useful Info
total_rac = assemble((u_1+u_2)*J*dx)
print("Time Step {} of {}".format(n,num_steps))
print("Current Volume: ", vol)
print("total Rac: ", total_rac)

# Save to file and plot solution
if n%output_freq == 0:
_u_1, _u_2 = u.split()
# Update current output indent
i += 1
_u_1.rename("ActiveGTPase","label")
_u_2.rename("InActiveGTPase","label")
disp.rename("Displacement","label")
# Writing Outputs
ActiveGTPase_file.write(_u_1,t)
InActiveGTPase_file.write(_u_2,t)
Displacement_file.write(disp,t)

# Update previous solution
u_n.assign(u)
disp_t2.assign(disp)

if __name__ == '__main__':
main()

``````

The FEnicSx translation is as follows:

``````#----------------------
# IMPORTS
#----------------------
import ufl
import numpy as np
import gmsh

from petsc4py import PETSc
from mpi4py import MPI

import dolfinx
from dolfinx import fem, mesh, plot, log, default_scalar_type, geometry
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_matrix, create_vector, create_matrix, assemble_vector, apply_lifting
from dolfinx.nls.petsc import NewtonSolver
import basix

#----------------------
# CREATE CIRCLE MESH
#----------------------

# Initialise the GMSH module
gmsh.initialize()

# Create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures
# Create Circular Mesh with area of A0 = 1
gmsh.model.occ.synchronize()

# Make membrane a physical surface for GMSH to recognise when generating the mesh
gdim = 2 # 2D Geometric Dimension of the Mesh
gmsh.model.addPhysicalGroup(gdim, [membrane], 0) # Dimension, Entity tag, Physical tag

# Generate 2D Mesh with uniform mesh size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.023)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.023)
gmsh.model.mesh.generate(gdim)

# Create a mesh from the GMSH model
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
domain.name = "initial_mesh"
gmsh.write("circle.msh")

# Finalize GMSH
gmsh.finalize()

#----------------------
# DEFINE COORDINATES
#----------------------

# Global Spatial Coordinates as a Vector-Valued Expression i.e. x[0] = x, x[1] = y
x_spatial = ufl.SpatialCoordinate(domain)

# Define Measure Over Cells of Mesh

# Define Measure Over Boundary of Mesh

# Time Discretisation Process
t = 0 # Initial Time
T = 30 # Final Time
dt_original = 0.1
dt = fem.Constant(domain, dt_original) # Time Step Size
Nsteps = int(T/dt_original) # Number of Time Steps

#---------------------------------
# LINEAR ELASTICITY FUNCTION SPACE
#---------------------------------

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, ))) # mesh, element

# Trial and Test Functions for Linear Elasticity Problem (Mechanics Vectors)
# Initial Conditions: Mechanical System is at Rest
u = ufl.TrialFunction(V)# Trial Function (Displacement)
v = ufl.TestFunction(V) # Test Function

u_ = fem.Function(V) # Function for Area Calculation
u_.name = "u_"

u_n = fem.Function(V) # Function for Previous Time Step (Displacement)
u_n.name = "u_n"

u_0 = fem.Constant(domain, 0.0)

u_n.sub(0).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_n.sub(1).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_n.sub(0).collapse()
u_n.sub(1).collapse()
u_n.x.scatter_forward()

u_.sub(0).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_.sub(1).interpolate(lambda x: np.full((x.shape[1],), u_0))
u_.sub(0).collapse()
u_.sub(1).collapse()
u_.x.scatter_forward()

#---------------------------------
# REACTION DIFFUSION FUNCTION SPACE
#---------------------------------

# Mixed-Element Function Space for Reaction Diffusion
P1 = basix.ufl.element(family="Lagrange", cell=domain.topology.cell_name(), degree=1)
element = basix.ufl.mixed_element([P1, P1])
C = fem.functionspace(mesh=domain, element=element)

# Trial and Test Functions for Reaction-Diffusion Problem (Concentrations)
c = fem.Function(C)
c.name = "c"

c_n = fem.Function(C)
c_n.name = "c_n"

c_A, c_I = ufl.split(c)
c_A_n, c_I_n = ufl.split(c_n)

w_A, w_I = ufl.TestFunction(C)

# Initial Species Concentration Conditions of Species I=2 and Species A=0
c_A_0 = fem.Constant(domain, 0.0)
c_I_0 = fem.Constant(domain, 2.0)

c.sub(0).interpolate(lambda x: np.full((x.shape[1],), c_A_0))
c.sub(1).interpolate(lambda x: np.full((x.shape[1],), c_I_0))
c_n.sub(0).interpolate(lambda x: np.full((x.shape[1],), c_A_0))
c_n.sub(1).interpolate(lambda x: np.full((x.shape[1],), c_I_0))

c.sub(0).collapse()
c.sub(1).collapse()
c_n.sub(0).collapse()
c_n.sub(1).collapse()

#------------------------------------
# LINEAR ELASTICITY MECHANICS PROBLEM
#------------------------------------

# Spatial Dimension
spa_dim = len(x_spatial)

# Current Time Step
# Identity Tensor
I = ufl.Identity(spa_dim)

F_inv = ufl.inv(F) # F^-1
F_inv_tra = F_inv.T # F^-T

# Jacobian Function
J = ufl.det(F)

# Previous Time Step
# Identity Tensor
I_n = ufl.Identity(spa_dim)

F_n_inv = ufl.inv(F_n) # F^-1
F_n_inv_tra = F_n_inv.T # F^-T

# Jacobian Function
J_n = ufl.det(F_n)

# Initial Cell Area (i.e. A0=1) and Tension Calculations (i.e. T0=0)
A0 = fem.assemble_scalar(fem.form(fem.Constant(domain, 1.0) * dx))
Area = fem.assemble_scalar(fem.form(J * dx))
Tension = Area - A0

# Calculate Normal Vectors along Boundary in Reference/Undeformed Configuration (might need minus to point outwards)
# n0 = ufl.FacetNormal(domain) # Stationary Domain Normal Vector
theta = ufl.operators.atan2(x_spatial[1], x_spatial[0]) # Unit Normal Vector on Reference Domain
n0 = ufl.as_vector([ufl.cos(theta), ufl.sin(theta)])
piola_transform = J * F_inv_tra * n0
n = piola_transform/ufl.sqrt(piola_transform[0]**2 + piola_transform[1]**2) # Unit Normal Direction of Membrane

# Mechanics Constants
youngs_modulus_E = 0.30 # kPa
poisson_ratio_nu = 0.45
stoke_drag_gamma = fem.Constant(domain, 20.0) # kPa*s/um^2 for final velocity of 0.17 um/s
lambda_ = fem.Constant(domain, youngs_modulus_E*poisson_ratio_nu/((1+poisson_ratio_nu)*(1-2*poisson_ratio_nu))) # ~0.93 kPa
mu_ = fem.Constant(domain, youngs_modulus_E/(2*(1+poisson_ratio_nu))) # ~0.1 kPa
beta_ = fem.Constant(domain, 10.0) # s

# Polymerisation Force Constants
s0 = fem.Constant(domain, 1.0)
s1 = fem.Constant(domain, 10.0)
f_R = fem.Constant(domain, 0.3) # nN/um

# Cauchy Stress Tensor
def sigma(u_current):
eps_current = epsilon(u_current)
stress = lambda_ * ufl.nabla_div(u_current) * ufl.Identity(len(u_current)) + 2 * mu_ * eps_current
# stress = lambda_ * ufl.tr(eps_current) * ufl.Identity(len(u_current)) + 2 * mu_ * eps_current
return stress

# Small-Scale Strain
def epsilon(u):
return small_strain

# Heaviside Polymerisation Force from Actin
polymerisation_force = f_R / (1 + ufl.exp(-2 * s1 * (c_A - s0)))

#------------------------------------
# LINEAR ELASTICITY WEAK FORM
#------------------------------------

# Linear Elasticity Weak Form
damping = (stoke_drag_gamma/dt) * ufl.inner((u - u_n), v) * dx # Damping Term #SHOULD IT BE NEGATIVE? OR IS GAMMA NEGATIVE?
elastic = ufl.inner(sigma(u), epsilon(v)) * dx
viscous = (beta_/dt) * ufl.inner(sigma(u - u_n), epsilon(v)) * dx
traction = - ufl.inner(polymerisation_force * n, v) * ds # Traction Term
#body = 0.0 # No source term

# Create Linear Problem and Solver
mechanics_formulation = viscous + elastic + traction + damping
mechanics_lhs = fem.form(ufl.lhs(mechanics_formulation))
mechanics_rhs = fem.form(ufl.rhs(mechanics_formulation))

lhs_bilinear_matrix = create_matrix(mechanics_lhs)
rhs_linear_vector = create_vector(mechanics_rhs)

mechanics_solver = PETSc.KSP().create(domain.comm)
mechanics_solver.setOperators(lhs_bilinear_matrix)

#----------------------------------
# REACTION DIFFUSION FUNCTION SPACE
#----------------------------------

class DeltaExpression:
def __init__(self, delta_0, area, area_init, delta_1):
self.delta_0 = delta_0
self.delta_1 = delta_1
self.area = area
self.area_init = area_init

def eval(self, x):
return np.full(x.shape[1], self.delta_0 + (self.area - self.area_init) * self.delta_1)

delta_0 = 2.0 #fem.Constant(domain, 2.0)
delta_1 = 30.0 #fem.Constant(domain, 30.0)
delta_expr = DeltaExpression(delta_0, Area, A0, delta_1)

# Create a Function to hold the values of delta
M = fem.functionspace(domain, ("CG", 2))
delta = fem.Function(M)
delta.interpolate(delta_expr.eval)

def __init__(self, b_0, b_1):
self.b_0 = b_0
self.b_1 = b_1

def eval(self, x):
return np.full(x.shape[1], self.b_0 + x[0] * self.b_1)

b_0 = 4.0 # Basal Activation Rate
b_1 = 0.1

b = fem.Function(M)

c_hill = fem.Constant(domain, 5.0) # Hill Function Magnitude
cA_0 = fem.Constant(domain, 1.0)
n = fem.Constant(domain, 6.0) # Hill's Coefficient

# Reaction Function (Identical for Both Inactivate and Active GTPase)
reaction_function = (b + c_hill * (c_A**n/(cA_0**n + c_A**n))) * c_I - delta * c_A

# Examine Isotropic Diffusion where D_A=D*I (Diffusivity Constant multiplied by the Identity Tensor)
D_a = fem.Constant(domain, 0.01)
D_i = fem.Constant(domain, 10.0)

D_A = D_a * I
D_A_eq = J * F_inv * D_A * F_inv_tra
D_I = D_i * I
D_I_eq = J * F_inv * D_I * F_inv_tra

u = fem.Function(V)

#------------------------------------
# REACTION DIFFUSION WEAK FORM
#------------------------------------

# Derive Weak/Variational Formulation of 2D Reaction-Diffusion PDE with Two Chemical Species
term1_A = J * ((c_A - c_A_n)/dt) * w_A * dx # Temporal Term
term2_A = ufl.dot(D_A_eq * ufl.grad(c_A), ufl.grad(w_A)) * dx # Diffusion+Flux+Source? (Mass Conservation)
term3_A = -J * reaction_function * w_A * dx # Reaction Term
term4_A = c_A * ((J - J_n)/dt) * w_A * dx # Dilution Term

Func_A = term1_A + term2_A + term3_A + term4_A

term1_I = J * ((c_I - c_I_n)/dt) * w_I * dx # Temporal Term
term3_I = J * reaction_function * w_I * dx # Reaction Term
term4_I = c_I * ((J - J_n)/dt) * w_I * dx # Dilution Term

Func_I = term1_I + term2_I + term3_I + term4_I

reaction_diffusion_formulation = Func_A + Func_I

# Non-Linear Reaction-Diffusion Problem and Solver Formulation
reaction_problem = NonlinearProblem(reaction_diffusion_formulation, c, bcs=[])

reaction_solver = NewtonSolver(mesh_comm, reaction_problem)
reaction_solver.atol = 1e-12
reaction_solver.rtol = 1e-12
reaction_solver.convergence_criterion = "incremental"
reaction_solver.report = True

#----------------------
# SOLVE COUPLED PROBLEM
#----------------------

# Create XDMF file for writing
xdmf = dolfinx.io.XDMFFile(mesh_comm, "jarod_move2.xdmf", "w")
xdmf.write_mesh(domain)

# Write initial function to XDMF file
xdmf.write_function(u_n, 0.0)
xdmf.write_function(c.sub(0), 0.0)
xdmf.write_function(c.sub(1), 0.0)

for i in range(Nsteps):

# Update time step
t += dt

# Assemble Jacobian and residual
with rhs_linear_vector.localForm() as loc_1:
loc_1.set(0)
assemble_vector(rhs_linear_vector, mechanics_rhs)

lhs_bilinear_matrix.zeroEntries()
assemble_matrix(lhs_bilinear_matrix, mechanics_lhs)
lhs_bilinear_matrix.assemble()
mechanics_solver.setOperators(lhs_bilinear_matrix)

# # Step 1: Solve Linear Elasticity Problem
mechanics_solver.solve(rhs_linear_vector, u.vector)

u.x.scatter_forward()
u_.x.array[:] = u.x.array #Might NEED [:] on the right side too

# Step 2: Compute new area and tension
Area = fem.assemble_scalar(fem.form(J * dx))
Tension = Area - A0
delta.area = Area

# Step 3: Solve Reaction-Diffusion PDE Problem
num_its, converged_reaction = reaction_solver.solve(c)
assert converged_reaction

c.x.scatter_forward()

c_n.x.array[:] = c.x.array # Update Concentrations at previous time step
u_n.x.array[:] = u.x.array # Update Displacement at previous time step

total_rac = fem.assemble_scalar(fem.form((c_A + c_I) * J * dx))
active_rac = fem.assemble_scalar(fem.form((c_A) * J * dx))
inactive_rac = fem.assemble_scalar(fem.form((c_I) * J * dx))

print(f'Original Area: {A0}')
print(f'Area: {Area}')
print(f'Tension: {Tension}')
print(f'Total Rac: {total_rac}')
print(f'Active Rac: {active_rac}')
print(f'Inactive Rac: {inactive_rac}')

xdmf.write_function(u_n, float(t))
xdmf.write_function(c.sub(0), float(t))
xdmf.write_function(c.sub(1), float(t))

# Close Save Files
xdmf.close()
``````

Legacy FEnicS was installed from the docker using the latest build. FEnicSx was installed via micromamba and is v0.8.0.

I assume a potential area where the problem could stem from is how the function spaces are set up and whether the values are being updated properly in the time loop as I’ve noticed total_rac is not conserved and slowly increases in the FEnicSx code. In this case, are my time-dependent variables (Area, Jacobian J and J_n, delta) being updated properly and therefore updating the problem (linear problem matrices and vectors, and reaction-diffusion J, J_n, and reaction_function) to be solved at each time step properly?

Additionally, in the FEnicSx code are you supposed to use classes instead of Expression now?

Cheers,
Volkan

Ended up fixing this in a less than ideal way. Last ditch effort was to just recalculate weak forms and values in the time loop. Unfortunately, this is a lot slower and not ideal for readability.

If anyone knows how FEnicSx can automatically update values for a Jacobian and functions which are updated every time loop that would be greatly appreciated.

Here is the new time loop:

``````#Output Labels
_c_A, _c_I = c_n.split()

# Create XDMF file for writing
xdmf = dolfinx.io.XDMFFile(mesh_comm, "cell_move2.xdmf", "w")
xdmf.write_mesh(domain)

# Write initial function to XDMF file
xdmf.write_function(u_n, 0.0)
xdmf.write_function(_c_A, 0.0)
xdmf.write_function(_c_I, 0.0)

for i in tqdm(range(Nsteps), desc="Solving PDE"):

# Update time step
t += dt

# Heaviside Polymerisation Force from Actin
polymerisation_force = f_R / (1 + ufl.exp(-2 * s1 * (c.sub(0) - s0)))

# Linear Elasticity Weak Form
damping = (stoke_drag_gamma/dt) * ufl.inner((u - u_n), v) * dx # Damping Term #SHOULD IT BE NEGATIVE? OR IS GAMMA NEGATIVE?
elastic = ufl.inner(sigma(u), epsilon(v)) * dx
viscous = (beta_/dt) * ufl.inner(sigma(u - u_n), epsilon(v)) * dx
traction = - ufl.dot(polymerisation_force*n, v) * ds # Traction Term
#body = 0.0 # No source term

# Create Linear Problem and Solver
mechanics_formulation = viscous + elastic + traction + damping
mechanics_lhs = fem.form(ufl.lhs(mechanics_formulation))
mechanics_rhs = fem.form(ufl.rhs(mechanics_formulation))

lhs_bilinear_matrix = assemble_matrix(mechanics_lhs, bcs=[])
lhs_bilinear_matrix.assemble()
rhs_linear_vector = assemble_vector(mechanics_rhs)

mechanics_solver = PETSc.KSP().create(domain.comm)
mechanics_solver.setOperators(lhs_bilinear_matrix)

# # Step 1: Solve Linear Elasticity Problem
mechanics_solver.solve(rhs_linear_vector, u_.vector)

u_.x.scatter_forward()

# Current Time Step
# Identity Tensor
I = ufl.variable(ufl.Identity(spa_dim))

F_inv = ufl.variable(ufl.inv(F)) # F^-1
F_inv_tra = ufl.variable(ufl.transpose(F_inv)) # F^-T

# Jacobian Function
J = ufl.variable(ufl.det(F))

# Step 2: Compute new area and tension
Area = fem.assemble_scalar(fem.form(J * dx))
Tension = Area - A0

delta = delta_0 + Tension * delta_1 # THIS WAS THE FUCKING ISSUE????

# Reaction Function (Identical for Both Inactivate and Active GTPase)
reaction_function = (b + c_hill * (c_A**n_power/(cA_0**n_power + c_A**n_power))) * c_I - delta * c_A
D_A_eq = J * F_inv * D_A * F_inv_tra
D_I_eq = J * F_inv * D_I * F_inv_tra

# Derive Weak/Variational Formulation of 2D Reaction-Diffusion PDE with Two Chemical Species
term1_A = J * ((c_A - c_A_n)/dt) * w_A * dx # Temporal Term
term2_A = ufl.inner(D_A_eq * ufl.grad(c_A), ufl.grad(w_A)) * dx # Diffusion+Flux+Source? (Mass Conservation)
term3_A = -J * reaction_function * w_A * dx # Reaction Term
term4_A = c_A * ((J - J_n)/dt) * w_A * dx # Dilution Term

Func_A = term1_A + term2_A + term3_A + term4_A

term1_I = J * ((c_I - c_I_n)/dt) * w_I * dx # Temporal Term
term3_I = J * reaction_function * w_I * dx # Reaction Term
term4_I = c_I * ((J - J_n)/dt) * w_I * dx # Dilution Term

Func_I = term1_I + term2_I + term3_I + term4_I

reaction_diffusion_formulation = Func_A + Func_I
# Non-Linear Reaction-Diffusion Problem and Solver Formulation
reaction_problem = NonlinearProblem(reaction_diffusion_formulation, c, bcs=[])

reaction_solver = NewtonSolver(domain.comm, reaction_problem)
reaction_solver.atol = 1e-8
reaction_solver.rtol = 1e-8
reaction_solver.convergence_criterion = "incremental"
reaction_solver.report = True
# Step 3: Solve Reaction-Diffusion PDE Problem
num_its, converged_reaction = reaction_solver.solve(c)
assert converged_reaction

c.x.scatter_forward()

c_n.x.array[:] = c.x.array # Update Concentrations at previous time step
u_n.x.array[:] = u_.x.array # Update Displacement at previous time step

# Previous Time Step
# Identity Tensor
I_n = ufl.variable(ufl.Identity(spa_dim))

F_n_inv = ufl.variable(ufl.inv(F_n)) # F^-1
F_n_inv_tra = ufl.variable(F_n_inv.T) # F^-T

# Jacobian Function
J_n = ufl.variable(ufl.det(F_n))

# total_rac = fem.assemble_scalar(fem.form((c_A + c_I) * J * dx))
# active_rac = fem.assemble_scalar(fem.form((c_A) * J * dx))
# inactive_rac = fem.assemble_scalar(fem.form((c_I) * J * dx))

# print(f'Original Area: {A0}')
# print(f'Area: {Area}')
# print(f'Tension: {Tension}')
# print(f'Total Rac: {total_rac}')
# print(f'Active Rac: {active_rac}')
# print(f'Inactive Rac: {inactive_rac}')

xdmf.write_function(u_n, float(t))
xdmf.write_function(c.sub(0), float(t))
xdmf.write_function(c.sub(1), float(t))

# Close Save Files
xdmf.close()
``````

You have given us quite a lot of code to parse through, which makes this quite time-consuming.

I guess you are updating `delta` in the wrong way in the initial post.
I would suggest reading through how to update constants over time:
https://jsdokken.com/FEniCS23-tutorial/src/form_compilation.html
and I think your handling of `delta` should be similar to what I have done for `u_exact` in:
https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html
i.e. you should update `delta_expr` in the loop, `delta_expr.area = area`
and then do `delta.interpolate(delta_expr)`

1 Like

Hey dokken,

Thanks so much for replying to this!

Yes, definitely don’t/didn’t expect anyone to look at this obtusely long code, sorry about that.

It seems this was the fix I needed for my problem.

I will leave the lines I updated here for anyone looking for a mixed problem in the future.

Updated the solver:

``````# Linear Elasticity Weak Form
damping = (stoke_drag_gamma/dt) * ufl.inner((u - u_n), v) * dx # Damping Term #SHOULD IT BE NEGATIVE? OR IS GAMMA NEGATIVE?
elastic = ufl.inner(sigma(u), epsilon(v)) * dx
viscous = (beta_/dt) * ufl.inner(sigma(u - u_n), epsilon(v)) * dx
traction = - ufl.dot(polymerisation_force * n, v) * ds # Traction Term
#body = 0.0 # No source term

# Create Linear Problem and Solver
mechanics_formulation = viscous + elastic + traction + damping
mechanics_lhs = fem.form(ufl.lhs(mechanics_formulation))
mechanics_rhs = fem.form(ufl.rhs(mechanics_formulation))

petsc_options = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
mechanics_problem = dolfinx.fem.petsc.LinearProblem(mechanics_lhs, mechanics_rhs, u=u_, bcs=[], petsc_options=petsc_options)
``````

Delta expression as per your instructions:

``````class DeltaExpression:
def __init__(self, delta_0, area, area_init, delta_1):
self.delta_0 = delta_0
self.delta_1 = delta_1
self.area = area
self.area_init = area_init

def eval(self, x):
return np.full(x.shape[1], self.delta_0 + (self.area - self.area_init) * self.delta_1)

delta_0 = fem.Constant(domain, 2.0)
delta_1 = fem.Constant(domain, 50.0)
delta_expr = DeltaExpression(delta_0, Area, A0, delta_1)

M = fem.functionspace(domain, ("CG", 2))
delta = fem.Function(M)
delta.interpolate(delta_expr.eval)
``````

Deleted the class for b and replaced it with a constant term, as well as renaming n to n_power (was using n for the normal vector):

``````c_hill = fem.Constant(domain, 5.0) # Hill Function Magnitude
cA_0 = fem.Constant(domain, 1.0)
n_power = fem.Constant(domain, 6.0) # Hill's Coefficient

b_0 = fem.Constant(domain, 4.0) # Basal Activation Rate
b_1 = fem.Constant(domain, 0.1)
b = b_0 + b_1 * x_spatial[0] # Applies an Internal Chemical Gradient Perturbation Function in the X-direction

# Reaction Function (Identical for Both Inactivate and Active GTPase)
reaction_function = (b + c_hill * (c_A**n_power/(cA_0**n_power + c_A**n_power))) * c_I - delta * c_A
``````

And finally an overhaul of the time loop update:

``````#Output Labels
_c_A, _c_I = c_n.split()

# Create XDMF file for writing
xdmf = dolfinx.io.XDMFFile(mesh_comm, "cell_move.xdmf", "w")
xdmf.write_mesh(domain)

# Write initial function to XDMF file
xdmf.write_function(u_n, 0.0)
xdmf.write_function(_c_A, 0.0)
xdmf.write_function(_c_I, 0.0)

for i in tqdm(range(Nsteps), desc="Solving PDE"):

# Update time step
t += dt

# Step 1: Solve Linear Elasticity Problem
mechanics_problem.solve()
u_.x.scatter_forward()

# Step 2: Compute new area and tension
Area = fem.assemble_scalar(fem.form(J * dx))
delta_expr.area = Area
delta.interpolate(delta_expr.eval)

# Step 3: Solve Reaction-Diffusion PDE Problem
num_its, converged_reaction = reaction_solver.solve(c)
assert converged_reaction

c.x.scatter_forward()

# Step 4: Update Concentrations and Displacements
c_n.x.array[:] = c.x.array # Update Concentrations at previous time step
u_n.x.array[:] = u_.x.array # Update Displacement at previous time step

# Step 5: Write to File
xdmf.write_function(u_n, float(t))
xdmf.write_function(_c_A, float(t))
xdmf.write_function(_c_I, float(t))

# Close Save Files
xdmf.close()
``````

Thanks as always