Non linear solver in FEniCSx. Comparison to SNES in legacy FEniCS

1

Hello FEniCSx family

In legacy FEniCS I am using the SNES solver as (definition of the solver and Non Linear Variational Problem):

snes_solver_parameters = {"nonlinear_solver": "snes",
                          "symmetric": True,
                          "snes_solver": {"maximum_iterations": 10,
                                          "report": True,
                                          "line_search": "bt",
                                          "linear_solver": "mumps",
                                          "method": "newtonls",
                                          "absolute_tolerance": 1e-9,
                                          "relative_tolerance": 1e-7,
                                          "error_on_nonconvergence": False}}

problem = NonlinearVariationalProblem(Res, w, bc, J=Jacobian)
solver_problem = NonlinearVariationalSolver(problem)
solver_problem.parameters.update(snes_solver_parameters)

I migrated my code in legacy FEniCS to FEniCSx. Here, I made it work with the following Newton solver:

problem = dolfinx.fem.petsc.NonlinearProblem(Res, u, bc, Jacobian)
solver_problem = dolfinx.nls.petsc.NewtonSolver(mesh.comm,problem)
solver_problem.convergence_criterion = "residual"
solver_problem.atol = 1e-9
solver_problem.rtol = 1e-7
solver_problem.max_it = 10
solver_problem.report = True
solver_problem.error_on_nonconvergence = False

However, I observe that the performance and the stability of the latter solver is much worse than the former SNES. Thus, my question:

What’s the closest solver/implementation in new FEniCSx to the SNES configuration I first report (my code in FEniCS legacy)?

Many thanks.

2

You can adapt dolfinx/test_nonlinear_assembler.py at 160ed13eb476df99944072aec70bd46a6fcb9bcf · FEniCS/dolfinx · GitHub

3

Many thanks.

With the second solver (NewtonSolver) in dolfinx, I also observe that if I launch a same simulation several times the behavior of the residual, and overall the convergence behavior, is not the same. Moreover, sometimes I get NaN, sometimes not.
This is unexpected for my because my codes in legacy FEniCS provided the same numerical behavior while solving, for all and any repetition of the simulation.

I cannot really understand the reason for such a behavior.

4

As you haven’t provided your code adapted to DOLFINx, I cannot comment on why that happens.

A Newton solver isn’t magical, it solves J du = - F(u_k), u_k+1=u_k+du,
As described in
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html

Therefore there is no reason in general for the result not being deterministic.

5

Many thanks for the comments, Dokken.
Here the full DOLFINx code:

from __future__ import print_function
import numpy as np
from numpy import array
import  dolfinx
import ufl
from mpi4py import MPI
import dolfinx.io
import math
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import petsc4py
import matplotlib.pyplot as plt

DOLFIN_EPS = 3.0e-16

name = "Output.xdmf" # Name of the output file
name1 = "Output.h5"

j,k,l,m = ufl.indices(4)

T = 6e-0
T_ramp=0.0000001
steps = 0

# Time parameters
tiempo = 0.0
tot_steps = 100;
dt = T / tot_steps

#Import mesh & subdomains from .msh file:
filename="Mesh"
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh(filename+'.msh', MPI.COMM_WORLD, gdim=3)

metadata = {"quadrature_degree": 2}
dx = ufl.Measure('dx')(domain=mesh, subdomain_data=cell_markers, metadata=metadata) #Integration measures for each subdomain, to integrate the constitutive equations on each subdomain.

DTT = dolfinx.fem.Constant(mesh,dt)


TT = dolfinx.fem.TensorFunctionSpace(mesh,('DG',0)) #To interpolate the stress tensors and stress fields at the end, not shape functions to solve the FEM problem.
VV = dolfinx.fem.VectorFunctionSpace(mesh,('DG',0)) #To eventually interpolate the nodal results for magnetic fields, not shape functions to solve the FEM problem.
SS = dolfinx.fem.FunctionSpace(mesh,('DG',1)) #To interpolate the phase field.

# FUNCTION SPACES
element = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
P1 = dolfinx.fem.FunctionSpace(mesh, element)
P3 = dolfinx.fem.FunctionSpace(mesh,('CG',1)) #Phase-field

# FUNCTIONS IN THE FUNCTION SPACE V (Function - Trial - Test: to produce the weak form later).
u = dolfinx.fem.Function(P1) #Unkowns.
u0 = dolfinx.fem.Function(P1) #Unkowns.
d_u = ufl.TrialFunction(P1)
v_u = ufl.TestFunction(P1)

d, d0, d_d, v_d = dolfinx.fem.Function(P3), dolfinx.fem.Function(P3), ufl.TrialFunction(P3), ufl.TestFunction(P3)
d0n = dolfinx.fem.Function(P3)
    
HHH0 = dolfinx.fem.Function(P3) #History energy.

