Form compiles differently inside vs outside time loop [3D elasticity]

[Deepak_Dhariwal]

Hi,

I am trying to solve a 3D elasticity problem with several additional field variables for eigenstrains (a simpler MWE enclosed at the end). Geometry is hemisphere with tetrahedral meshing, compressed uniaxially against a rigid wall.

Results of MWE when form is compiled in the loop (CONSISTENT)

Results of MWE when form is compiled outside the loop and boundary conditions updated after each iteration (ERROR)

I am enclosing MWE with just two eigenstrains. I want to draw your attention to the latter part of code which I am modifying to compile the forms outside of loop.

Common block of code:

from numpy import infty, iinfo, int32, argwhere, min, max, array

from dolfin import *


# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

parameters.parse()



## Geometrical parameters
R = 100  # nm
mesh_size = 8
compression = R/50.0
time_step   = 5e-2


## material parameters
AM_surface_energy = 0.20     # J/m^2 = Pa.m = Pa'.nm (no change in magnitude in new scaling)
AM_length_scale   = 15.0     # nm
MM_surface_energy = 0.00     # J/m^2 = Pa.m = Pa'.nm (no change in magnitude in new scaling)
MM_length_scale   = 0.00     # nm

phi               = 5.0e-3   # 1 MPa = 1e-3 * Pa'

contact_penalty   = 100
phase_penalty     = 100

mob               = 10e-3



mesh_folder = "./meshes/" + "tetrahedral/"
mesh_file = "mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".xdmf"


mesh = Mesh()

#db_mesh_folder = "/Users/dhariwal/Library/CloudStorage/OneDrive-VirginiaTech/Zirconia/Particle_SIMULATION/Code/single_variant/meshes/"
#mesh_file = "mesh3D_R_" + str(int(self.R)) + "_size_" + str(int(mesh_size)) + ".xdmf"

with XDMFFile(mesh_folder + mesh_file) as infile:
    infile.read(mesh)


meshSIZE = 0.5*(mesh.hmax() + mesh.hmin())


# basis for 'vector' displacement
P1 = VectorElement("Lagrange", mesh.ufl_cell(), 1)

# basis for 'scalar' monovariant phase field
P2 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

W1 = FunctionSpace(mesh, P1)
W2 = FunctionSpace(mesh, P2)



# lattice parameters
alpha = 1.0243 
beta = -0.0247
gamma = 0.9563
delta = 0.058


e_v1 = array([[gamma - 1, beta, beta],
        [beta, alpha - 1, delta],
        [beta, delta, alpha - 1]
        ])

e_v2 = array([[gamma - 1, - beta, - beta],
        [- beta, alpha - 1, delta],
        [- beta, delta, alpha - 1]
        ])


shape_strains_list = [as_tensor(e_v1), as_tensor(e_v2)]

num_variants =  2

# single variant
element = MixedElement([P1, P2, P2])


# create the space of solutions over the domain 
W = FunctionSpace(mesh, element)

# a generic function from W
w  = Function(W)
dw = TrialFunction(W)
wx = TestFunction(W)



compliance_cubic = array([
    [169.0, 138.0, 138.0,  0.0,  0.0,  0.0],
    [138.0, 169.0, 138.0,  0.0,  0.0,  0.0],
    [138.0, 138.0, 169.0,  0.0,  0.0,  0.0],
    [  0.0,   0.0,   0.0, 40.0,  0.0,  0.0],
    [  0.0,   0.0,   0.0,  0.0, 40.0,  0.0],
    [  0.0,   0.0,   0.0,  0.0,  0.0, 40.0]
])

compliance_ortho = array([
    [223.0, 129.0,  99.0,   0.0,    0.0,  27.0],
    [129.0, 241.0, 125.0,   0.0,    0.0,  -9.0],
    [ 99.0, 125.0, 200.0,   0.0,    0.0,   4.0],
    [  0.0,   0.0,   0.0,  76.0,   -4.0,   0.0],
    [  0.0,   0.0,   0.0,  -4.0,   45.0,   0.0],
    [ 27.0,  -9.0,   4.0,   0.0,    0.0,   90.0]
])


Compliance_cubic = as_tensor(compliance_cubic)
Compliance_ortho = as_tensor(compliance_ortho)



# Definition of normal vectors to the elastic body and to the rigid half-space
nRIGID = Constant((0.0, 0.0, 1.0))
n = FacetNormal(mesh)


