Rate-dependent pull-in instability

Rate-dependent pull-in instability#

  • Code for 3D large deformation electro-elasticity with u-p formulation

  • Electro-viscoelastic pull-in instability of a 3D VHB block

  • Uses quadrature representation of internal variable: viscous deformation tensor Cv.

Units#

  • Length: mm

  • Mass: kg

  • Time: s

  • Charge: nC

  • Force: mN

  • Stress: kPa

  • Energy: microJ

  • Electric potential: kV

Software:#

  • Dolfinx v0.8.0

Import modules#

# Import FEnicSx/dolfinx
import dolfinx

# For numerical arrays
import numpy as np

# For MPI-based parallelization
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

# PETSc solvers
from petsc4py import PETSc

# specific functions from dolfinx modules
from dolfinx import fem, mesh, io, plot, log
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, Expression )
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import VTXWriter, XDMFFile


# specific functions from ufl modules
import ufl
from ufl import (TestFunctions, TrialFunction, Identity, grad, det, div, dev, inv, tr, sqrt, conditional ,\
                 gt, dx, inner, derivative, dot, ln, split, exp, eq, cos, acos, ge, le, outer)

# basix finite elements
import basix
from basix.ufl import element, mixed_element, quadrature_element

# Matplotlib for plotting
import matplotlib.pyplot as plt
plt.close('all')

# For timing the code
from datetime import datetime


# Set level of detail for log messages (integer)
# Guide:
# CRITICAL  = 50, // errors that may lead to data corruption
# ERROR     = 40, // things that HAVE gone wrong
# WARNING   = 30, // things that MAY go wrong later
# INFO      = 20, // information of general interest (includes solver info)
# PROGRESS  = 16, // what's happening (broadly)
# TRACE     = 13, // what's happening (in detail)
# DBG       = 10  // sundry
#
log.set_log_level(log.LogLevel.WARNING)

Define geometry#

# A 3-D cube
length = 1.0 # mm
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0,0.0,0.0], [length,length,length]],\
                         [1,1,1], mesh.CellType.tetrahedron)

x = ufl.SpatialCoordinate(domain)

# Identify the planar boundaries of the  box mesh
#
def xBot(x):
    return np.isclose(x[0], 0)
def xTop(x):
    return np.isclose(x[0], length)
def yBot(x):
    return np.isclose(x[1], 0)
def yTop(x):
    return np.isclose(x[1], length)
def zBot(x):
    return np.isclose(x[2], 0)
def zTop(x):
    return np.isclose(x[2], length)
    
# Mark the sub-domains
boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]

# build collections of facets on each subdomain and mark them appropriately.
facet_indices, facet_markers = [], [] # initalize empty collections of indices and markers.
fdim = domain.topology.dim - 1 # geometric dimension of the facet (mesh dimension - 1)
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(domain, fdim, locator) # an array of all the facets in a 
                                                         # given subdomain ("locator")
    facet_indices.append(facets)                         # add these facets to the collection.
    facet_markers.append(np.full_like(facets, marker))   # mark them with the appropriate index.

# Format the facet indices and markers as required for use in dolfinx.
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
# 
# Add these marked facets as "mesh tags" for later use in BCs.
facet_tags = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

Print out the unique facet index numbers

top_imap = domain.topology.index_map(2)      # index map of 2D entities in domain (facets)
values = np.zeros(top_imap.size_global)      # an array of zeros of the same size as number of 2D entities
values[facet_tags.indices]=facet_tags.values # populating the array with facet tag index numbers
print(np.unique(facet_tags.values))          # printing the unique indices

[1 2 3 4 5 6]

Visualize reference configuration and boundary facets

import pyvista
pyvista.set_jupyter_backend('html')
from dolfinx.plot import vtk_mesh
pyvista.start_xvfb()

# initialize a plotter
plotter = pyvista.Plotter()

# Add the mesh.
topology, cell_types, geometry = plot.vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter.add_mesh(grid, show_edges=True) 

labels = dict(zlabel='Z', xlabel='X', ylabel='Y')
plotter.add_axes(**labels)

plotter.screenshot("results/cube_mesh.png")

from IPython.display import Image
Image(filename='results/cube_mesh.png')
../_images/6daff8b51b4539eb0f3252485dd47d1829c276757cfe442965c29fc18c01f568.png

Un-comment this cell to see an interactive visualization of the mesh#

# plotter.show()

Define boundary and volume integration measure#

# Define the boundary integration measure "ds" using the facet tags,
# also specify the number of surface quadrature points.
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tags, metadata={'quadrature_degree':4})

# Define the volume integration measure "dx" 
# also specify the number of volume quadrature points.
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 4})

# Create facet to cell connectivity required to determine boundary facets.
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim)
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

#  Define facet normal
n = ufl.FacetNormal(domain)

Material parameters#

# Equilibrium material parameters
#
rho = Constant(domain, PETSc.ScalarType(1e-6))           # 1000 kg/m^3 = 1e-6 kg/mm^3
#
Geq_0   = Constant(domain, PETSc.ScalarType(15.36))      # Shear modulus, kPa
Kbulk   = Constant(domain, PETSc.ScalarType(1e3*Geq_0))  # Bulk modulus, kPa
lambdaL = Constant(domain, PETSc.ScalarType(5.99))        # Arruda-Boyce locking stretch

# Viscoelasticity parameters from Wang et al (2016) for VHB4910
#
Gneq_1  = Constant(domain, PETSc.ScalarType(26.06))    #  Non-equilibrium shear modulus, kPa
tau_1   = Constant(domain, PETSc.ScalarType(0.6074))   #  relaxation time, s
#
Gneq_2  = Constant(domain, PETSc.ScalarType(26.53))    #  Non-equilibrium shear modulus, kPa
tau_2   = Constant(domain, PETSc.ScalarType(6.56))     #  relaxation time, s
#
Gneq_3  = Constant(domain, PETSc.ScalarType(10.83))    #  Non-equilinrium shear modulus, kPa
tau_3   = Constant(domain, PETSc.ScalarType(61.25))    #  relaxation time, s

# Electrostatic  parameters
vareps_0 = Constant(domain, PETSc.ScalarType(8.85E-3)) #  permittivity of free space pF/mm
vareps_r = Constant(domain, PETSc.ScalarType(4.8))     #  relative permittivity, dimensionless
vareps   = vareps_r*vareps_0         #  permittivity of the material

Function spaces#