#############################################################################
#BOUNDARY CONDITIONS
def bnd_top(x):
    tol=1e-4
    return np.isclose(x[1],30) #& np.greater(x[0],-3.5-tol) & np.less(x[0],6.5+tol)#\
        #& np.greater(x[1],30-tol) & np.less(x[1],30+tol) #np.isclose(x[1],30)
    #mesh.coordinates()[0, :].max()

def bnd_down(x):
    tol=1e-4
    return np.isclose(x[1],0)#&\
      #np.greater(x[0],-3.5-tol) & np.less(x[0],6.5+tol)

def bnd_top_pto(x):
    tol=1e-4
    return np.isclose(x[1],30) #& np.isclose(x[0],-3.5)

def bnd_down_pto(x):
    tol=1e-4
    return np.isclose(x[1],0) #& np.isclose(x[0],-3.5)

def bnd_downpoint(x):
    return np.close(x[1],0) & np.close(x[0],3)

def bnd_ref_pot(x):
    return np.isclose(x[1],-225) | np.isclose(x[1],255) | np.isclose(x[0],-150) | np.isclose(x[0],153)

def bnd_crack(x):
    tol=1e-5
    return np.greater(x[0],-3.5-tol) & np.less(x[0],-1.5) &\
      np.greater(x[1],14.975-tol) & np.less(x[1],15.025+tol)

def bnd_all(x):
    return True

def bnd_all2(x):
    return not bnd_crack(x)

def air(x):
    tol=1e-5
    return np.less(x[0],-3.5-tol)\
      | np.greater(x[0],6.5+tol)\
      | np.greater(x[1],30+tol)\
      | np.less(x[1],0-tol)

def bnd_horizontal_sym(x):
    tol = 0.001
    return np.greater(x[1],15-tol) & np.less(x[1],15+tol) &\
      np.greater(x[0],-3.5) & np.less(x[0],6.5)
      
def bnd_vertical_sym(x):
    return np.isclose(x[0],3)


##Multipoint constraint
def multip_const_up(x):
    tol=1e-5
    return np.greater(x[0],-150+tol) & np.less(x[0],153-tol)\
      & np.greater(x[1],30+tol) & np.less(x[1],255-tol)
      
def multip_const_down(x):
    tol=1e-5
    return np.greater(x[0],-150+tol) & np.less(x[0],153-tol)\
      & np.less(x[1],0-tol) & np.greater(x[1],-225+tol)

#######################################

displac_max = 170*2/2

class Mydisplacementup:
    def __init__(self):
        self.t = 0.0
        self.bias_up = 0.0
    def eval(self,x):
        #Add some spatial/temporal variation here:
        #x.shape[1]: number of columns.
        if self.t <= T_ramp:
            return np.full(x.shape[1], displac_max *(self.t-T_ramp)/(T-T_ramp))
        else:
            return np.full(x.shape[1], displac_max *(self.t-T_ramp)/(T-T_ramp))# + self.bias_up)
    def expres(self):
        if self.t <= T_ramp:
            return displac_max *(self.t-T_ramp)/(T-T_ramp)
        else:
            return displac_max *(self.t-T_ramp)/(T-T_ramp)# + self.bias_up

V_real = dolfinx.fem.FunctionSpace(mesh, ("CG",1))
load1_ = Mydisplacementup()   
load1_.t = 0.0
load1 = dolfinx.fem.Function(V_real)
load1.interpolate(load1_.eval)

class Mydisplacementdown:
    def __init__(self):
        self.t = 0.0
        self.bias_down = 0.0
    def eval(self,x):
        #Add some spatial/temporal variation here:
        #x.shape[1]: number of columns.
        if self.t <= T_ramp:
            return np.full(x.shape[1], 0.0)
        else:
            return np.full(x.shape[1], -displac_max *(self.t-T_ramp)/(T-T_ramp))# + self.bias_down)
    def expres(self):
        if self.t <= T_ramp:
            return 0.0
        else:
            return -displac_max *(self.t-T_ramp)/(T-T_ramp)# + self.bias_down
        

V_real = dolfinx.fem.FunctionSpace(mesh, ("CG",1))
load2_ = Mydisplacementdown()   
load2_.t = 0.0
load2 = dolfinx.fem.Function(V_real)
load2.interpolate(load2_.eval)

Taux = 0.2 #1.2*T_ramp

class Myprogressive:
    def __init__(self):
        self.t = 0.0
    def eval(self,x):
        #Add some spatial/temporal variation here:
        #x.shape[1]: number of columns.
        if self.t <= T_ramp:
            return np.full(x.shape[1], 0.0)
        else:
            if self.t <= Taux:
                #return np.full(x.shape[1], 0.0)
                return np.full(x.shape[1], (self.t-1.0*T_ramp)/(Taux-1.0*T_ramp)) #
            else:
                return np.full(x.shape[1], 1.0)

V_real = dolfinx.fem.FunctionSpace(mesh, ("CG",2))
progressive_ = Myprogressive()   
progressive_.t = 0.0
progressive = dolfinx.fem.Function(V_real)
progressive.interpolate(progressive_.eval)