##################################
# Phase field definitions
##################################


def trans_strain(c):
    strain = [c[i] * shape_strains_list[i] for i in range(num_variants)]
    return sum(strain) 


def mixture_compliance(c):
    return ((1 - sum(c)) * Compliance_cubic) + (sum(c) * Compliance_ortho)


##################################
# Displacement field definitions
##################################


# matrix-to-vector transform
def strain3voigt(e):
    return as_vector(
        [e[0, 0], 
        e[1, 1], 
        e[2, 2], 
        e[1, 2] * 2, 
        e[0, 2] * 2, 
        e[0, 1] * 2]
                    )

# vector-to-matrix transform
def voigt3stress(s):
    return as_tensor(
        [[s[0], s[5], s[4]], 
        [s[5], s[1], s[3]], 
        [s[4], s[3], s[2]]]
                    )


def epsilon(u, strain_approx=None):

    if strain_approx is None:
        e = sym(grad(u))
        return e 
    
    elif strain_approx =='lagrangian':
        I = Identity(3)
        Fgrad = I + grad(u)
        C = Fgrad.T * Fgrad
        E = 0.5*(C - I)
        return E
    
    elif strain_approx =='hencky_order_22':
        I = Identity(3)
        Fgrad = I + grad(u)
        C = Fgrad.T * Fgrad
        a = C*C - I
        b = (C*C) + 4*C + I
        H = 1.5 * a * inv(b)
        return H

def epsilon_elastic(u, c):
    return epsilon(u) - trans_strain(c)

def Vstrain(u, c):
    return strain3voigt(epsilon_elastic(u, c))

def Vsigma(u, c):
    return dot(mixture_compliance(c), Vstrain(u, c))

def sigma(u, c):
    return voigt3stress(Vsigma(u, c))

def dilatation(c):
    return det(trans_strain(c))


def gap(u): # Definition of gap function
    x = SpatialCoordinate(mesh)
    return x[2] + u[2]

def maculay(x): # Definition of Maculay bracket
    return (x+abs(x))/2

def pN(u, c):
    return dot(dot(nRIGID, sigma(u, c)), nRIGID) 

def convert_to_integer(a):
    "Convert to a 32-bit int. Raise exception if this will cause an overflow"
    if abs(a) <= iinfo(int32).max:
        return int32(a)
    else:
        raise OverflowError("Cannot safely convert to a 32-bit int.")


def psi(w):

    u = split(w)[0]
    c = split(w)[1:]


    def elastic(u, c):
        # elastic_energy = 0.5 * dilatation(c) * inner(sigma(u, c), epsilon_elastic(u, c))
        
        # approximate energy without volume dilatation (det F = 1)
        elastic_energy = 0.5 * inner(sigma(u, c), epsilon_elastic(u, c))
        
        return elastic_energy

    def interface(c):

        energy_prefactor = (4.0/DOLFIN_PI) * (AM_surface_energy/AM_length_scale)
        AM_interface = energy_prefactor * inner(1 - sum(c), sum(c))
        energy_cost = AM_interface


        if MM_length_scale > 0.0:
            MM_interface = (4.0/DOLFIN_PI) * (MM_surface_energy/MM_length_scale)
            energy_cost = AM_interface + MM_interface
        
        return energy_cost, energy_prefactor

    def chemical(c):
        chemical_energy = phi*sum(c)
        return chemical_energy


    def contact_pen(u):
        pen = (0.5 * contact_penalty *inner(maculay(-gap(u)),maculay(-gap(u))))
        return pen

    def phase_pen(c):   
        penp1 = [inner(maculay(c[i] - 1), maculay(c[i] - 1)) for i in range(num_variants)]
        penp0 = [inner(maculay(-c[i]), maculay(-c[i])) for i in range(num_variants)]

        penp = sum(penp0) + sum(penp1)
                    
        return 0.5 * phase_penalty * penp
    

    def dissipation(u, c):

        prefactor = interface(c)[1]

        diss = 0.0

        for i in range(num_variants):

            X_AM = inner(sigma(u, c), shape_strains_list[i]) \
                    - 0.5 * inner(voigt3stress(dot((Compliance_ortho - Compliance_cubic), Vstrain(u, c))), epsilon_elastic(u, c)) \
                    - phi - prefactor*(1.0 - 2*sum(c))

            cdot = mob * X_AM
            
            diss = diss + (X_AM*cdot)

        return diss


    # return elastic(u, c), contact_pen(u), phase_pen(c)
    return elastic(u, c) + interface(c)[0] + chemical(c), contact_pen(u), phase_pen(c), dissipation(u, c)