U2 = element("Lagrange", domain.basix_cell(), 2, shape=(3,))  # For displacement
P1 = element("Lagrange", domain.basix_cell(), 1)  # For pressure and electric potential
P0 = quadrature_element(domain.basix_cell(), degree=2, scheme="default")
# Note: it seems that for the current version of dolfinx, 
# only degree=2 quadrature elements actually function properly 
# in e.g. visualization interpolations and problem solution.
T0 = basix.ufl.blocked_element(P0, shape=(3, 3)) # for Cv
#
TH = mixed_element([U2, P1, P1])  # Taylor-Hood style mixed element
ME = functionspace(domain, TH)    # Total space for all DOFs
#
V1 = functionspace(domain, P1) # Scalar function space.
V2 = functionspace(domain, U2) # Vector function space
V3 = functionspace(domain, T0) # Tensor function space
#
# Define actual functions with the required DOFs
w = Function(ME)
u, p, phi = split(w) # displacement u, presssure p, and electric potential phi

# A copy of functions to store values in the previous step
w_old  = Function(ME)
u_old, p_old, phi_old = split(w_old)   

# Define test functions   
u_test, p_test, phi_test = TestFunctions(ME)    

# Define trial functions needed for automatic differentiation
dw = TrialFunction(ME)

# Define a tensor-valued function for Cv.
Cv_1_old = Function(V3) 
Cv_2_old = Function(V3) 
Cv_3_old = Function(V3) 

# Functions for storing the velocity and acceleration at prev. step
v_old = Function(V2)
a_old = Function(V2)

# Initial conditions: 

# A function for constructing the identity matrix.
#
# To use the interpolate() feature, this must be defined as a 
# function of x.
def identity(x):
    values = np.zeros((domain.geometry.dim*domain.geometry.dim,
                      x.shape[1]), dtype=np.float64)
    values[0] = 1
    values[4] = 1
    values[8] = 1
    return values

# interpolate the identity onto the tensor-valued Cv function.
Cv_1_old.interpolate(identity)  
Cv_2_old.interpolate(identity) 
Cv_3_old.interpolate(identity)

Subroutines for kinematics and constitutive equations#

#------------------------------------------------------------- 
# Utility subroutines
#-------------------------------------------------------------
 
# Subroutine for a "safer" sqrt() function which avoids a divide by zero 
# when differentiated. 
def safe_sqrt(x):
    return sqrt(x + 1.0e-16)


#------------------------------------------------------------- 
# Subroutines for kinematics
#-------------------------------------------------------------

# Deformation gradient 
def F_calc(u):
    Id = Identity(3) 
    F = Id + grad(u) 
    return F

# Subrountine for computing the effective stretch
def lambdaBar_calc(u):
    F = F_calc(u)
    J = det(F)
    Fbar = J**(-1/3)*F
    Cbar = Fbar.T*Fbar
    I1 = tr(Cbar)
    lambdaBar = safe_sqrt(I1/3.0)
    return lambdaBar

#------------------------------------------------------------- 
# Subroutines for computing the viscous flow update
#-------------------------------------------------------------

# Subroutine for computing the viscous stretch Cv at the end of the step.
def Cv_update(u, Cv_old, tau_r):
    F = F_calc(u)
    J = det(F)
    C = F.T*F
    # new update from Shawn here
    Cbar = J**(-2./3.)*C
    fac = Cv_old + (dk/tau_r)*Cbar
    Cv_new = (det(fac))**(-1./3.) * fac
    return Cv_new 
 
#------------------------------------------------------------- 
# Subroutines for calculating the electric field and displacement
#-------------------------------------------------------------

# Referential electric displacement 
def Dmat_calc(u, phi):
    F = F_calc(u)
    J = det(F)
    C = F.T*F
    e_R  = - grad(phi) # referential electric field
    Dmat = vareps * J* inv(C)*e_R
    return Dmat

#------------------------------------------------------------- 
# Subroutines for calculating the Cauchy stress
#-------------------------------------------------------------

# Subroutine for computing the zeta-function in the Arruda-Boyce model.
def zeta_calc(u):
    lambdaBar = lambdaBar_calc(u)
    # Use Pade approximation of Langevin inverse (A. Cohen, 1991)
    # This is sixth-order accurate.
    z = lambdaBar/lambdaL
    z = conditional(gt(z,0.99), 0.99, z) # Keep from blowing up
    beta = z*(3.0 - z**2.0)/(1.0 - z**2.0)
    zeta = (lambdaL/(3*lambdaBar))*beta
    return zeta

# Generalized shear modulus for Arruda-Boyce model
def Geq_AB_calc(u):
    zeta = zeta_calc(u)
    Geq_AB  = Geq_0 * zeta
    return Geq_AB

# Subroutine for calculating the  equilibrium Cauchy stress
def T_eq_calc(u,p):
    F = F_calc(u)
    J = det(F)
    Fbar = J**(-1./3.)*F
    Bbar = Fbar*Fbar.T
    Geq  = Geq_AB_calc(u)
    T_eq = (1/J)* Geq * dev(Bbar) - p * Identity(3)
    return T_eq


# Subroutine for calculating the electrotatic contribution to the Cauchy stress
def T_maxw_calc(u,phi):
    F = F_calc(u)
    e_R  = - grad(phi)    # referential electric field
    e_sp = inv(F.T)*e_R   # spatial electric field 
    # Spatial Maxwel stress
    T_maxw = vareps*(outer(e_sp,e_sp) - 1/2*(inner(e_sp,e_sp))*Identity(3))
    return T_maxw

# Subroutine for the non-equilibrium Cauchy stress.
def T_neq_calc(u, Cv, Gneq):
    F  = F_calc(u)
    J = det(F)
    C = F.T*F
    # Shawn update here
    Fbar = J**(-1./3.) * F
    Cbar = J**(-2./3.) * C
    T_neq = J**(-1.) * Gneq * (Fbar * inv(Cv) * Fbar.T - (1./3.) * inner(Cbar, inv(Cv)) * Identity(3) ) 
    return T_neq

#------------------------------------------------------------- 
# Subroutine for calculating the total Piola stress
#-------------------------------------------------------------
   