#DIRICHLET BOUNDARY CONDITIONS:
#DISPLACEMENT BC:
#d.o.f.s in the region where the bcs are to be prescribed,
bnd_top_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_top)
#dirichletbc:
bc_u_up=dolfinx.fem.dirichletbc(load1,bnd_top_dofs01,P1.sub(1))

#d.o.f.s in the region where the bcs are to be prescribed,
bnd_top_dofs00 = dolfinx.fem.locate_dofs_geometrical((P1.sub(0),P1),bnd_top_pto)
#function with bc value
U0, submap0 = P1.sub(0).collapse()
u_top_vect00 = dolfinx.fem.Function(U0)
u_top_vect00.x.array[:] = 0
#dirichletbc:
bc_u_uph=dolfinx.fem.dirichletbc(u_top_vect00,bnd_top_dofs00,P1.sub(0))

#d.o.f.s in the region where the bcs are to be prescribed,
bnd_down_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_down)
#function with bc value
#U1=P1.sub(1)
#u_down_vect01 = dolfinx.fem.Function(U1)
#u_down_vect01.x.array[:] = 0
#u_down_vect01 = dolfinx.fem.Function(U1)
#dirichletbc:
bc_u_down = dolfinx.fem.dirichletbc(load2,bnd_down_dofs01,P1.sub(1)) #u_down_vect01

#d.o.f.s in the region where the bcs are to be prescribed,
bnd_down_dofs00 = dolfinx.fem.locate_dofs_geometrical((P1.sub(0),P1),bnd_down_pto)
#function with bc value
U0, submap0 = P1.sub(0).collapse()
u_down_vect00 = dolfinx.fem.Function(U0)
u_down_vect00.x.array[:] = 0
#dirichletbc:
bc_u_downh = dolfinx.fem.dirichletbc(u_down_vect00,bnd_down_dofs00,P1.sub(0))

#SYMMETRY DISPLACEMENT:
#d.o.f.s in the region where the bcs are to be prescribed,
bnd_symaxis_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_horizontal_sym)
#function with bc value
U1, submap1 = P1.sub(1).collapse()
u_sym_vect01 = dolfinx.fem.Function(U1)
u_sym_vect01.x.array[:] = 0
#dirichletbc:
bc_sym_h =dolfinx.fem.dirichletbc(u_sym_vect01,bnd_symaxis_dofs01,P1.sub(1))

#CRACK DAMAGE:
#d.o.f.s in the region where the bcs are to be prescribed,
bnd_crack_dofs = dolfinx.fem.locate_dofs_geometrical((P3,P3),bnd_crack)
#function with bc value
damage_bounairvectP3 = dolfinx.fem.Function(P3)
#dirichletbc:
bc_crack = dolfinx.fem.dirichletbc(progressive,bnd_crack_dofs,P3)

bcu = [bc_u_up,bc_u_down,bc_u_uph,bc_u_downh]
bcd = [bc_crack]

## VECTOR AND TENSOR FIELDS DERIVED FROM SCALAR AND VECTORIAL POTENTIAL FIELDS gi AND u (magnetic potential and displacement fields).
dd = len(u)
I = ufl.Identity(dd)
F = I + ufl.grad(u)
CG = ufl.dot(F.T,F)
BG = ufl.dot(F,F.T)

F0 = I + ufl.grad(u0)

def deactivate(d):
    return ufl.conditional(ufl.lt(d,0.1),1,0)

def d_positive(d):
    return ufl.conditional(ufl.lt(d,0.0),0,d)

def d_penalty_neg(d):
    return ufl.conditional(ufl.lt(d,0.0),1,0)
def d_penalty_pos(d):
    return ufl.conditional(ufl.gt(d,1.0),1,0)

## CONSTITUTIVE RELATIONS
i, j = ufl.indices(2)

########
G = [(6.2e-3)*5,1e-5]

n_ = [1.0, 1.0]
b_ = [1, 1]

Gc = [(2.45e-2)*2, 0]
l = 0.08

#Degradation function for the magnetic contributions, e.g., Pmag, Pmagr, Pcoup, Br,...
def mag_degrad_function(d):
    return (1-d)**2
    #return (1-d)**4
    #return 1-3*d**2+2*d**3
    #return (1-6*d_positive(d)**2+8*d_positive(d)**3-3*d_positive(d)**4)

def Piso(F,i):
    return G[i]*F*(1+b_[i]/n_[i]*(ufl.tr(ufl.dot(F.T,F))-3))**(n_[i]-1)

poisson=[0.45,0.4]
beta_=[2*poisson[0]/(1-2*poisson[0]), 2*poisson[1]/(1-2*poisson[1])]
beta=[((1-(1-d_positive(d))**2) *0.125*1.5 + ((1-d_positive(d))**2) *beta_[0]), beta_[1]]