def sphere_boundaries(step):

    # Spatial coordinates of mesh
    x = SpatialCoordinate(mesh)

    # Definition of tolerance
    tol = DOLFIN_EPS # FEniCS tolerance - necessary for "boundary_markers" and the correct definition of the boundaries 


    bT = CompiledSubDomain('near(x[2], R, tol) && on_boundary', R = R, tol = tol)

    bP = CompiledSubDomain('sqrt(x[0]*x[0] + x[1]*x[1]) < meshSIZE \
                            && near(x[2], R, tol) \
                            && on_boundary', meshSIZE = meshSIZE, R=R, tol = tol)

    # Definition of contact zones
    bC = CompiledSubDomain('on_boundary && x[2] < compression + tol', R = R, tol=tol, compression = step) # Contact search - contact part of the boundary

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    

    bT.mark(boundaries, 2)      # top
    bP.mark(boundaries, 4)      # lateral_penalty (for stability)
    bC.mark(boundaries, 6)      # contact_zone

    print("Nodes at contact zone: ", boundaries.array()[boundaries.array() == 6])



    dom_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
    dom_markers.set_all(1)


    # Displacement boundary conditions
    bc1 = DirichletBC(W.sub(0).sub(2), Constant((-step)), boundaries, 2)
    bc2 = DirichletBC(W.sub(0).sub(0), Constant((0)), boundaries, 4)
    bc3 = DirichletBC(W.sub(0).sub(1), Constant((0)), boundaries, 4)


    bcs = [bc1,bc2,bc3]


    dx = Measure('dx', subdomain_data=dom_markers)
    ds = Measure('ds', subdomain_data=boundaries)

    return dx, ds, boundaries, bcs

When form is compiled inside time loop:

while t_curr <= time_end:

    dx, ds, boundaries, bcs = sphere_boundaries(d_step)

    print("Current displacement at boundary = ", d_step, " nm and time_step = ", time_step)



    # Total potential energy at the current step (linear elasticity, gradient-interface)
    PiT = psi(w)[0]*dx + psi(w)[1]*ds(6) + psi(w)[2]*dx + psi(w)[3]*dx

    # Total potential energy of previous time-step (from known solution)
    Pi0 = psi(w0)[0]*dx

    # Incremental energy (increment potential energy + increment in dissiaption)
    Pi = PiT - Pi0

    # Compute first variation of Pi (directional derivative about u in the direction of v)
    F = derivative(Pi, w, wx)

    # Compute Jacobian of F for iterative "newton_solver"
    J = derivative(F, w, dw)



    for bc in bcs:
        bc.apply(w0.vector())


    w.vector()[:] = w0.vector()[:]


    problem = NonlinearVariationalProblem(F, w, bcs, J=J)
    solver  = NonlinearVariationalSolver(problem)


    snes_solver_parameters = {"nonlinear_solver": "snes",              
                                "snes_solver": {"linear_solver": 'lu', #"lu",
                                                 "preconditioner": "pcmg", 
                                                 "line_search": "bt",
                                                 "maximum_iterations": 50,
                                                 "absolute_tolerance": 1E-3,
                                                 "relative_tolerance": 1E-3,
                                                 "report": True,
                                                 #"method":'vinewtonrsls',
                                                 "error_on_nonconvergence": False,
                                                 }}
        
    solver.parameters.update(snes_solver_parameters)

    iter, converged = solver.solve()

    if converged:

        #time_step = 1.5 * time_step

        t_curr = t_curr + time_step

        #u_converged = _u
        w0.assign(w)

    else:
        # new timestamp is set, go back to the previous
        t_curr = t_curr - time_step

        # update the step-size
        time_step = time_step*(0.5)

        # modified current time 
        t_curr = t_curr + time_step


    d_step = compression * t_curr        



_c = w.split(deepcopy=True)[1:]

mart = project(sum(_c), W2)

from vedo.dolfin import plot
plot(mart)

When form is compiled outside time loop and BCs updated after each iteration:

dx, ds, boundaries, bcs = sphere_boundaries(d_step)