# Subroutine for the total Piola stress.
def  Piola_calc(u, p, Cv_1, Cv_2, Cv_3, Gneq_1, Gneq_2, Gneq_3):
    F  = F_calc(u)
    J = det(F)
    
    T_eq = T_eq_calc(u,p)
    
    T_maxw = T_maxw_calc(u,phi)
    
    T_neq_1 = T_neq_calc(u, Cv_1, Gneq_1)
    
    T_neq_2 = T_neq_calc(u, Cv_2, Gneq_2)
        
    T_neq_3 = T_neq_calc(u, Cv_3, Gneq_3)
    
    T = T_eq + T_maxw + T_neq_1 + T_neq_2 + T_neq_3
    
    Piola = J*T*inv(F.T)
    
    return Piola

#---------------------------------------------------------------------
# Subroutine for updating  acceleration using the Newmark beta method:
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
#---------------------------------------------------------------------
def update_a(u, u_old, v_old, a_old):
    return (u-u_old-dk*v_old)/beta/dk**2 - (1-2*beta)/2/beta*a_old

#---------------------------------------------------------------------
# Subroutine for updating  velocity using the Newmark beta method
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
#---------------------------------------------------------------------
def update_v(a, u_old, v_old, a_old):
    return v_old + dk*((1-gamma)*a_old + gamma*a)

#---------------------------------------------------------------------
# alpha-method averaging function
#---------------------------------------------------------------------
def avg(x_old, x_new, alpha):
    return alpha*x_old + (1-alpha)*x_new

Evaluate kinematics and constitutive relations#

# Get acceleration and velocity at end of step
a_new = update_a(u, u_old, v_old, a_old)
v_new = update_v(a_new, u_old, v_old, a_old)

# get avg (u,p) fields for generalized-alpha method
u_avg  = avg(u_old, u, alpha)
p_avg  = avg(p_old, p, alpha)

# Kinematical quantities
F  = F_calc(u_avg)
J  = det(F)
lambdaBar = lambdaBar_calc(u_avg)

# update the Cv tensors
Cv_1 = Cv_update(u_avg, Cv_1_old, tau_1)
Cv_2 = Cv_update(u_avg, Cv_2_old, tau_2)
Cv_3 = Cv_update(u_avg, Cv_3_old, tau_3)

# Referential electric displacement
Dmat = Dmat_calc(u_avg, phi)

#  Evaulate the total Piola stress
Piola = Piola_calc(u_avg, p_avg, Cv_1, Cv_2, Cv_3, Gneq_1, Gneq_2, Gneq_3)

Weak forms#

# The weak form for the equilibrium equation
#
Res_1  =  inner( Piola, grad(u_test))*dx \
    + inner(rho*a_new, u_test)*dx 
              
# The auxiliary equation for the pressure
#
Res_2 = inner((p_avg/Kbulk + ln(J)/J) , p_test)*dx

#  The weak form for Gauss's equation
Res_3 = inner(Dmat, grad(phi_test))*dx 

# The total residual
Res = Res_1 + Res_2 + Res_3

# Automatic differentiation tangent:
a = derivative(Res, w, dw)

Set-up output files#

# results file name
results_name = "3D_pullin"

# Function space for projection of results
P1 = element("Lagrange", domain.basix_cell(), 1)
VV1 = fem.functionspace(domain, P1) # linear scalar function space
#
U1 = element("Lagrange", domain.basix_cell(), 1, shape=(3,)) 
VV2 = fem.functionspace(domain, U1) # linear Vector function space
#
T1 = element("Lagrange", domain.basix_cell(), 1, shape=(3,3)) 
VV3 = fem.functionspace(domain, T1) # linear tensor function space

# For visualization purposes, we need to re-project the stress tensor onto a linear function space before 
# we write it (and its components and the von Mises stress, etc) to the VTX file. 
#
# This is because the stress is a complicated "mixed" function of the (quadratic Lagrangian) displacements
# and the (quadrature representation) plastic strain tensor and scalar equivalent plastic strain. 
#
# First, define a function for setting up this kind of projection problem for visualization purposes:
def setup_projection(u, V):

    trial = ufl.TrialFunction(V)
    test  = ufl.TestFunction(V)   

    a = ufl.inner(trial, test)*dx
    L = ufl.inner(u, test)*dx

    projection_problem = dolfinx.fem.petsc.LinearProblem(a, L, [], \
        petsc_options={"ksp_type": "cg", "ksp_rtol": 1e-16, "ksp_atol": 1e-16, "ksp_max_it": 1000})
    
    return projection_problem

# Create a linear problem for projecting the stress tensor onto the linear tensor function space VV3.
#
tensor_projection_problem = setup_projection(Piola, VV3)
Piola_temp = tensor_projection_problem.solve()

# primary fields to write to output file
u_vis      = Function(VV2, name="disp")
p_vis      = Function(VV1, name="p")
phi_vis    = Function(VV1, name="phi")

# Mises stress
T     = Piola_temp*F.T/J
T0    = T - (1/3)*tr(T)*Identity(3)
Mises = sqrt((3/2)*inner(T0, T0))
Mises_vis= Function(VV1,name="Mises")
Mises_expr = Expression(Mises,VV1.element.interpolation_points())

# Cauchy stress components
T11 = Function(VV1)
T11.name = "T11"
T11_expr = Expression(T[0,0],VV1.element.interpolation_points())

T22 = Function(VV1)
T22.name = "T22"
T22_expr = Expression(T[1,1],VV1.element.interpolation_points())

T33 = Function(VV1)
T33.name = "T33"
T33_expr = Expression(T[2,2],VV1.element.interpolation_points())

# Stretch measure
lambdaBar_vis      = Function(VV1)
lambdaBar_vis.name = "lambdaBar"
lambdaBar_expr     = Expression(lambdaBar, VV1.element.interpolation_points())

# Volumetric deformation
J_vis      = Function(VV1)
J_vis.name = "J"
J_expr     = Expression(J, VV1.element.interpolation_points())

# set up the output VTX files.
file_results = VTXWriter(
    MPI.COMM_WORLD,
    "results/" + results_name + ".bp",
    [  # put the functions here you wish to write to output
        u_vis, p_vis, phi_vis, # DOF outputs
        Mises_vis, T11, T22, T33, # stress outputs
        lambdaBar_vis, J_vis, # Kinematical outputs
    ],
    engine="BP4",
)

def writeResults(t):
    
    # Update the output fields before writing to VTX.
    #
    u_vis.interpolate(w.sub(0))
    p_vis.interpolate(w.sub(1))
    phi_vis.interpolate(w.sub(2))
    #
    # re-project to smooth visualization of quadrature functions
    # before interpolating.
    Piola_temp = tensor_projection_problem.solve()
    Mises_vis.interpolate(Mises_expr)
    T11.interpolate(T11_expr)
    T22.interpolate(T22_expr)
    T33.interpolate(T33_expr)
    #
    lambdaBar_vis.interpolate(lambdaBar_expr)
    J_vis.interpolate(J_expr)
       
    # Finally, write output fields to VTX.
    #
    file_results.write(t)