kappa = [2*G[0]*(1+poisson[0])/3/(1-2*poisson[0]),2*G[1]*(1+poisson[1])/3/(1-2*poisson[1])]

strength_hs = [(130e-3)/2*2, 0]

def Pvol(F,i):
    return (G[i])*(-(ufl.conditional(ufl.lt(ufl.det(F),DOLFIN_EPS),DOLFIN_EPS,ufl.det(F))+DOLFIN_EPS)**(-beta[i])*ufl.inv(F).T)

def split_(F):
    return ufl.conditional(ufl.lt(ufl.det(F),1),1,1)

def psi(F,i):
    Psimech = G[i]/2/b_[i]*((1+b_[i]/n_[i]*(ufl.tr(CG)-3))**n_[i]-1)
    Psivol = G[i]/beta[i]*(ufl.det(F)**(-beta[i])-1)
    Psi = Psimech + Psivol*split_(F)
    return Psi

def psi_vol(F,i):
    Psivol = G[i]/beta[i]*(det(F)**(-beta[i])-1)
    return Psivol*split_(F)
def psi_mechiso(F,i):
    Psimech = G[i]/2*(ufl.tr(CG)-3)
    return Psimech

def HHH(F,HHH0):
    return ufl.conditional(ufl.lt(HHH0,psi(F,0)),psi(F,0),HHH0)

def Jacobiano(F):
    return ufl.det(F)

#### INTEGRATION MEASURES ####
tol = 1e-3
cdim = mesh.topology.dim
cells_bot = dolfinx.mesh.locate_entities(mesh, cdim, lambda x: (np.greater(x[1],22.5-tol) & np.less(x[1],24.5+tol)))
facet_tag = dolfinx.mesh.meshtags(mesh, 2, cells_bot, 1)
dxx = ufl.Measure("dx", domain=mesh, subdomain_data=facet_tag) #,
print('Area 1:',dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dxx(1))))

# VARIATIONAL PROBLEM - WEAK FORM
Res_u = ( ((1.0 - d0 )**2 + 0.003) * ufl.inner(Piso(F, 0), ufl.grad(v_u))+\
    ((1.0 - d0)**2 + 0.003)*ufl.inner(Pvol(F, 0), ufl.grad(v_u))) * dx

mu_visc = 0.1

Res_d = ufl.det(F)*(3/8 * Gc[0] * l * ufl.inner(ufl.dot(ufl.inv(F),ufl.grad(d)), ufl.dot(ufl.inv(F),ufl.grad(v_d))) + 3/8* v_d * Gc[0] / l + ufl.inner(d, v_d) * ( 2.0 * HHH(F,HHH0))\
        - 2.0 * HHH(F,HHH0) * v_d \
        - (1-d)*v_d*3**(5/4)*Gc[0]*kappa[0]/2/l/strength_hs[0]*ufl.tr(F-I)/(3+ufl.as_tensor((F-I)[j,k]*(F-I)[j,k]))**(5/4)\
        + 2*10**4*Gc[0]*d*v_d*d_penalty_neg(d) + 2*10**4*Gc[0]*(d-1)*v_d*d_penalty_pos(d)\
        + mu_visc/DTT*ufl.inner((d-d0n),v_d)) * dx

# Compute directional derivative about w in the direction of du (Jacobian)
Jacobian_u = ufl.derivative(Res_u, u, d_u)
Jacobian_d = ufl.derivative(Res_d, d, d_d)


# LAGRANGIAN FIELDS: magnetic induction field. To project the B-field on each subdomain according to its magnetic and mechanical properties.
def split_projectP(Piso,Pvol,F,d,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem = dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner(((1-d)**2)*(Piso(F,0)+Pvol(F,0)),v)*dx)
    f_tot = problem.solve()
    return f_tot


def split_projectCauchy(Piso,Pvol,F,d,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem =  dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner((1/det(F))*ufl.dot(((1-d)**2)*(Piso(F,0)+Pvol(F,0)),F.T),v)*dx) #Modify depending on the subdomains.
    f_tot = problem.solve()
    return f_tot


def split_projectPmech(Piso,Pvol,F,d,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem =  dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner(((1-d)**2)*(Piso(F,0)+Pvol(F,0)),v)*dx) #Modify depending on the subdomains.
    f_tot = problem.solve()
    return f_tot

def split_projectPmechiso(Piso,Pvol,F,d,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem =  dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner(((1-d)**2)*(Piso(F,0)),v)*dx) #Modify depending on the subdomains.
    f_tot = problem.solve()
    return f_tot

def split_projectPmechvol(Piso,Pvol,F,d,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem =  dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner(((1-d)**2)*(Pvol(F,0)),v)*dx) #Modify depending on the subdomains.
    f_tot = problem.solve()
    return f_tot


