Dear FEniCS community,
I’m currently working on the implementation of the Three Network Model for nonlinear viscoelasticity.
My code is structured following two similar FEniCS examples:
- Linear Viscoelasticity – COMET-FEniCS documentation
- Discussion on a viscoelastic problem – FEniCS Discourse
I’ve formulated a three-field variational problem to solve for the displacement field u and the two internal deformation gradient fields FvA and FvB.
However, I suspect there’s a problem with how I’m initializing or updating FvA and FvB, because when I print their values for debugging, both appear as zero matrices throughout the domain.
Has anyone encountered a similar issue or could suggest how to correctly initialize and evolve internal tensor fields like FvA and FvB in a mixed function space?
Thanks in advance for your help!
Here is the code:
from dolfinx import *
from matplotlib import pyplot as plt
import meshio
from dolfinx.io import gmshio
import pyvista
import numpy as np
from mpi4py import MPI
import math
import os
import sys
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import pyvista
from dolfinx import mesh, fem, plot, default_scalar_type, log
from basix.ufl import element, mixed_element
import ufl
from dolfinx.fem import dirichletbc, locate_dofs_topological, Function
from dolfinx.mesh import locate_entities_boundary
length = 10
width = 4
depth = 2
num_ele_along_depth = 4
size_ele = depth/num_ele_along_depth
n = [int(length/size_ele), int(width/size_ele), int(depth/size_ele)]
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([length, width, depth])], n, cell_type=mesh.CellType.hexahedron)
tdim = domain.topology.dim # Topological dimension
dim = domain.topology.dim - 1 # Facets are one dimension lower than the domain
######### Define the function space #########
#############################################
element_u = element("CG", domain.ufl_cell().cellname(), 2, shape=(tdim,)) # Displacement field
element_FvA = element("DG", domain.ufl_cell().cellname(), 0, shape=(tdim,tdim)) # FvA field
element_FvB = element("DG", domain.ufl_cell().cellname(), 0, shape=(tdim,tdim)) # FvB field
mixedelement = mixed_element([element_u, element_FvA, element_FvB]) # Mixed element for the function space
M = fem.functionspace(domain, mixedelement)
V, _ = M.sub(0).collapse() #Displacement
print("Number of DoFs in displacement:", V.dofmap.index_map.size_global) #DoFs 9328
V1, V1_to_M1 = M.sub(1).collapse() #Fv field A
print("Number of DoFs in Fv field A:", V1.dofmap.index_map.size_global)
V2, V2_to_M1 = M.sub(2).collapse() #Fv field B
print("Number of DoFs in Fv field B:", V2.dofmap.index_map.size_global)
# Coordinate X = 0
def boundary_x0(x):
return np.isclose(x[0], 0.0)
facets_x0 = locate_entities_boundary(domain, dim, boundary_x0)
dofs_x0 = locate_dofs_topological((M.sub(0), V), dim, facets_x0)
zero_displacement = Function(V)
zero_displacement.x.array[:] = 0.0
bc_x0 = dirichletbc(zero_displacement, dofs_x0, M.sub(0))
bcs = [bc_x0]
####### VARIATIONAL FORMULATION #######
#######################################
v = ufl.TestFunctions(M)
δu, δFvA, δFvB = v
m = fem.Function(M)
u, FvA, FvB = ufl.split(m)
print(f"u: {u.ufl_shape}")
print(f"FvA: {FvA.ufl_shape}")
print(f"FvB: {FvB.ufl_shape}")
# Saving the function at the previous time step
m_old = fem.Function(M) # to store all previous fields
u_old, Fv_oldA, Fv_oldB= ufl.split(m_old) # to access individual components
print(f"u_old: {u_old.ufl_shape}")
print(f"FvA_old: {Fv_oldA.ufl_shape}")
print(f"FvB_old: {Fv_oldB.ufl_shape}")
# ####### VARIATIONAL FORMULATION #######
# #######################################
#Traction vector
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
#Body force vector
B = fem.Constant(domain, default_scalar_type((0, -1, 0)))
#Material parameters
Bulk_Modulus=default_scalar_type(1000) #Considered as compressible materials
K=fem.Constant(domain, Bulk_Modulus)
#Material parameters for network A
MUA= fem.Constant(domain, default_scalar_type(15.58)) #Shear modulus network A
LAMBDAL= fem.Constant(domain, default_scalar_type(6.28)) #Lambda lock -- the same for all networks
TAUHATA= fem.Constant(domain, default_scalar_type(11.57)) #tauhat
MA=fem.Constant(domain, default_scalar_type(9.51))
#Material parameters for network B
MUBI=fem.Constant(domain, default_scalar_type(47.25)) #Initial Shear modulus network B
MUBF=fem.Constant(domain, default_scalar_type(1.32)) #Final Shear modulus network B
BETA=fem.Constant(domain, default_scalar_type(4.19)) #BETA
TAUHATB=fem.Constant(domain, default_scalar_type(39.56))
MB=fem.Constant(domain, default_scalar_type(14.67))
#Material parameters for network C
MUC=fem.Constant(domain, default_scalar_type(6.27))
State_value_n= fem.functionspace(domain, ("Lagrange", 1))
State_value_d= fem.functionspace(domain, ("Lagrange", 1))
V_viss = fem.functionspace(domain, ("Lagrange", 1))
def inv_Langevin(ED_chainstretch):
# Padè approximation for the inverse Langevin function
ED_limited = ufl.min_value(ED_chainstretch, 0.99)
L_inv = (ED_limited*(3-ED_limited**2))/(1-ED_limited**2)
chainstretch_expr = fem.Expression(ED_limited, V_viss.element.interpolation_points())
chainstretch_fn = fem.Function(V_viss)
chainstretch_fn.interpolate(chainstretch_expr)
print("Max chain stretch:", np.max(np.abs(chainstretch_fn.x.array)))
return L_inv
# Cauchy stress tensor for the Eight Chain model
def cauchy_stress_EC(FF_local, shear_modulus, lambdalock, K_local):
# Identity tensor
I_local = ufl.variable(ufl.Identity(3))
C_local = ufl.variable(ufl.transpose(FF_local) * FF_local)
J_local = ufl.max_value(ufl.variable(ufl.det(FF_local)), 1e-8)
b_local=ufl.variable(FF_local*ufl.transpose(FF_local))
b_local_dist = ufl.variable((J_local**(-2/3))*b_local)
Ib1_local = ufl.variable(ufl.tr(b_local))
Ib1_dist_local = ufl.variable(ufl.tr(b_local_dist))
# Invariants of deformation tensors
Ic1_local = ufl.variable(ufl.tr(C_local))
Ic1_dist_local = ufl.variable((J_local**(-2/3))*ufl.tr(C_local))
eff_dist_chainstrech= ufl.max_value(ufl.variable(ufl.sqrt(Ic1_dist_local/3)), 1e-8)
chainstretch_expr2 = fem.Expression(eff_dist_chainstrech, V_viss.element.interpolation_points())
chainstretch_fn2 = fem.Function(V_viss)
chainstretch_fn2.interpolate(chainstretch_expr2)
print("Max chain stretch 2:", np.max(np.abs(chainstretch_fn2.x.array)))
###########################################
# Calcola l'argomento di inv_Langevin
langevin_arg_n = ufl.variable(eff_dist_chainstrech / lambdalock)
# Calcola anche il valore della funzione inv_Langevin sull'argomento
langevin_val_n = ufl.variable(inv_Langevin(langevin_arg_n))
# Interpola entrambi in funzioni FEM per il debug numerico
langevin_arg_expr_n = fem.Expression(langevin_arg_n, State_value_n.element.interpolation_points(), domain.comm)
langevin_val_expr_n = fem.Expression(langevin_val_n, State_value_n.element.interpolation_points(), domain.comm)
langevin_arg_fn_n = fem.Function(State_value_n)
langevin_val_fn_n = fem.Function(State_value_n)
langevin_arg_fn_n.interpolate(langevin_arg_expr_n)
langevin_val_fn_n.interpolate(langevin_val_expr_n)
# Stampa i valori estremi
print("[DEBUG] Langevin argument NUMERATORE (eff_dist_chainstrech / lambdalock):")
print(" min =", np.min(langevin_arg_fn_n.x.array), "max =", np.max(langevin_arg_fn_n.x.array))
print("[DEBUG] NUMERATORE inv_Langevin(argument):")
print(" min =", np.min(langevin_val_fn_n.x.array), "max =", np.max(langevin_val_fn_n.x.array))
#################################
#################################
# Calcola l'argomento di inv_Langevin
langevin_arg_d = ufl.variable(1/ lambdalock)
# Calcola anche il valore della funzione inv_Langevin sull'argomento
langevin_val_d = ufl.variable(inv_Langevin(langevin_arg_d))
# Interpola entrambi in funzioni FEM per il debug numerico
langevin_arg_expr_d = fem.Expression(langevin_arg_d, State_value_d.element.interpolation_points(), domain.comm)
langevin_val_expr_d = fem.Expression(langevin_val_d, State_value_d.element.interpolation_points(), domain.comm)
langevin_arg_fn_d = fem.Function(State_value_d)
langevin_val_fn_d = fem.Function(State_value_d)
langevin_arg_fn_d.interpolate(langevin_arg_expr_d)
langevin_val_fn_d.interpolate(langevin_val_expr_d)
# Stampa i valori estremi
print("[DEBUG] Langevin argument DENOMINATORE (eff_dist_chainstrech / lambdalock):")
print(" min =", np.min(langevin_arg_fn_d.x.array), "max =", np.max(langevin_arg_fn_d.x.array))
print("[DEBUG] DENOMINATORE inv_Langevin(argument):")
print(" min =", np.min(langevin_val_fn_d.x.array), "max =", np.max(langevin_val_fn_d.x.array))
#################################
numeratore = inv_Langevin(eff_dist_chainstrech/lambdalock)
denominatore = (inv_Langevin(1/lambdalock))
part1=(shear_modulus/(J_local*eff_dist_chainstrech))*((numeratore)/(denominatore))
part2=part1*(b_local_dist-((1/3)*Ib1_dist_local*I_local))
return part2 + K_local*(J_local-1)*I_local
State_value = fem.functionspace(domain, ("Lagrange", 1))
State_tensor= fem.functionspace(domain, ("Lagrange", 1, (tdim,tdim)))
State_vector = fem.functionspace(domain, ("Lagrange", 1, (tdim,)))
#Inizializzazione FvA
initial_stateFviA = ufl.Identity(3)
initial_stateFviA_exp = fem.Expression(initial_stateFviA, V1.element.interpolation_points(), domain.comm)
initial_stateFviA_fn = fem.Function(V1)
initial_stateFviA_fn.interpolate(initial_stateFviA_exp)
m_old.x.array[V1_to_M1] = initial_stateFviA_fn.x.array
#m.x.array[V1_to_M1] = initial_stateFviA_fn.x.array
#Inizializzazione FvA
initial_stateFviB = ufl.Identity(3) # V2: spazio per FvB
initial_stateFviB_exp = fem.Expression(initial_stateFviB, V2.element.interpolation_points(), domain.comm)
initial_stateFviB_fn = fem.Function(V2)
initial_stateFviB_fn.interpolate(initial_stateFviB_exp)
m_old.x.array[V2_to_M1] = initial_stateFviB_fn.x.array
print(np.isnan(m_old.x.array).any())
#m.x.array[V2_to_M1] = initial_stateFviB_fn.x.array
state_muB_in = MUBI #Initial shear modulus of network B
state_muB_current = fem.Constant(domain, default_scalar_type(0)) #Current state variable for the shear modulus of network B
tau_C=fem.Constant(domain, default_scalar_type(0.001))
freq = 0.05 # Time step size
deltaT = freq #time step size
d = len(u) # Spatial dimension
I = ufl.variable(ufl.Identity(d)) # Identity tensor
F_new = ufl.variable(I + ufl.grad(u))
##F=F_el*F_vi
FelasticA_new = ufl.variable(F_new*ufl.inv(FvA))
FelasticB_new = ufl.variable(F_new*ufl.inv(FvB))
Cv_Anew = ufl.variable(ufl.transpose(FvA)*(FvA))
Cv_Aold = ufl.variable(ufl.transpose(Fv_oldA)*(Fv_oldA)) #CvA at the previous time step
Cv_Bnew = ufl.variable(ufl.transpose(FvB)*(FvB))
Cv_Bold = ufl.variable(ufl.transpose(Fv_oldB)*(Fv_oldB)) #CvB at the previous time step
#Network A
Cauchy_stress_A = cauchy_stress_EC(FelasticA_new, MUA, LAMBDAL, K)
dev_stressA =(Cauchy_stress_A-((1/3)*ufl.tr(Cauchy_stress_A)*I))
eps = 1e-8
tau_A = ufl.sqrt(ufl.tr(ufl.dot(dev_stressA,dev_stressA))) #Frobenius norm
hydr_pressure_A = -ufl.tr(Cauchy_stress_A)/3
RampA = (hydr_pressure_A+ufl.sqrt(hydr_pressure_A**2))/2
argA= ufl.conditional(ufl.gt(tau_A/(TAUHATA + RampA) - tau_C, 0.0),
tau_A/(TAUHATA + RampA) - tau_C,
0.0)
deltagammaA =(argA)**MA #Delta gamma A
f_A = 2*ufl.transpose(Fv_oldA)*deltagammaA*(dev_stressA/tau_A) * Fv_oldA
#Network B
state_muB_current = - BETA*(state_muB_in-MUBF)*(deltagammaA)*deltaT + state_muB_in
Cauchy_stress_B= cauchy_stress_EC(FelasticB_new, state_muB_current, LAMBDAL, K)
dev_stressB=(Cauchy_stress_B-((1/3)*ufl.tr(Cauchy_stress_B)*I))
tau_B=ufl.sqrt(ufl.tr(ufl.dot(dev_stressB,dev_stressB))) #Frobenius norm
hydr_pressure_B=-ufl.tr(Cauchy_stress_B)/3 #Ho levato qui il meno
RampB=(hydr_pressure_B+ufl.sqrt(hydr_pressure_A**2))/2
argB= ufl.conditional(ufl.gt(tau_B/(TAUHATB + RampB) - tau_C, 0.0),
tau_B/(TAUHATB + RampB) - tau_C,
0.0)
deltagammaB =(argB)**MB #Delta gamma B
f_B = 2*ufl.transpose(Fv_oldB)*deltagammaB*(dev_stressB/tau_B) * Fv_oldB
# Network C
CS_NetworkC = cauchy_stress_EC(F_new, MUC, LAMBDAL, K)
deltaCvA=(Cv_Anew-Cv_Aold)
deltaCvB=(Cv_Bnew-Cv_Bold)
def δe(test_function_local):
return ufl.sym(ufl.grad(test_function_local)) #Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata) # Elementarily volume
F_int_A= ufl.inner(Cauchy_stress_A, δe(δu))*dx + ufl.inner(deltaCvA,δFvA)*dx - deltaT*ufl.inner(f_A, δFvA)*dx
F_int_B= ufl.inner(Cauchy_stress_B, δe(δu))*dx + ufl.inner(deltaCvB,δFvB)*dx - deltaT*ufl.inner(f_B, δFvB)*dx
F_int_C= ufl.inner(CS_NetworkC, δe(δu))*dx
F_ext= - ufl.inner(δu, B) * dx - ufl.inner(δu, T) * ds
Variational= F_int_A + F_int_C + F_ext
log.set_log_level(log.LogLevel.INFO)
def bc_vec_moving(x, rate, time):
x = np.asarray(x)
dim = x.shape[0]
num_points = x.shape[1]
disp = np.zeros((dim, num_points), dtype=PETSc.ScalarType)
disp[0, :] = rate * time # Spostamento lungo x
return disp
def boundary_xL(x):
return np.isclose(x[0], length)
facets_xL = locate_entities_boundary(domain, dim, boundary_xL)
dofs_xL = locate_dofs_topological((M.sub(0), V), dim, facets_xL)
######## VISUALIZATION ########
#Function space for visualization
V_vis = fem.functionspace(domain, ("Lagrange", 1, (3,)))
TS = fem.functionspace(domain, ("DG", 0,(3, 3))) #True stress tensor
print(f"Number of DoFs in True Stress Tensor: {TS.dofmap.index_map.size_global}")
uh_vis_list = []
true_stress_vis_list =[]
with io.XDMFFile(domain.comm, "Viscoplasticity/results_1.xdmf", "w") as outfile:
outfile.write_mesh(domain) # write mesh once
t_end = 0.20 # Total simulation time
t = 0.10
step = 1
list_bc_m = []
while t <= t_end:
print(f"Time step {step}, time: {t:.2f}")
####
#### DEBUG NUMERICO u
arg_expr_n = fem.Expression(u, State_vector.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_vector)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] u:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("u new at first cell:", arg_fn_n.x.array[:3].reshape(3))
####
#### DEBUG NUMERICO F_new
arg_expr_n = fem.Expression(F_new, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] F_new:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("F new tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO det(F_new)
detF_expr = fem.Expression(ufl.det(F_new), State_value.element.interpolation_points(), domain.comm)
detF_fn = fem.Function(State_value)
detF_fn.interpolate(detF_expr)
print(f"[t = {t:.2f}] det(F): min={np.min(detF_fn.x.array)}, max={np.max(detF_fn.x.array)}")
####
#### DEBUG NUMERICO FelasticA_new
arg_expr_n = fem.Expression(FelasticA_new, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] FelasticA_new:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("Felastic new tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO FvA
arg_expr_n = fem.Expression(FvA, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] FvA:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("FvA tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO FvA_old
arg_expr_n = fem.Expression(Fv_oldA, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] FvA_old:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("FvA_old tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO Cv_Anew
arg_expr_n = fem.Expression(Cv_Anew, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] Cv_Anew :")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("Cv_Anew tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO Cv_Aold
arg_expr_n = fem.Expression(Cv_Aold, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] Cv_Aold:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("Cv_Aold tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO Cauchy_stress_A
arg_expr_n = fem.Expression(Cauchy_stress_A, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] Cauchy_stress_A:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print("Cauchy_stress_A tensor at first cell:", arg_fn_n.x.array[:9].reshape(3, 3))
####
#### DEBUG NUMERICO tau_A
arg_expr_n = fem.Expression(tau_A, State_value.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_value)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] tau_A:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
####
#### DEBUG NUMERICO hydr_pressure_A
arg_expr_n = fem.Expression(hydr_pressure_A, State_value.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_value)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] hydr_pressure_A:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
####
#### DEBUG NUMERICO RampA
arg_expr_n = fem.Expression(RampA, State_value.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_value)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] RampA:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
####
#### DEBUG NUMERICO argA
arg_expr_n = fem.Expression(argA, State_value.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_value)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] argA:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
####
#### DEBUG NUMERICO deltagammaA
arg_expr_n = fem.Expression(deltagammaA, State_value.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_value)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] deltagammaA:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
print(ufl.shape(deltagammaA))
####
#### DEBUG NUMERICO f_A
arg_expr_n2 = fem.Expression(f_A, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n2 = fem.Function(State_tensor)
arg_fn_n2.interpolate(arg_expr_n2)
print("[DEBUG] f_A:")
print(" min =", np.min(arg_fn_n2.x.array), "max =", np.max(arg_fn_n2.x.array))
####
#### DEBUG NUMERICO CS_NetworkC
arg_expr_n = fem.Expression(CS_NetworkC, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] CS_NetworkC:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
CS_TNM = Cauchy_stress_A + Cauchy_stress_B + CS_NetworkC
#####
# DEBUG NUMERICO deltaCvA
arg_expr_n = fem.Expression(deltaCvA, State_tensor.element.interpolation_points(), domain.comm)
arg_fn_n = fem.Function(State_tensor)
arg_fn_n.interpolate(arg_expr_n)
print("[DEBUG] deltaCvA:")
print(" min =", np.min(arg_fn_n.x.array), "max =", np.max(arg_fn_n.x.array))
##RISOLUZIONE
problem = NonlinearProblem(Variational, m, bcs)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8 #Sets the absolute tolerance (stopping criterion when the residual is small)
solver.rtol = 1e-8 #Sets the relative tolerance (stopping criterion based on relative change)
solver.max_it = 500 #Sets the maximum number of iterations
solver.convergence_criterion = "incremental" #Sets the convergence criterion to incremental
num_its, converged = solver.solve(m)
#num_its stores the number of Newton iterations required for convergence.
#converged is True if the solver successfully found a solution, otherwise False
if not converged:
print(f"Solver failed at step {step}. Iterations: {num_its}")
assert (converged) #Test if a condition returns True, if not, it raises an AssertionError
print(f"Time step {step}, Number of iterations {num_its}, Displacement: {displacement_rate*t}")
m.x.scatter_forward()
# Update the state variables
#initial_stateFviA = StateFviA_new
initial_stateFviA_expr = fem.Expression(FvA, State_tensor.element.interpolation_points(), domain.comm)
initialstateFviA = fem.Function(State_tensor)
initialstateFviA.interpolate(initial_stateFviA_expr)
Fv_oldA.x.array[:] = initialstateFviA.x.array[:]
assert not np.isnan(initial_stateFviA.x.array).any(), "NaN in StateFviA"
#initial_stateFviB = StateFviB_new
initial_stateFviB_expr = fem.Expression(FvB, State_tensor.element.interpolation_points(), domain.comm)
initialstateFviB = fem.Function(State_tensor)
initialstateFviB.interpolate(initial_stateFviB_expr)
Fv_oldB.x.array[:] = initialstateFviB.x.array[:]
assert not np.isnan(initial_stateFviB.x.array).any(), "NaN in StateFviB"
#state_muB_in = Stateshear_modulusB_T
state_muB_in_expr = fem.Expression(state_muB_current, State_value.element.interpolation_points(), domain.comm)
state_muB_in = fem.Function(State_value)
state_muB_in.interpolate(state_muB_in_expr)
assert not np.isnan(state_muB_in.x.array).any(), "NaN in muB"
#
TrueStress_expr = fem.Expression(CS_TNM, TS.element.interpolation_points())
stress_tensor_fn = fem.Function(TS)
stress_tensor_fn.interpolate(TrueStress_expr)
print("True stress tensor at first cell:", stress_tensor_fn.x.array[:9].reshape(3, 3))
true_stress_copy = fem.Function(TS)
true_stress_copy.x.array[:] = stress_tensor_fn.x.array[:]
true_stress_vis_list.append(true_stress_copy)
cs_array = stress_tensor_fn.x.array[:]
assert not np.isnan(cs_array).any(), "NaN in CS_TNM"
max_stress = np.max(np.abs(cs_array))
if max_stress > 1e6:
print(f"⚠️ Warning: High stress detected: {max_stress:.2e}")
print("Saved true stress solution")
# Update the previous solution
# Save a copy of the solution to keep history
u, FvA, FvB = m.split()
uh = m.sub(0).collapse()
uh_vis = fem.Function(V_vis)
uh_vis.interpolate(uh)
uh_copy = fem.Function(V_vis)
uh_copy.x.array[:] = uh_vis.x.array[:] # copia dei dati
uh_vis_list.append(uh_copy)
print("Saved displacement solution")
print("Max displacement norm:", np.max(np.abs(uh_vis.x.array)))
time_step = step
t += deltaT
step += 1
# WRITE XDMFFile
uh.name = "Displacement"
stress_tensor_fn.name = "TrueStress"
outfile.write_function(uh_vis, time_step)
outfile.write_function(stress_tensor_fn, time_step)
outfile.close()
print("Simulation completed")
Here is the output from my debug prints:
[DEBUG] u:
min = 0.0 max = 0.0
u new at first cell: [0. 0. 0.]
[DEBUG] F_new:
min = 0.0 max = 1.0
F new tensor at first cell: [[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[t = 0.10] det(F): min=1.0, max=1.0
[DEBUG] FelasticA_new:
min = nan max = nan
Felastic new tensor at first cell: [[nan nan nan]
[nan nan nan]
[nan nan nan]]
[DEBUG] FvA:
min = 0.0 max = 0.0
FvA tensor at first cell: [[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[DEBUG] FvA_old:
min = 0.0 max = 1.0
FvA_old tensor at first cell: [[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[DEBUG] Cv_Anew :
min = 0.0 max = 0.0
Cv_Anew tensor at first cell: [[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[DEBUG] Cv_Aold:
min = 0.0 max = 1.0
Cv_Aold tensor at first cell: [[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
[DEBUG] Cauchy_stress_A:
min = nan max = nan
Cauchy_stress_A tensor at first cell: [[nan nan nan]
[nan nan nan]
[nan nan nan]]
[DEBUG] tau_A:
min = nan max = nan
[DEBUG] hydr_pressure_A:
min = nan max = nan
[DEBUG] RampA:
min = nan max = nan
[DEBUG] argA:
min = 0.0 max = 0.0
[DEBUG] deltagammaA:
min = 0.0 max = 0.0
()
[DEBUG] f_A:
min = nan max = nan
[DEBUG] CS_NetworkC:
min = 0.0 max = 0.0
[DEBUG] deltaCvA:
min = -1.0 max = 0.0
[2025-05-30 18:09:48.376] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-05-30 18:09:48.377] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-05-30 18:09:52.047] [info] Column ghost size increased from 0 to 0
[2025-05-30 18:09:52.310] [info] PETSc Krylov solver starting to solve system.
[2025-05-30 18:09:54.668] [info] PETSc Krylov solver starting to solve system.
[2025-05-30 18:09:56.433] [info] Newton iteration 2: r (abs) = inf (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
[2025-05-30 18:09:56.700] [info] PETSc Krylov solver starting to solve system.