Infrastructure for pulling out time history data (force, displacement, etc.)#

# # computing the reaction force using the stress field
# traction = dot(Piola_temp, n)
# Force    = dot(traction, n)*ds(4)
# rxnForce = fem.form(Force) 

# # infrastructure for evaluating functions at a certain point efficiently
pointForEval = np.array([length, length, length])

bb_tree = dolfinx.geometry.bb_tree(domain,domain.topology.dim)
cell_candidates = dolfinx.geometry.compute_collisions_points(bb_tree, pointForEval)
colliding_cells = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates, pointForEval).array

Boundary condtions#

# Constant for applied electric potential
phi_cons = Constant(domain,PETSc.ScalarType(phiRamp(0)))

# Find the specific DOFs which will be constrained.
xBot_u1_dofs = fem.locate_dofs_topological(ME.sub(0).sub(0), facet_tags.dim, facet_tags.find(1))
yBot_u2_dofs = fem.locate_dofs_topological(ME.sub(0).sub(1), facet_tags.dim, facet_tags.find(3))
zBot_u3_dofs = fem.locate_dofs_topological(ME.sub(0).sub(2), facet_tags.dim, facet_tags.find(5))
yTop_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(4))
yBot_phi_dofs = fem.locate_dofs_topological(ME.sub(2), facet_tags.dim, facet_tags.find(3))

# boundaries = [(1, xBot),(2,xTop),(3,yBot),(4,yTop),(5,zBot),(6,zTop)]

# building Dirichlet BCs
bcs_1 = dirichletbc(0.0, xBot_u1_dofs, ME.sub(0).sub(0))  # u1 fix - xBot
bcs_2 = dirichletbc(0.0, yBot_u2_dofs, ME.sub(0).sub(1))  # u2 fix - yBot
bcs_3 = dirichletbc(0.0, zBot_u3_dofs, ME.sub(0).sub(2))  # u3 fix - zBot
#
bcs_4 = dirichletbc(phi_cons, yTop_phi_dofs, ME.sub(2))  # phi ramp - yTop
bcs_5 = dirichletbc(0.0, yBot_phi_dofs, ME.sub(2))       # phi ground - yBot

bcs = [bcs_1, bcs_2, bcs_3, bcs_4, bcs_5]

Define the nonlinear variational problem#

# Set up nonlinear problem
problem = NonlinearProblem(Res, w, bcs, a)

# the global newton solver and params
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.atol = 1e-8
solver.max_it = 50
solver.report = True

#  The Krylov solver parameters.
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly" # "preonly" works equally well
opts[f"{option_prefix}pc_type"] = "lu" # do not use 'gamg' pre-conditioner
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"] = 30
ksp.setFromOptions()

Start calculation loop#

# Give the step a descriptive name
step = "Actuate"

# Variables for storing time history
totSteps = 1000000
timeHist0 = np.zeros(shape=[totSteps])
timeHist1 = np.zeros(shape=[totSteps]) 

# Iinitialize a counter for reporting data
ii=0

#  Set up temporary "helper" functions and expressions 
#  for updating the internal variables.
#
# For the Cv tensors:
Cv_1_temp = Function(V3)
Cv_1_expr = Expression(Cv_1, V3.element.interpolation_points())
#
Cv_2_temp = Function(V3)
Cv_2_expr = Expression(Cv_2, V3.element.interpolation_points())
#
Cv_3_temp = Function(V3)
Cv_3_expr = Expression(Cv_3, V3.element.interpolation_points())
#
# and also for the velocity and acceleration.
v_temp = Function(V2)
a_temp = Function(V2)
#
v_expr = Expression(v_new,V2.element.interpolation_points())
a_expr = Expression(a_new,V2.element.interpolation_points())

# Write initial state to file
writeResults(t=0.0)    

# print a message for simulation startup
print("------------------------------------")
print("Simulation Start")
print("------------------------------------")
# Store start time 
startTime = datetime.now()

# Time-stepping solution procedure loop
while (round(t + dt, 9) <= Ttot):
  
    # increment time
    t += dt 
    
    # increment counter
    ii += 1
    
    # update time variables in time-dependent BCs 
    phi_cons.value = phiRamp(t)
    
    # Solve the problem
    try:
      (iter, converged) = solver.solve(w)
    except: # Break the loop if solver fails
      break
        
    # Collect results from MPI ghost processes
    w.x.scatter_forward()

    # Print progress of calculation         
    if ii%1 == 0:
      now = datetime.now()
      current_time = now.strftime("%H:%M:%S")
      print("Step: {} | Increment: {} | Iterations: {}".format(step, ii, iter))
      print("dt: {} | Simulation Time: {} s | Percent of total time: {}%".format(round(dt,4), round(t,4), round(100*t/Ttot,4)))
      print() 
      
    # Write output to file
    writeResults(t)
       
    # Store time history variables at this time
    timeHist0[ii] = w.sub(0).sub(1).eval([length, length, length],colliding_cells[0])[0] # time history of displacement at a point
    timeHist1[ii] = w.sub(2).eval([length, length, length],colliding_cells[0])[0] # time history of voltage phi at a point
    
    # update internal variables 
    #
    # interpolate the values of the internal variables into their "temp" functions
    Cv_1_temp.interpolate(Cv_1_expr)
    Cv_2_temp.interpolate(Cv_2_expr)
    Cv_3_temp.interpolate(Cv_3_expr)
    #
    v_temp.interpolate(v_expr)
    a_temp.interpolate(a_expr)
    
    # Update DOFs for next step
    w_old.x.array[:] = w.x.array
    
    # update the old values of internal variables for next step
    Cv_1_old.x.array[:] = Cv_1_temp.x.array[:]
    Cv_2_old.x.array[:] = Cv_2_temp.x.array[:]
    Cv_3_old.x.array[:] = Cv_3_temp.x.array[:]
    #
    v_old.x.array[:] = v_temp.x.array[:]
    a_old.x.array[:] = a_temp.x.array[:]
    

    
# close the output file.
file_results.close()
         
