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_radius = 1/np.sqrt(np.pi) # Radius of Circle
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
domain = Circle(Point(0.0,0.0), mesh_radius)
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)
def_grad = I + nabla_grad(disp_t1) # Deformation Gradient Tensor
J = det(def_grad) # Jacobian
def_grad_inv = inv(def_grad) # Inverse Deformation Gradient Tensor
def_grad_inv_T = transpose(inv(def_grad)) # Inverse Transpose Deformation Gradient Tensor
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())
def_grad_n = I_n + grad(disp_t2)
J_n = det(def_grad_n)
# Chemical Gradient
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_Def = J*def_grad_inv_T*Normal_Ref
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):
return 0.5*(nabla_grad(disp_func) + nabla_grad(disp_func).T)
#-------------------------------------------------------------------------
# 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_diffusion = J*Da*dot(def_grad_inv*def_grad_inv_T*grad(u_1), grad(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_diffusion = J*Di*dot(def_grad_inv*def_grad_inv_T*grad(u_2), grad(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
membrane = gmsh.model.occ.addDisk(0, 0, 0, np.sqrt(1/np.pi), np.sqrt(1/np.pi)) # x, y, z, x-radius, y-radius
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)
metadata = {"quadrature_degree": 6}
# Define Measure Over Cells of Mesh
dx = ufl.Measure('dx', domain=domain, metadata=metadata)
# Define Measure Over Boundary of Mesh
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
# 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)
# Deformation Gradient Tensor
F = I + ufl.nabla_grad(u_) # Nabla_grad (Jared's Code) or grad?
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)
# Deformation Gradient Tensor
F_n = I_n + ufl.grad(u_n)
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):
small_strain = 0.5 * (ufl.nabla_grad(u) + ufl.transpose(ufl.nabla_grad(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)
class BGradExpression:
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_grad = BGradExpression(b_0=b_0, b_1=b_1)
b = fem.Function(M)
b.interpolate(b_grad.eval)
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 = J*D_a*ufl.dot(F_inv*F_inv_tra* ufl.grad(c_A), ufl.grad(w_A)) * dx
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
# term2_I = J*D_i*ufl.dot(F_inv*F_inv_tra* ufl.grad(c_I), ufl.grad(w_I)) * dx
term2_I = ufl.dot(D_I_eq * ufl.grad(c_I), ufl.grad(w_I)) * dx # Diffusion+Flux+Source?
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)
rhs_linear_vector.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# # 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