def split_projectHHH(F,H,HHH0,VSpace):
    u = ufl.TrialFunction(VSpace)
    v = ufl.TestFunction(VSpace)
    problem =  dolfinx.fem.petsc.LinearProblem(ufl.inner(u,v)*dx,\
        ufl.inner(HHH(F,H,HHH0,0),v)*dx) #Modify depending on the subdomains.
    f_tot = problem.solve()
    return f_tot

# SOLVER
# MECHANIC PROBLEM:
problem_u = dolfinx.fem.petsc.NonlinearProblem(Res_u, u, bcu, Jacobian_u)
solver_problem_u = dolfinx.nls.petsc.NewtonSolver(mesh.comm,problem_u)
solver_problem_u.convergence_criterion = "residual" #"residual" #"incremental"
solver_problem_u.atol = 1e-12
solver_problem_u.rtol = 1e-12
solver_problem_u.max_it = 10
solver_problem_u.report = True
solver_problem_u.error_on_nonconvergence = False
#PHASE FIELD PROBLEM
problem_d = dolfinx.fem.petsc.NonlinearProblem(Res_d, d, bcd, Jacobian_d)
solver_problem_d = dolfinx.nls.petsc.NewtonSolver(mesh.comm,problem_d)
solver_problem_d.max_it = 10
ksp = solver_problem_d.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
solver_problem_d1.error_on_nonconvergence = False

import os
if os.path.exists(name):
    os.remove(name)
    os.remove(name1)

# Save results into an .xdmf
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
file_results = dolfinx.io.XDMFFile(MPI.COMM_WORLD,name,'w')

#write the mesh and subdomains info
file_results.write_mesh(mesh)
cell_markers.name = "Phases"
file_results.write_meshtags(cell_markers)

#Subdomains volumetric fraction
print("Vol frac MRE: ", dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dx(1)))/dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dx(1)+1*dx(2))))
print("Vol frac Free space: ", dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dx(2)))/dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dx(1)+1*dx(2))))

#Point where the results will be writen in the txt files.

ttt=np.array([0])
fy=np.array([0]);fy2=np.array([0]);
fyy=np.array([0]);fyy1=np.array([0]); fyy1Cauchy=np.array([0]) #P_11=np.array([0])
v_d_base=np.array([0]);v_d_up=np.array([0]);displacem=np.array([0])
logstrain=np.array([0])
engstrain=np.array([0])

fyyPiola3=np.array([0])

point0 = np.array([[1.5,30,0]])
point5 = np.array([[1.5,0,0]])
point1 = np.array([1.5,15,0.0])
disp0 = np.zeros(3);
disp5 = np.zeros(3)