# End analysis
print("-----------------------------------------")
print("End computation")                 
# Report elapsed real time for the analysis
endTime = datetime.now()
elapseTime = endTime - startTime
print("------------------------------------------")
print("Elapsed real time:  {}".format(elapseTime))
print("------------------------------------------")

------------------------------------
Simulation Start
------------------------------------
Step: Actuate | Increment: 1 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0004 s | Percent of total time: 0.4%

Step: Actuate | Increment: 2 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0008 s | Percent of total time: 0.8%

Step: Actuate | Increment: 3 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0012 s | Percent of total time: 1.2%

Step: Actuate | Increment: 4 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0016 s | Percent of total time: 1.6%

Step: Actuate | Increment: 5 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.002 s | Percent of total time: 2.0%

Step: Actuate | Increment: 6 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0024 s | Percent of total time: 2.4%

Step: Actuate | Increment: 7 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0028 s | Percent of total time: 2.8%

Step: Actuate | Increment: 8 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0032 s | Percent of total time: 3.2%

Step: Actuate | Increment: 9 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0036 s | Percent of total time: 3.6%

Step: Actuate | Increment: 10 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.004 s | Percent of total time: 4.0%

Step: Actuate | Increment: 11 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0044 s | Percent of total time: 4.4%

Step: Actuate | Increment: 12 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0048 s | Percent of total time: 4.8%

Step: Actuate | Increment: 13 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0052 s | Percent of total time: 5.2%

Step: Actuate | Increment: 14 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0056 s | Percent of total time: 5.6%

Step: Actuate | Increment: 15 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.006 s | Percent of total time: 6.0%

Step: Actuate | Increment: 16 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0064 s | Percent of total time: 6.4%

Step: Actuate | Increment: 17 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0068 s | Percent of total time: 6.8%

Step: Actuate | Increment: 18 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0072 s | Percent of total time: 7.2%

Step: Actuate | Increment: 19 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0076 s | Percent of total time: 7.6%

Step: Actuate | Increment: 20 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.008 s | Percent of total time: 8.0%

Step: Actuate | Increment: 21 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0084 s | Percent of total time: 8.4%

Step: Actuate | Increment: 22 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0088 s | Percent of total time: 8.8%

Step: Actuate | Increment: 23 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0092 s | Percent of total time: 9.2%

Step: Actuate | Increment: 24 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0096 s | Percent of total time: 9.6%

Step: Actuate | Increment: 25 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.01 s | Percent of total time: 10.0%

Step: Actuate | Increment: 26 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0104 s | Percent of total time: 10.4%

Step: Actuate | Increment: 27 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0108 s | Percent of total time: 10.8%

Step: Actuate | Increment: 28 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0112 s | Percent of total time: 11.2%

Step: Actuate | Increment: 29 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0116 s | Percent of total time: 11.6%

Step: Actuate | Increment: 30 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.012 s | Percent of total time: 12.0%

Step: Actuate | Increment: 31 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0124 s | Percent of total time: 12.4%

Step: Actuate | Increment: 32 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0128 s | Percent of total time: 12.8%

Step: Actuate | Increment: 33 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0132 s | Percent of total time: 13.2%

Step: Actuate | Increment: 34 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0136 s | Percent of total time: 13.6%

Step: Actuate | Increment: 35 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.014 s | Percent of total time: 14.0%

Step: Actuate | Increment: 36 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0144 s | Percent of total time: 14.4%

Step: Actuate | Increment: 37 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0148 s | Percent of total time: 14.8%

Step: Actuate | Increment: 38 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0152 s | Percent of total time: 15.2%

Step: Actuate | Increment: 39 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0156 s | Percent of total time: 15.6%

Step: Actuate | Increment: 40 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.016 s | Percent of total time: 16.0%

Step: Actuate | Increment: 41 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0164 s | Percent of total time: 16.4%

Step: Actuate | Increment: 42 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0168 s | Percent of total time: 16.8%

Step: Actuate | Increment: 43 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0172 s | Percent of total time: 17.2%

Step: Actuate | Increment: 44 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0176 s | Percent of total time: 17.6%

Step: Actuate | Increment: 45 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.018 s | Percent of total time: 18.0%

Step: Actuate | Increment: 46 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0184 s | Percent of total time: 18.4%

Step: Actuate | Increment: 47 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0188 s | Percent of total time: 18.8%

Step: Actuate | Increment: 48 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0192 s | Percent of total time: 19.2%

Step: Actuate | Increment: 49 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0196 s | Percent of total time: 19.6%

Step: Actuate | Increment: 50 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.02 s | Percent of total time: 20.0%

Step: Actuate | Increment: 51 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0204 s | Percent of total time: 20.4%

Step: Actuate | Increment: 52 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0208 s | Percent of total time: 20.8%

Step: Actuate | Increment: 53 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0212 s | Percent of total time: 21.2%

Step: Actuate | Increment: 54 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0216 s | Percent of total time: 21.6%

Step: Actuate | Increment: 55 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.022 s | Percent of total time: 22.0%

Step: Actuate | Increment: 56 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0224 s | Percent of total time: 22.4%

Step: Actuate | Increment: 57 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0228 s | Percent of total time: 22.8%

Step: Actuate | Increment: 58 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0232 s | Percent of total time: 23.2%

Step: Actuate | Increment: 59 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0236 s | Percent of total time: 23.6%

Step: Actuate | Increment: 60 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.024 s | Percent of total time: 24.0%

Step: Actuate | Increment: 61 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0244 s | Percent of total time: 24.4%

Step: Actuate | Increment: 62 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0248 s | Percent of total time: 24.8%

Step: Actuate | Increment: 63 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0252 s | Percent of total time: 25.2%

Step: Actuate | Increment: 64 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0256 s | Percent of total time: 25.6%

Step: Actuate | Increment: 65 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.026 s | Percent of total time: 26.0%

Step: Actuate | Increment: 66 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0264 s | Percent of total time: 26.4%

Step: Actuate | Increment: 67 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0268 s | Percent of total time: 26.8%

Step: Actuate | Increment: 68 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0272 s | Percent of total time: 27.2%

Step: Actuate | Increment: 69 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0276 s | Percent of total time: 27.6%

Step: Actuate | Increment: 70 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.028 s | Percent of total time: 28.0%

Step: Actuate | Increment: 71 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0284 s | Percent of total time: 28.4%

Step: Actuate | Increment: 72 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0288 s | Percent of total time: 28.8%

Step: Actuate | Increment: 73 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0292 s | Percent of total time: 29.2%

Step: Actuate | Increment: 74 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0296 s | Percent of total time: 29.6%