for bc in bcs:
    bc.apply(w0.vector())

w.vector()[:] = w0.vector()[:]


# Total potential energy at the current step (linear elasticity, gradient-interface)
PiT = psi(w)[0]*dx + psi(w)[1]*ds(6) + psi(w)[2]*dx + psi(w)[3]*dx

# Total potential energy of previous time-step (from known solution)
Pi0 = psi(w0)[0]*dx

# Incremental energy (increment potential energy + increment in dissiaption)
Pi = PiT - Pi0

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, w, wx)

# Compute Jacobian of F for iterative "newton_solver"
J = derivative(F, w, dw)




while t_curr <= time_end:

    print("Current displacement at boundary = ", d_step, " nm and time_step = ", time_step)

    problem = NonlinearVariationalProblem(F, w, bcs, J=J)
    solver  = NonlinearVariationalSolver(problem)

    solver.parameters.update(snes_solver_parameters)


    iter, converged = solver.solve()

    if converged:

        #time_step = 1.5 * time_step

        t_curr = t_curr + time_step

        #u_converged = _u
        w0.assign(w)

    else:
        # new timestamp is set, go back to the previous
        t_curr = t_curr - time_step

        # update the step-size
        time_step = time_step*(0.5)

        # modified current time 
        t_curr = t_curr + time_step


    d_step = compression * t_curr

    dx, ds, boundaries, bcs = sphere_boundaries(d_step) 



_c = w.split(deepcopy=True)[1:]

mart = project(sum(_c), W2)

from vedo.dolfin import plot
plot(mart)

I am not able to understand from where such discrepancy in results might be arising?

Thanks,
Deepak

[dokken]

The code you have supplied is way too big for most people to parse. Please try to make your example a minimal example, removing all things that are not required to reproduce the problem.

I would guess the issue comes down to you redefining dx, ds etc inside the loop.

@Deepak_Dhariwal

  d_step = compression * t_curr

    dx, ds, boundaries, bcs = sphere_boundaries(d_step) 