# Solve for each value using the previous solution as a starting point
while (tiempo < T):
    #Display report simulation
    from dolfinx import log
    log.set_log_level(log.LogLevel.INFO)
    ##Update time-dependent parameters:
    #dolfinx.fem Classes:
    load2_.t = tiempo
    load2.interpolate(load2_.eval)
    load1_.t = tiempo
    load1.interpolate(load1_.eval)
    progressive_.t = tiempo
    progressive.interpolate(progressive_.eval)

    DTT.value = dt
    
    if tiempo <= T_ramp:
        print('Applying mag. field')
    else:
        print('Mechanical load')
    
    print("dt: ",dt)
    print("Step: ",steps)
    print("Tiempo: ",tiempo,"/",T)

    iter = 0

    for aux in range(3):
        iter += 1
        print('Iteration (staggered scheme):', iter, 'Monolithic u')
        solver_problem_u.solve(u)
        print('Iteration (staggered scheme):', iter, 'Staggered d')
        solver_problem_d.solve(d)

        #Update history variables
        u0.vector.array[submap0] = u.vector.array[submap0]
        u0.vector.array[submap1] = u.vector.array[submap1]
        d0.vector.array[:] = d.vector.array
    
    HHH_expres = dolfinx.fem.Expression(HHH(F,HHH0), P3.element.interpolation_points());
    HHH0.interpolate(HHH_expres)


    log.set_log_level(log.LogLevel.OFF)
    
    # Write to .xdmf results file
    uVector=dolfinx.fem.Function(VV); uVector.name = "Displacement"; u_expr = dolfinx.fem.Expression(u, VV.element.interpolation_points()); uVector.interpolate(u_expr)
    file_results.write_function(uVector,tiempo)

    dField=dolfinx.fem.Function(SS); dField.name = "d - phase field"; d_expr = dolfinx.fem.Expression(d, SS.element.interpolation_points()); dField.interpolate(d_expr)
    file_results.write_function(dField,tiempo)

    d0Field=dolfinx.fem.Function(SS); d0Field.name = "d0 - phase field"; d0_expr = dolfinx.fem.Expression(d0, SS.element.interpolation_points()); d0Field.interpolate(d0_expr)
    file_results.write_function(d0Field,tiempo)

    HHH0Vector=dolfinx.fem.Function(SS); HHH0Vector.name = "History field"; HHH0_expr = dolfinx.fem.Expression(HHH0, SS.element.interpolation_points()); HHH0Vector.interpolate(HHH0_expr)
    file_results.write_function(HHH0Vector,tiempo)

    FTensor=dolfinx.fem.Function(TT); FTensor.name = "F"; F_expr = dolfinx.fem.Expression(F, TT.element.interpolation_points()); FTensor.interpolate(F_expr)
    file_results.write_function(FTensor,tiempo)

    PTensor = split_projectP(Piso,Pvol,F,d,TT); PTensor.name = "P"
    file_results.write_function(PTensor,tiempo)

    PmechTensor = split_projectPmech(Piso,Pvol,F,d,TT); PmechTensor.name = "Pmech"
    file_results.write_function(PmechTensor,tiempo)

    PmechisoTensor = split_projectPmechiso(Piso,Pvol,F,d,TT); PmechisoTensor.name = "Pmechiso"
    file_results.write_function(PmechisoTensor,tiempo)

    PmechvolTensor = split_projectPmechvol(Piso,Pvol,F,d,TT); PmechvolTensor.name = "Pmechvol"
    file_results.write_function(PmechvolTensor,tiempo)

    JacobianDet = dolfinx.fem.Function(SS); Jacobiandet_expres = dolfinx.fem.Expression(Jacobiano(F), SS.element.interpolation_points()); JacobianDet.name = "J"; JacobianDet.interpolate(Jacobiandet_expres);
    file_results.write_function(JacobianDet,tiempo)

    #REACTION FORCE
    fyyPiola3_t = dolfinx.fem.assemble_scalar(dolfinx.fem.form((1 / ufl.det(F)) * ufl.dot(ufl.dot((((1 - d) ** 2 + 0.0003) * (Piso(F, 0) + Pvol(F, 0))), F.T),ufl.as_vector([0, 1, 0]))[1] * dxx(1))) / 2 * 1

    ttt = np.append(ttt, tiempo)

    tree = dolfinx.geometry.BoundingBoxTree(mesh, 2)
    cell_candidates0 = dolfinx.geometry.compute_collisions(tree, point0)
    colliding_cells0 = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates0, point0)
    cells0 = colliding_cells0.links(0)
    disp0 = uVector.eval(point0, cells0[:1])
    v_d_up = np.append(v_d_up, disp0[1])
    
    cell_candidates5 = dolfinx.geometry.compute_collisions(tree, point5)
    colliding_cells5 = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates5, point5)
    cells5 = colliding_cells5.links(0)
    disp5 = uVector.eval(point5, cells5[:1])
    v_d_base = np.append(v_d_base, disp5[1])

    if tiempo < T_ramp:
        load1.bias_up = disp0[1]
        load2.bias_down = disp5[1]

        offset_up_ = disp0[1]
        offset_down_ = disp5[1]

        offset_stress = fyyPiola3_t

        print("offset up:",offset_up_)
        print("offset down:",offset_down_)

        displacem = np.append(displacem,0)
        logstrain = np.append(logstrain,0)
        engstrain = np.append(engstrain,0)

        fyyPiola3 = np.append(fyyPiola3,0) # N/mm

    else:
        displacem = np.append(displacem,(disp0[1])-(disp5[1]) + np.absolute(offset_up_) + np.absolute(offset_down_))
        logstrain = np.append(logstrain,ufl.ln((30+((disp0[1]) - (disp5[1]) + np.absolute(offset_up_) + np.absolute(offset_down_)))/(30+offset_up_-offset_down_)))
        engstrain = np.append(engstrain,(((disp0[1])-(disp5[1]) + np.absolute(offset_up_) + np.absolute(offset_down_))/(30+offset_up_-offset_down_)))

        fyyPiola3=np.append(fyyPiola3,fyyPiola3_t - offset_stress) # N/mm

    np.savetxt("Cauchy11_reaccion[kPa].txt", fyy1Cauchy)
    np.savetxt("displacement[mm].txt", displacem)
    np.savetxt("logstrain[-].txt", logstrain)
    np.savetxt("engstrain[-].txt", engstrain)
    np.savetxt("reaction_force[N_mm].txt", fyyPiola3)

    plt.figure()
    plt.plot(displacem[:],fyyPiola3[:])
    plt.grid();plt.xlabel("Displacement [mm]");plt.ylabel("Reaction force [N/mm]")
    plt.savefig('Reaction_force_integrando_dxx_arriba.png')
    plt.close()

    #ADAPTATIVE TIME STEP
    phase_field_diffFS = dolfinx.fem.FunctionSpace(mesh,("DG",1))
    phase_field_diff = dolfinx.fem.Function(phase_field_diffFS)
    phase_field_diff_express = dolfinx.fem.Expression(d-d0n, phase_field_diffFS.element.interpolation_points());
    phase_field_diff.interpolate(phase_field_diff_express)
    if phase_field_diff.vector.array.max() > 5e-4*100*3:
        dt = np.maximum(dt/1.2,0.08/15/5.0)
        tiempo +=dt
    else:
        dt = np.minimum(1.2*dt,0.08/5.0)
        tiempo += dt

    d0n.x.array[:] = d.x.array

    # Define boundary conditions
    steps += 1

    file_results.close()