Step: Actuate | Increment: 75 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.03 s | Percent of total time: 30.0%

Step: Actuate | Increment: 76 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0304 s | Percent of total time: 30.4%

Step: Actuate | Increment: 77 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0308 s | Percent of total time: 30.8%

Step: Actuate | Increment: 78 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0312 s | Percent of total time: 31.2%

Step: Actuate | Increment: 79 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0316 s | Percent of total time: 31.6%

Step: Actuate | Increment: 80 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.032 s | Percent of total time: 32.0%

Step: Actuate | Increment: 81 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0324 s | Percent of total time: 32.4%

Step: Actuate | Increment: 82 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0328 s | Percent of total time: 32.8%

Step: Actuate | Increment: 83 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0332 s | Percent of total time: 33.2%

Step: Actuate | Increment: 84 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0336 s | Percent of total time: 33.6%

Step: Actuate | Increment: 85 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.034 s | Percent of total time: 34.0%

Step: Actuate | Increment: 86 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0344 s | Percent of total time: 34.4%

Step: Actuate | Increment: 87 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0348 s | Percent of total time: 34.8%

Step: Actuate | Increment: 88 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0352 s | Percent of total time: 35.2%

Step: Actuate | Increment: 89 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0356 s | Percent of total time: 35.6%

Step: Actuate | Increment: 90 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.036 s | Percent of total time: 36.0%

Step: Actuate | Increment: 91 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0364 s | Percent of total time: 36.4%

Step: Actuate | Increment: 92 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0368 s | Percent of total time: 36.8%

Step: Actuate | Increment: 93 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0372 s | Percent of total time: 37.2%

Step: Actuate | Increment: 94 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0376 s | Percent of total time: 37.6%

Step: Actuate | Increment: 95 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.038 s | Percent of total time: 38.0%

Step: Actuate | Increment: 96 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0384 s | Percent of total time: 38.4%

Step: Actuate | Increment: 97 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0388 s | Percent of total time: 38.8%

Step: Actuate | Increment: 98 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0392 s | Percent of total time: 39.2%

Step: Actuate | Increment: 99 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0396 s | Percent of total time: 39.6%

Step: Actuate | Increment: 100 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.04 s | Percent of total time: 40.0%

Step: Actuate | Increment: 101 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0404 s | Percent of total time: 40.4%

Step: Actuate | Increment: 102 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0408 s | Percent of total time: 40.8%

Step: Actuate | Increment: 103 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0412 s | Percent of total time: 41.2%

Step: Actuate | Increment: 104 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0416 s | Percent of total time: 41.6%

Step: Actuate | Increment: 105 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.042 s | Percent of total time: 42.0%

Step: Actuate | Increment: 106 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0424 s | Percent of total time: 42.4%

Step: Actuate | Increment: 107 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0428 s | Percent of total time: 42.8%

Step: Actuate | Increment: 108 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0432 s | Percent of total time: 43.2%

Step: Actuate | Increment: 109 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0436 s | Percent of total time: 43.6%

Step: Actuate | Increment: 110 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.044 s | Percent of total time: 44.0%

Step: Actuate | Increment: 111 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0444 s | Percent of total time: 44.4%

Step: Actuate | Increment: 112 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0448 s | Percent of total time: 44.8%

Step: Actuate | Increment: 113 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0452 s | Percent of total time: 45.2%

Step: Actuate | Increment: 114 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0456 s | Percent of total time: 45.6%

Step: Actuate | Increment: 115 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.046 s | Percent of total time: 46.0%

Step: Actuate | Increment: 116 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0464 s | Percent of total time: 46.4%

Step: Actuate | Increment: 117 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0468 s | Percent of total time: 46.8%

Step: Actuate | Increment: 118 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0472 s | Percent of total time: 47.2%

Step: Actuate | Increment: 119 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0476 s | Percent of total time: 47.6%

Step: Actuate | Increment: 120 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.048 s | Percent of total time: 48.0%

Step: Actuate | Increment: 121 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0484 s | Percent of total time: 48.4%

Step: Actuate | Increment: 122 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0488 s | Percent of total time: 48.8%

Step: Actuate | Increment: 123 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0492 s | Percent of total time: 49.2%

Step: Actuate | Increment: 124 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0496 s | Percent of total time: 49.6%

Step: Actuate | Increment: 125 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.05 s | Percent of total time: 50.0%

Step: Actuate | Increment: 126 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0504 s | Percent of total time: 50.4%

Step: Actuate | Increment: 127 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0508 s | Percent of total time: 50.8%

Step: Actuate | Increment: 128 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0512 s | Percent of total time: 51.2%

Step: Actuate | Increment: 129 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0516 s | Percent of total time: 51.6%

Step: Actuate | Increment: 130 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.052 s | Percent of total time: 52.0%

Step: Actuate | Increment: 131 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0524 s | Percent of total time: 52.4%

Step: Actuate | Increment: 132 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0528 s | Percent of total time: 52.8%

Step: Actuate | Increment: 133 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0532 s | Percent of total time: 53.2%

Step: Actuate | Increment: 134 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0536 s | Percent of total time: 53.6%

Step: Actuate | Increment: 135 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.054 s | Percent of total time: 54.0%

Step: Actuate | Increment: 136 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0544 s | Percent of total time: 54.4%

Step: Actuate | Increment: 137 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0548 s | Percent of total time: 54.8%

Step: Actuate | Increment: 138 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0552 s | Percent of total time: 55.2%

Step: Actuate | Increment: 139 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0556 s | Percent of total time: 55.6%

Step: Actuate | Increment: 140 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.056 s | Percent of total time: 56.0%

Step: Actuate | Increment: 141 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0564 s | Percent of total time: 56.4%

Step: Actuate | Increment: 142 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0568 s | Percent of total time: 56.8%

Step: Actuate | Increment: 143 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0572 s | Percent of total time: 57.2%

Step: Actuate | Increment: 144 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0576 s | Percent of total time: 57.6%

Step: Actuate | Increment: 145 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.058 s | Percent of total time: 58.0%

Step: Actuate | Increment: 146 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0584 s | Percent of total time: 58.4%

Step: Actuate | Increment: 147 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0588 s | Percent of total time: 58.8%

Step: Actuate | Increment: 148 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0592 s | Percent of total time: 59.2%

Step: Actuate | Increment: 149 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0596 s | Percent of total time: 59.6%

Step: Actuate | Increment: 150 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.06 s | Percent of total time: 60.0%