I would strongly suggest avoiding this, and rather updating the values of the once defined dom_markers, boundaries`, instead of redefining them.

[Deepak_Dhariwal]

Thank you for your suggestions, @dokken

I have shortened the code, removed extra comments and long whitespaces. This code is enclosed at the end of this post. I couldn’t change the energy functional because it is needed for initiation of transformation and comparison.

Redefining boundaries: This anisotropic sphere is compressed against the wall, I am using penalty regularization (gap > 0, pressure < 0). With every incremental boundary displacement (or, compression) the contact surface changes. Surface integral of contact penalty in total potential energy expression needs an updated boundary surface. See this:

psi(w)[0]*dx + psi(w)[1]*ds(6) + psi(w)[2]*dx + psi(w)[3]*dx

Boundary marker (6) is for the contact surface.

I have made changes per your suggestion NOT to redefine boundaries. The result is still the same:

(breaking into two parts for relatively easy parsing)

Part 1. material parameters, geometry, function definitions)

from time import time
from numpy import infty, iinfo, int32, argwhere, min, max, array

from dolfin import *


# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

parameters.parse()



## Geometrical parameters
R = 100  # nm
mesh_size = 8
compression = R/50.0
time_step   = 5e-2


## material parameters
AM_surface_energy = 0.20     # J/m^2 = Pa.m = Pa'.nm (no change in magnitude in new scaling)
AM_length_scale   = 15.0     # nm
MM_surface_energy = 0.00     # J/m^2 = Pa.m = Pa'.nm (no change in magnitude in new scaling)
MM_length_scale   = 0.00     # nm

phi               = 5.0e-3   # 1 MPa = 1e-3 * Pa'

contact_penalty   = 100
phase_penalty     = 100

mob               = 10e-3


mesh_folder = "./meshes/" + "tetrahedral/"
mesh_file = "mesh3D_R_" + str(int(R)) + "_size_" + str(int(mesh_size)) + ".xdmf"

mesh = Mesh()

with XDMFFile(mesh_folder + mesh_file) as infile:
    infile.read(mesh)

meshSIZE = 0.5*(mesh.hmax() + mesh.hmin())


# basis for 'vector' displacement
P1 = VectorElement("Lagrange", mesh.ufl_cell(), 1)

# basis for 'scalar' monovariant phase field
P2 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

W1 = FunctionSpace(mesh, P1)
W2 = FunctionSpace(mesh, P2)


# lattice parameters
alpha, beta, gamma, delta = 1.0243, -0.0247, 0.9563, 0.058 


e_v1 = array([[gamma - 1, beta, beta],
        [beta, alpha - 1, delta],
        [beta, delta, alpha - 1]
        ])

e_v2 = array([[gamma - 1, - beta, - beta],
        [- beta, alpha - 1, delta],
        [- beta, delta, alpha - 1]
        ])


shape_strains_list = [as_tensor(e_v1), as_tensor(e_v2)]
num_variants =  2

# single variant
element = MixedElement([P1, P2, P2])

# create the space of solutions over the domain 
W = FunctionSpace(mesh, element)

# a generic function from W
w  = Function(W)
dw = TrialFunction(W)
wx = TestFunction(W)


Compliance_cubic = as_tensor(array([
    [169.0, 138.0, 138.0,  0.0,  0.0,  0.0],
    [138.0, 169.0, 138.0,  0.0,  0.0,  0.0],
    [138.0, 138.0, 169.0,  0.0,  0.0,  0.0],
    [  0.0,   0.0,   0.0, 40.0,  0.0,  0.0],
    [  0.0,   0.0,   0.0,  0.0, 40.0,  0.0],
    [  0.0,   0.0,   0.0,  0.0,  0.0, 40.0]
]))

Compliance_ortho = as_tensor(array([
    [223.0, 129.0,  99.0,   0.0,    0.0,  27.0],
    [129.0, 241.0, 125.0,   0.0,    0.0,  -9.0],
    [ 99.0, 125.0, 200.0,   0.0,    0.0,   4.0],
    [  0.0,   0.0,   0.0,  76.0,   -4.0,   0.0],
    [  0.0,   0.0,   0.0,  -4.0,   45.0,   0.0],
    [ 27.0,  -9.0,   4.0,   0.0,    0.0,   90.0]
]))


# Definition of normal vectors to the elastic body and to the rigid half-space
nRIGID = Constant((0.0, 0.0, 1.0))
n = FacetNormal(mesh)


### Phase field definitions

def trans_strain(c):
    strain = [c[i] * shape_strains_list[i] for i in range(num_variants)]
    return sum(strain) 

def mixture_compliance(c):
    return ((1 - sum(c)) * Compliance_cubic) + (sum(c) * Compliance_ortho)



### Displacement field definitions

def strain3voigt(e):
    return as_vector(
        [e[0, 0], 
        e[1, 1], 
        e[2, 2], 
        e[1, 2] * 2, 
        e[0, 2] * 2, 
        e[0, 1] * 2]
                    )

def voigt3stress(s):
    return as_tensor(
        [[s[0], s[5], s[4]], 
        [s[5], s[1], s[3]], 
        [s[4], s[3], s[2]]]
                    )

def epsilon(u):
    e = sym(grad(u))
    return e 
    
def epsilon_elastic(u, c):
    return epsilon(u) - trans_strain(c)

def Vstrain(u, c):
    return strain3voigt(epsilon_elastic(u, c))

def Vsigma(u, c):
    return dot(mixture_compliance(c), Vstrain(u, c))

def sigma(u, c):
    return voigt3stress(Vsigma(u, c))

def dilatation(c):
    return det(trans_strain(c))



def gap(u): # Definition of gap function
    x = SpatialCoordinate(mesh)
    return x[2] + u[2]

def maculay(x): # Definition of Maculay bracket
    return (x+abs(x))/2

def convert_to_integer(a):
    "Convert to a 32-bit int. Raise exception if this will cause an overflow"
    if abs(a) <= iinfo(int32).max:
        return int32(a)
    else:
        raise OverflowError("Cannot safely convert to a 32-bit int.")



##################################
# Total potential energy definition
##################################

def psi(w):

    u = split(w)[0]
    c = split(w)[1:]


    def elastic(u, c):
        # elastic_energy = 0.5 * dilatation(c) * inner(sigma(u, c), epsilon_elastic(u, c))
        
        # approximate energy without volume dilatation (det F = 1)
        elastic_energy = 0.5 * inner(sigma(u, c), epsilon_elastic(u, c))
        
        return elastic_energy

    def interface(c):

        energy_prefactor = (4.0/DOLFIN_PI) * (AM_surface_energy/AM_length_scale)
        AM_interface = energy_prefactor * inner(1 - sum(c), sum(c))
        energy_cost = AM_interface


        if MM_length_scale > 0.0:
            MM_interface = (4.0/DOLFIN_PI) * (MM_surface_energy/MM_length_scale)
            energy_cost = AM_interface + MM_interface
        
        return energy_cost, energy_prefactor

    def chemical(c):
        chemical_energy = phi*sum(c)
        return chemical_energy


    def contact_pen(u):
        pen = (0.5 * contact_penalty *inner(maculay(-gap(u)),maculay(-gap(u))))
        return pen

    def phase_pen(c):   
        penp1 = [inner(maculay(c[i] - 1), maculay(c[i] - 1)) for i in range(num_variants)]
        penp0 = [inner(maculay(-c[i]), maculay(-c[i])) for i in range(num_variants)]

        penp = sum(penp0) + sum(penp1)
                    
        return 0.5 * phase_penalty * penp
    

    def dissipation(u, c):

        prefactor = interface(c)[1]

        diss = 0.0

        for i in range(num_variants):

            X_AM = inner(sigma(u, c), shape_strains_list[i]) \
                    - 0.5 * inner(voigt3stress(dot((Compliance_ortho - Compliance_cubic), Vstrain(u, c))), epsilon_elastic(u, c)) \
                    - phi - prefactor*(1.0 - 2*sum(c))

            cdot = mob * X_AM
            
            diss = diss + (X_AM*cdot)

        return diss


    # return elastic(u, c), contact_pen(u), phase_pen(c)
    return elastic(u, c) + interface(c)[0] + chemical(c), contact_pen(u), phase_pen(c), dissipation(u, c)

Part 2: boundaries definition, time_stepping loop

def sphere_boundaries(step):

    # Spatial coordinates of mesh
    x = SpatialCoordinate(mesh)

    # Definition of tolerance
    tol = DOLFIN_EPS # FEniCS tolerance - necessary for "boundary_markers" and the correct definition of the boundaries 


    bT = CompiledSubDomain('near(x[2], R, tol) && on_boundary', R = R, tol = tol)

    bP = CompiledSubDomain('sqrt(x[0]*x[0] + x[1]*x[1]) < meshSIZE \
                            && near(x[2], R, tol) \
                            && on_boundary', meshSIZE = meshSIZE, R=R, tol = tol)

    # Definition of contact zones
    bC = CompiledSubDomain('on_boundary && x[2] < compression + tol', R = R, tol=tol, compression = -step.values()) # Contact search - contact part of the boundary

    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    boundaries.set_all(0)
    

    bT.mark(boundaries, 2)      # top
    bP.mark(boundaries, 4)      # lateral_penalty (for stability)
    bC.mark(boundaries, 6)      # contact_zone

    print("Nodes at contact zone: ", boundaries.array()[boundaries.array() == 6])



    dom_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
    dom_markers.set_all(1)


    # Displacement boundary conditions
    bc1 = DirichletBC(W.sub(0).sub(2), step, boundaries, 2)
    bc2 = DirichletBC(W.sub(0).sub(0), Constant((0)), boundaries, 4)
    bc3 = DirichletBC(W.sub(0).sub(1), Constant((0)), boundaries, 4)


    bcs = [bc1,bc2,bc3]


    dx = Measure('dx', subdomain_data=dom_markers)
    ds = Measure('ds', subdomain_data=boundaries)

    return dx, ds, boundaries, bcs



time_end = 1.0
t_curr   = 0.0

t_curr = t_curr + time_step
d_step = compression * t_curr


# previous solution
w0 = Function(W)
w0.vector()[:] = 0.0

bdry_disp = Constant((-d_step))
dx, ds, boundaries, bcs = sphere_boundaries(bdry_disp)

# Total potential energy at the current step (linear elasticity, gradient-interface)
PiT = psi(w)[0]*dx + psi(w)[1]*ds(6) + psi(w)[2]*dx + psi(w)[3]*dx

# Total potential energy of previous time-step (from known solution)
Pi0 = psi(w0)[0]*dx

# Incremental energy (increment potential energy + increment in dissiaption)
Pi = PiT - Pi0

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, w, wx)

# Compute Jacobian of F for iterative "newton_solver"
J = derivative(F, w, dw)


for bc in bcs:
    bc.apply(w0.vector())


snes_solver_parameters = {"nonlinear_solver": "snes",              
                                "snes_solver": {"linear_solver": 'lu', #"lu",
                                                 "preconditioner": "pcmg", 
                                                 "line_search": "bt",
                                                 "maximum_iterations": 50,
                                                 "absolute_tolerance": 1E-3,
                                                 "relative_tolerance": 1E-3,
                                                 "report": True,
                                                 #"method":'vinewtonrsls',
                                                 "error_on_nonconvergence": False,
                                                 }}


while t_curr <= time_end:
    
    print("Current displacement at boundary = ", d_step, " nm and time_step = ", time_step)
    
    w.vector()[:] = w0.vector()[:]
    problem = NonlinearVariationalProblem(F, w, bcs, J=J)
    solver  = NonlinearVariationalSolver(problem)

    solver.parameters.update(snes_solver_parameters)

    iter, converged = solver.solve()

    if converged:
        t_curr = t_curr + time_step
        w0.assign(w)

    else:
        t_curr = t_curr - time_step
        time_step = time_step*(0.5)
        t_curr = t_curr + time_step

    d_step = compression * t_curr
    bdry_disp.assign(-d_step)        



_c = w.split(deepcopy=True)[1:]
mart = project(sum(_c), W2)

from vedo.dolfin import plot
plot(mart)

Side Question: It makes total sense when you asked to actively avoid redefining boundaries in the loop when form is compiled outside, how would the code give almost correct results when form is compiled inside the loop where boundaries are always redefined in every iteration and only then the form is calculated ??

I appreciate your time and kind help.

Sincerely,
Deepak

[dokken]

@Deepak_Dhariwal

I have made changes per your suggestion NOT to redefine boundaries. The result is still the same:

That is not suprising as you would like these values to update.
What I meant is that you should not redefine dx, ds and boundaries. You should instead change the values of the already existing object boundaries and dom_marker.
Note that you would have to update the BCs, as mentioned in

Changing dirichlet BC region

In most cases, you would have to redefine the bu_u object, as it applies caching when applying the first boundary condition, see: Bitbucket This can also be found in the documentation: https://fenicsproject.org/olddocs/dolfin/latest/cpp/d1/d95/classdolfin_1_1DirichletBC.html#details

I am going to give you a simple example, where the markers are updated, without the form and measure being re-defined

from dolfin import *

mesh = UnitSquareMesh(32, 32)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
ds_c = Measure("ds", domain=mesh, subdomain_data=boundaries)
L = 0.8
marker = CompiledSubDomain("on_boundary && x[0]<L", L=L)
marker.mark(boundaries, 1)
form = 1*ds_c(1)


print(f"L={marker.get_property('L')}, val: {assemble(form)}")

for i in [0.1, 0.2, 0.4, 0.8]:
    marker.set_property("L", i)
    boundaries.set_all(0)

    marker.mark(boundaries, 1)
    print(f"L={marker.get_property('L')}, val: {assemble(form)}")

gives

fenics@dokken-XPS-9320:/root/shared$ python3 mwe.py 
L=0.8, val: 2.5625
L=0.1, val: 1.1875
L=0.2, val: 1.375
L=0.4, val: 1.75
L=0.8, val: 2.5625

[Deepak_Dhariwal]

Thanks very much @dokken. I have made changes per your suggestions and getting consistent results. I am marking this form compilation as solved.

On a different note, as I am generating results, the iterative methods (GMRES) aren’t proving helpful. Instead, using one of the direct solvers (LU), the results seem to converge expectedly, and solution values are also resonable.

Should I have to fine tune my computations to move away from direct solvers, what is an ideal approach in your view?

So far, my reference has been what’s explained in this answer

I am using the following parameter set for SNES solver.

    snes_solver_parameters = {"nonlinear_solver": "snes",              
                              "snes_solver": {"linear_solver": 'gmres', #"lu",
                                              "preconditioner": "pcmg", 
                                               "line_search": "bt",
                                               "maximum_iterations": 50,
                                               "absolute_tolerance": 1E-3,
                                               "relative_tolerance": 1E-3,
                                               "report": True,
                                               #"method":'vinewtonrsls',
                                               "error_on_nonconvergence": False}}

(tolerances are kept low for the tine being to get a approximate sense of solution)

Thanks a lot.

Deepak