print("Finished")
6

There are alot of things going on in this code, and for anyone to be able to help you, you need to do some debugging first to narrow down the cause of the issue.

  1. Do you have problems with the first solve (in time). If not, it means that your updating solutions scheme does not work. If you have problems after one step, i.e. First time you go through

It means that something in your setup further up is not correct, and you can distegard everything after this.

  1. check that your time dependent parameters are correct in the first loop prior to solving.

  2. you use a lot of projections. Is there a reason for not interpolating into a suitable space instead?

7

I would say there is something wrong with the time-updating Dirichlet BC, because they seem to no apply properly in the simulation , i.e., :

class Mydisplacementup:
    def __init__(self):
        self.t = 0.0
        self.bias_up = 0.0
    def eval(self,x):
        #Add some spatial/temporal variation here:
        #x.shape[1]: number of columns.
        return np.full(x.shape[1], displac_max  * self.t/T)

def bnd_top(x):
    tol=1e-4
    return np.isclose(x[1],30)

V_real = dolfinx.fem.FunctionSpace(mesh, ("CG",1))
load1_ = Mydisplacementup()   
load1_.t = 0.0
load1 = dolfinx.fem.Function(V_real)
load1.interpolate(load1_.eval)

#d.o.f.s in the region where the bcs are to be prescribed,
bnd_top_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_top)
#dirichletbc:
bc_u_up=dolfinx.fem.dirichletbc(load1,bnd_top_dofs01,P1.sub(1))

#And, to update the time-dependent Dirichtlet BC after the first iteration:
tiempo +=dt
load1_.t = tiempo
load1.interpolate(load1_.eval)

Note that I defined, in the same manner, class Mydisplacementup:, to prescribe a similar displacement Dirichtlet BC (i.e., load2) on another edge (2D simulation).
It looks weird: one Dirichlet BC seems to be working, but the other on the other edge, not. The mesh is symmetric, and it is about applying equal, but opposite sign, time-updating displacement BCs.
Perhaps issue with locate_dofs_geometrical?

The motivation to do it that way, i.e., using the functions split_project that I define previously, comes from the need to work with multiple subdomains. This simulation, however, was simplified to a single domain though.
Anyhow, there might be a better with to project the conjugate fields on a complete mesh with different subdomains according to specific properties of such subdomains.

8

[Update related to the not working application of the time-updating PCs:]

I have detected that they do work when using a CG Function Space (P1 in my code) with degree 1, but do not work with degree 2.
Moreover, I think it is related to the order of the Function Space used to interpolate the Expression defined within the Class object, i.e., V_real in my code.

9

This here does not make sense to me, what is the difference between the two?

I would still suggest trying to make a minimal code illustrating this issue, rather than the complex example, or a non complete example (make sure that the code is executable by anyone copy-pasting it).

10

Here a minimal code illustrating the issue:

from __future__ import print_function
import numpy as np
from numpy import array
import  dolfinx
import ufl
from mpi4py import MPI
import dolfinx.io
import math
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import petsc4py
import matplotlib.pyplot as plt

name = "Output.xdmf"

T = 2e-0  # final time

# Time parameters
tiempo = 0.0
tot_steps = 100;
dt = T / tot_steps

#Import mesh & subdomains from .msh file:
filename="Mesh"
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh(filename+'.msh', MPI.COMM_WORLD, gdim=3)
metadata = {"quadrature_degree": 2}
dx = ufl.Measure('dx')(domain=mesh, subdomain_data=cell_markers, metadata=metadata)

VV = dolfinx.fem.VectorFunctionSpace(mesh,('DG',0))
element = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
P1 = dolfinx.fem.FunctionSpace(mesh, element)
u = dolfinx.fem.Function(P1)
d_u = ufl.TrialFunction(P1)
v_u = ufl.TestFunction(P1)

def bnd_top(x):
    tol=1e-4
    return np.isclose(x[1],15)
def bnd_down(x):
    tol=1e-4
    return np.isclose(x[1],-15)

#DIRICHLET BOUNDARY CONDITIONS:
displac_max = 10
class Mydisplacementup:
    def __init__(self):
        self.t = 0.0
    def eval(self,x):
        return np.full(x.shape[1], displac_max *self.t/T)

V_real = dolfinx.fem.FunctionSpace(mesh, ("CG",2))
load1_ = Mydisplacementup()   
load1_.t = 0.0
load1 = dolfinx.fem.Function(V_real)
load1.interpolate(load1_.eval)

class Mydisplacementdown:
    def __init__(self):
        self.t = 0.0
    def eval(self,x):
        return np.full(x.shape[1], -displac_max *self.t/T)

load2_ = Mydisplacementdown()   
load2_.t = 0.0
load2 = dolfinx.fem.Function(V_real)
load2.interpolate(load2_.eval)