Step: Actuate | Increment: 151 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0604 s | Percent of total time: 60.4%

Step: Actuate | Increment: 152 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0608 s | Percent of total time: 60.8%

Step: Actuate | Increment: 153 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0612 s | Percent of total time: 61.2%

Step: Actuate | Increment: 154 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0616 s | Percent of total time: 61.6%

Step: Actuate | Increment: 155 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.062 s | Percent of total time: 62.0%

Step: Actuate | Increment: 156 | Iterations: 3
dt: 0.0004 | Simulation Time: 0.0624 s | Percent of total time: 62.4%

Step: Actuate | Increment: 157 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0628 s | Percent of total time: 62.8%

Step: Actuate | Increment: 158 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0632 s | Percent of total time: 63.2%

Step: Actuate | Increment: 159 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0636 s | Percent of total time: 63.6%

Step: Actuate | Increment: 160 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.064 s | Percent of total time: 64.0%

Step: Actuate | Increment: 161 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0644 s | Percent of total time: 64.4%

Step: Actuate | Increment: 162 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0648 s | Percent of total time: 64.8%

Step: Actuate | Increment: 163 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0652 s | Percent of total time: 65.2%

Step: Actuate | Increment: 164 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0656 s | Percent of total time: 65.6%

Step: Actuate | Increment: 165 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.066 s | Percent of total time: 66.0%

Step: Actuate | Increment: 166 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0664 s | Percent of total time: 66.4%

Step: Actuate | Increment: 167 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0668 s | Percent of total time: 66.8%

Step: Actuate | Increment: 168 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0672 s | Percent of total time: 67.2%

Step: Actuate | Increment: 169 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0676 s | Percent of total time: 67.6%

Step: Actuate | Increment: 170 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.068 s | Percent of total time: 68.0%

Step: Actuate | Increment: 171 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0684 s | Percent of total time: 68.4%

Step: Actuate | Increment: 172 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0688 s | Percent of total time: 68.8%

Step: Actuate | Increment: 173 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0692 s | Percent of total time: 69.2%

Step: Actuate | Increment: 174 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0696 s | Percent of total time: 69.6%

Step: Actuate | Increment: 175 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.07 s | Percent of total time: 70.0%

Step: Actuate | Increment: 176 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0704 s | Percent of total time: 70.4%

Step: Actuate | Increment: 177 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0708 s | Percent of total time: 70.8%

Step: Actuate | Increment: 178 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0712 s | Percent of total time: 71.2%

Step: Actuate | Increment: 179 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0716 s | Percent of total time: 71.6%

Step: Actuate | Increment: 180 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.072 s | Percent of total time: 72.0%

Step: Actuate | Increment: 181 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0724 s | Percent of total time: 72.4%

Step: Actuate | Increment: 182 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0728 s | Percent of total time: 72.8%

Step: Actuate | Increment: 183 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0732 s | Percent of total time: 73.2%

Step: Actuate | Increment: 184 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0736 s | Percent of total time: 73.6%

Step: Actuate | Increment: 185 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.074 s | Percent of total time: 74.0%

Step: Actuate | Increment: 186 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0744 s | Percent of total time: 74.4%

Step: Actuate | Increment: 187 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0748 s | Percent of total time: 74.8%

Step: Actuate | Increment: 188 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0752 s | Percent of total time: 75.2%

Step: Actuate | Increment: 189 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0756 s | Percent of total time: 75.6%

Step: Actuate | Increment: 190 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.076 s | Percent of total time: 76.0%

Step: Actuate | Increment: 191 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0764 s | Percent of total time: 76.4%

Step: Actuate | Increment: 192 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0768 s | Percent of total time: 76.8%

Step: Actuate | Increment: 193 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0772 s | Percent of total time: 77.2%

Step: Actuate | Increment: 194 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0776 s | Percent of total time: 77.6%

Step: Actuate | Increment: 195 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.078 s | Percent of total time: 78.0%

Step: Actuate | Increment: 196 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0784 s | Percent of total time: 78.4%

Step: Actuate | Increment: 197 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0788 s | Percent of total time: 78.8%

Step: Actuate | Increment: 198 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0792 s | Percent of total time: 79.2%

Step: Actuate | Increment: 199 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0796 s | Percent of total time: 79.6%

Step: Actuate | Increment: 200 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.08 s | Percent of total time: 80.0%

Step: Actuate | Increment: 201 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0804 s | Percent of total time: 80.4%

Step: Actuate | Increment: 202 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0808 s | Percent of total time: 80.8%

Step: Actuate | Increment: 203 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0812 s | Percent of total time: 81.2%

Step: Actuate | Increment: 204 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0816 s | Percent of total time: 81.6%

Step: Actuate | Increment: 205 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.082 s | Percent of total time: 82.0%

Step: Actuate | Increment: 206 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0824 s | Percent of total time: 82.4%

Step: Actuate | Increment: 207 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0828 s | Percent of total time: 82.8%

Step: Actuate | Increment: 208 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0832 s | Percent of total time: 83.2%

Step: Actuate | Increment: 209 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0836 s | Percent of total time: 83.6%

Step: Actuate | Increment: 210 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.084 s | Percent of total time: 84.0%

Step: Actuate | Increment: 211 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0844 s | Percent of total time: 84.4%

Step: Actuate | Increment: 212 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0848 s | Percent of total time: 84.8%

Step: Actuate | Increment: 213 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0852 s | Percent of total time: 85.2%

Step: Actuate | Increment: 214 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0856 s | Percent of total time: 85.6%

Step: Actuate | Increment: 215 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.086 s | Percent of total time: 86.0%

Step: Actuate | Increment: 216 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0864 s | Percent of total time: 86.4%

Step: Actuate | Increment: 217 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0868 s | Percent of total time: 86.8%

Step: Actuate | Increment: 218 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0872 s | Percent of total time: 87.2%

Step: Actuate | Increment: 219 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0876 s | Percent of total time: 87.6%

Step: Actuate | Increment: 220 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.088 s | Percent of total time: 88.0%

Step: Actuate | Increment: 221 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0884 s | Percent of total time: 88.4%

Step: Actuate | Increment: 222 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0888 s | Percent of total time: 88.8%

Step: Actuate | Increment: 223 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0892 s | Percent of total time: 89.2%

Step: Actuate | Increment: 224 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0896 s | Percent of total time: 89.6%

Step: Actuate | Increment: 225 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.09 s | Percent of total time: 90.0%

Step: Actuate | Increment: 226 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0904 s | Percent of total time: 90.4%

Step: Actuate | Increment: 227 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0908 s | Percent of total time: 90.8%

Step: Actuate | Increment: 228 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0912 s | Percent of total time: 91.2%

Step: Actuate | Increment: 229 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0916 s | Percent of total time: 91.6%

Step: Actuate | Increment: 230 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.092 s | Percent of total time: 92.0%

Step: Actuate | Increment: 231 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0924 s | Percent of total time: 92.4%

Step: Actuate | Increment: 232 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0928 s | Percent of total time: 92.8%

Step: Actuate | Increment: 233 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0932 s | Percent of total time: 93.2%

Step: Actuate | Increment: 234 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0936 s | Percent of total time: 93.6%

Step: Actuate | Increment: 235 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.094 s | Percent of total time: 94.0%

Step: Actuate | Increment: 236 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0944 s | Percent of total time: 94.4%

Step: Actuate | Increment: 237 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0948 s | Percent of total time: 94.8%

Step: Actuate | Increment: 238 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0952 s | Percent of total time: 95.2%

Step: Actuate | Increment: 239 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.0956 s | Percent of total time: 95.6%

Step: Actuate | Increment: 240 | Iterations: 4
dt: 0.0004 | Simulation Time: 0.096 s | Percent of total time: 96.0%

Step: Actuate | Increment: 241 | Iterations: 5
dt: 0.0004 | Simulation Time: 0.0964 s | Percent of total time: 96.4%

Step: Actuate | Increment: 242 | Iterations: 5
dt: 0.0004 | Simulation Time: 0.0968 s | Percent of total time: 96.8%

Step: Actuate | Increment: 243 | Iterations: 7
dt: 0.0004 | Simulation Time: 0.0972 s | Percent of total time: 97.2%

-----------------------------------------
End computation
------------------------------------------
Elapsed real time:  0:00:01.071366
------------------------------------------

Plot results#

# set plot font to size 18
font = {'size'   : 18}
plt.rc('font', **font)

# Get array of default plot colors
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

# Plot the normalized dimensionless quantity for $\phi$ used in Wang et al. 2016
# versus stretch in the vertical direction.
#
normVolts = timeHist1/(length * np.sqrt(float(Geq_0)/float(vareps))) 
normVolts = normVolts[0:ii]
#
stretch = timeHist0/length + 1.0
stretch = stretch[0:ii]
#
plt.plot(normVolts, stretch, c=colors[0], linewidth=2.0)
plt.grid(linestyle="--", linewidth=0.5, color='b')
ax = plt.gca()
#
ax.set_ylabel(r'$\lambda$')
# ax.set_ylim([0.6,1.05])
# ax.set_yticks([0.2, 0.4, 0.6, 0.8, 1.0])
#
# ax.set_xlabel(r'$(\phi/\ell_0)/\sqrt{G_0/\varepsilon} $')
ax.set_xlabel(r'Normalized potential, $\frac{\phi/\ell_0}{\sqrt{G_0/\varepsilon}}$',size=18)
# ax.set_xlim([0.0,0.8])
# ax.set_xticks([0.0, 0.25, 0.5, 0.75, 1.0])
#
# from matplotlib.ticker import AutoMinorLocator,FormatStrFormatter
# ax.xaxis.set_minor_locator(AutoMinorLocator())
# ax.yaxis.set_minor_locator(AutoMinorLocator())
# plt.show()


plt.text(normVolts[ii-1], stretch[ii-1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')



fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/cube_pullin_Im175.png", dpi=600)
../_images/d37724a93b6c991bd2cec4f76583edb8b3682315e9c0d53762c39209021c1e84.png

Save data for later, then plot all the data at once#

  • I commented out the viscous stress contributions to create the equilibrium simulation

# np.savetxt('rateEq.csv', np.vstack((normVolts,stretch)).T, delimiter=',')

# retrieve data from prior runs
rateEq = np.genfromtxt('rateEq.csv', delimiter=',')
rate1m3 = np.genfromtxt('rate1e-3.csv', delimiter=',')
rate1m2 = np.genfromtxt('rate1e-2.csv', delimiter=',')
rate1m1 = np.genfromtxt('rate1e-1.csv', delimiter=',')
rate1m0 = np.genfromtxt('rate1e-0.csv', delimiter=',')
rate1p1 = np.genfromtxt('rate1e1.csv', delimiter=',')


fig = plt.figure()
plt.plot(rateEq[:,0],rateEq[:,1], c=colors[0], linewidth=2.0, label=r'Equil.')
plt.plot(rate1m3[:,0],rate1m3[:,1], c=colors[1], linewidth=2.0, label=r'$10^{-3}$/s')
plt.plot(rate1m2[:,0],rate1m2[:,1], c=colors[2], linewidth=2.0, label=r'$10^{-2}$/s')
plt.plot(rate1m1[:,0],rate1m1[:,1], c=colors[3], linewidth=2.0, label=r'$10^{-1}$/s')
plt.plot(rate1m0[:,0],rate1m0[:,1], c=colors[4], linewidth=2.0, label=r'$10^{0}$/s')
plt.plot(rate1p1[:,0],rate1p1[:,1], c=colors[5], linewidth=2.0, label=r'$10^{1}$/s')
plt.grid(linestyle="--", linewidth=0.5, color='b')
plt.legend(loc="lower left")
ax = plt.gca()
#
ax.set_ylabel(r'$\lambda$')
ax.set_xlabel(r'Normalized potential, $\frac{\phi/\ell_0}{\sqrt{G_0/\varepsilon}}$',size=18)
ax.set_ylim([0.35,1.05])
ax.set_xlim([0.0,1.65])


plt.text(rateEq[-1,0],rateEq[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')
plt.text(rate1m3[-1,0],rate1m3[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')
plt.text(rate1m2[-1,0],rate1m2[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')
plt.text(rate1m1[-1,0],rate1m1[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')
plt.text(rate1m0[-1,0],rate1m0[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')
plt.text(rate1p1[-1,0],rate1p1[-1,1], r'X', color='k',\
           horizontalalignment='center',verticalalignment='center')


fig = plt.gcf()
fig.set_size_inches(7,5)
plt.tight_layout()
plt.savefig("results/cube_pullin_RateSweep.png", dpi=600)
../_images/46b4080c348343ebf3c5bc60c0840cee10b69052a2f742e51e7b6e9eff69a9fc.png