#DISPLACEMENT BC:
bnd_top_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_top)
bc_u_up = dolfinx.fem.dirichletbc(load1,bnd_top_dofs01,P1.sub(1))

bnd_down_dofs01 = dolfinx.fem.locate_dofs_geometrical((P1.sub(1),P1),bnd_down)
bc_u_down = dolfinx.fem.dirichletbc(load2,bnd_down_dofs01,P1.sub(1))

bcu = [bc_u_up,bc_u_down]

dd = len(u)
I = ufl.Identity(dd)
F = I + ufl.grad(u)
CG = ufl.dot(F.T,F)
## CONSTITUTIVE RELATIONS
G = [(6.2e-3)*10,0]
n_ = [1.0, 0]
b_ = [1, 0]
def Piso(F,i):
    return G[i]*F*(1+b_[i]/n_[i]*(ufl.tr(ufl.dot(F.T,F))-3))**(n_[i]-1)
poisson=[0.495,0]; beta=[2*poisson[0]/(1-2*poisson[0]), 2*poisson[1]/(1-2*poisson[1])]; kappa = [2*G[0]*(1+poisson[0])/3/(1-2*poisson[0]),2*G[1]*(1+poisson[1])/3/(1-2*poisson[1])]
def Pvol(F,i):
    return (G[i])*(-ufl.det(F)**(-beta[i])*ufl.inv(F).T)

Res_u = (ufl.inner(Piso(F, 0), ufl.grad(v_u)) + ufl.inner(Pvol(F, 0), ufl.grad(v_u))) * dx
Jacobian_u = ufl.derivative(Res_u, u, d_u)

problem_u = dolfinx.fem.petsc.NonlinearProblem(Res_u, u, bcu, Jacobian_u)
solver_problem_u = dolfinx.nls.petsc.NewtonSolver(mesh.comm,problem_u)
solver_problem_u.convergence_criterion = "residual" #"residual" #"incremental"
solver_problem_u.atol = 1e-12
solver_problem_u.rtol = 1e-12
solver_problem_u.max_it = 10
solver_problem_u.report = True
solver_problem_u.error_on_nonconvergence = False

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
file_results = dolfinx.io.XDMFFile(MPI.COMM_WORLD,name,'w')
file_results.write_mesh(mesh); cell_markers.name = "Phases"; file_results.write_meshtags(cell_markers)

while (tiempo < T):
    from dolfinx import log
    log.set_log_level(log.LogLevel.INFO)

    load1_.t = tiempo
    load1.interpolate(load1_.eval)
    load2_.t = tiempo
    load2.interpolate(load2_.eval)
            
    solver_problem_u.solve(u)

    uVector=dolfinx.fem.Function(VV); uVector.name = "Displacement"; u_expr = dolfinx.fem.Expression(u, VV.element.interpolation_points()); uVector.interpolate(u_expr)
    file_results.write_function(uVector,tiempo)

    tiempo += dt
    file_results.close()

This does not work (i.e., correct application of the BCs) as is, see that there is written element = ufl.VectorElement("CG", mesh.ufl_cell(), 2).

However, it does apply correctly the time-updating displacement BCs when using element = ufl.VectorElement("CG", mesh.ufl_cell(), 1).

Note: the mesh is loaded from a .msh file. I do not get to attach it here.

11

Well, you are using massive underintegration

when using a second order element:

from ufl.algorithms.compute_form_data \
    import estimate_total_polynomial_degree
import ufl.algorithms

# ...

print(estimate_total_polynomial_degree(ufl.algorithms.expand_derivatives(Jacobian_u)))
print(estimate_total_polynomial_degree(ufl.algorithms.expand_derivatives(Res_u)))

returning

17 
11

for P2 and

2
2

for P1.

Secondly, your boundary condition does not make sense

As the second argument should be the V_real space, which I recommend you getting from


V_real, _ = P1.sub(0).collapse()
# Define load 1 and load 2
# ...

# DISPLACEMENT BC:
bnd_top_dofs01 = dolfinx.fem.locate_dofs_geometrical(
    (P1.sub(1), V_real), bnd_top)
bc_u_up = dolfinx.fem.dirichletbc(load1, bnd_top_dofs01)

bnd_down_dofs01 = dolfinx.fem.locate_dofs_geometrical(
    (P1.sub(1), V_real), bnd_down)
bc_u_down = dolfinx.fem.dirichletbc(load2, bnd_down_dofs01)

bcu = [bc_u_up, bc_u_down]

I am not able to give you any further guidance before you:

  1. Use a built in mesh (Say a square or a cube)
  2. Remove all code not needed to reproduce the result (you should be able to run a single solve (in time), instead of solving for multiple steps
  3. Remove all unused imports.
12

Now it works. Apparently it was the incorrect definition of de dofs. Using V_real as argument of locate_dofs_geometrical seems to solve it. Many thanks.