# Dynamic Linear elasticity problem

Hi,
I am interested in simulating the dynamic extension for a linear elastic material. Its a simple geometry which is to be pulled from the top, whose extension is a linear function of time.

I am using the equations of linear elasticity to figure out the stress and strain generated. But I am unable to get the desired results, infact the extensions are not even equal to those for boundary conditions

Though the loss function is in orders of 10^-5 even then there is visible error. Thus I am struggling with genrating required results.

Is there is problem with configuration, batch size, learning rate or with the physics?

Scripts:

1. config
``````# Copyright (c) 2023, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

defaults :
- modulus_default
- arch:
- fully_connected
- scheduler: tf_exponential_lr
- loss: sum
- _self_

# for time varied problems
jit: false

arch:
fully_connected:
layer_size: 128
nr_layers: 4

scheduler:
decay_rate: 0.95
decay_steps: 10000

training:
rec_validation_freq: 100
rec_results_freq : 100
rec_constraint_freq: 100
max_steps : 2000

batch_size:
panel_left: 4400
panel_right: 4400
panel_bottom: 1600
panel_top: 1600
lr_interior: 20000

``````
1. the running file
``````
# import the dependicies
from ast import And
from sympy import Symbol, Eq
import numpy as np

import modulus.sym
from modulus.sym.hydra import instantiate_arch, ModulusConfig
from modulus.sym.models.layers import Activation
from modulus.sym.solver import Solver
from modulus.sym.domain import Domain
from modulus.sym.geometry.primitives_2d import Rectangle, Circle
from modulus.sym.domain.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
)
from modulus.sym.domain.validator import PointwiseValidator
from modulus.sym.domain.inferencer import PointwiseInferencer
from modulus.sym.key import Key
from modulus.sym.node import Node
from modulus.sym.utils.io.vtk import VTKUniformGrid
from modulus.sym.domain.inferencer import PointVTKInferencer

from GE import Model2D
"""
Custom model specific for out problem

PDE solving classes are to scripted here
"""

# material specifics
E = 210e9 # Pa; young's modulus
nu = 0.33 # poisons ratio

# model specifics
beta = 1/8100 # N; damaage capturing term
r1 = 1e4 # Pa; Linear term for /delta^{diss}
r2 = 1e4 # Pa-s; Qudratic term for /delta^{diss}

"""
Working variables: Rate dependent gradient {in exp decay of damage}
/phi_{0} = elatic energy stored

/phi = total energy for

/delta^{diss} = damage dissapiation term {Pa/s}

f = damage function

/beta = regularization term
"""

# create run; setup the config filesconda
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))  # Pa
mu_real = E / (2 * (1 + nu))  # Pa
lambda_ = lambda_ / mu_real  # Dimensionless
mu = 1.0  # Dimensionless

u_max = 0.01 # m, say for the sake of example here

# Using LE for now
Equations = Model2D(E = E, nu = nu, lambda_=lambda_, mu = mu, rho = 1, beta = beta, flag = 1, r1 = r1, r2 = r2)

"""
Neural Network design:

INPUT: ((x,y), t) {space-time coordinates}
OUTPUT: ((u,v), d) {extensions and damage}

Number of hidden layers and number of neurons to decided after tweaking with values

LE setup:

"""
# asumming 2D: planar stress state
net = instantiate_arch(
# INPUTs
input_keys=[
Key("x"),
Key("y"),
Key("t")],
# OUTPUTs
output_keys=[
Key("u"),
Key("v"),
Key("d"),
],
cfg=cfg.arch.fully_connected, # change size from the config.yaml
activation_fn=Activation.TANH,

)
nodes = Equations.make_nodes() + [net.make_node(name="virtual_work_model")]

# Generate the sample geometry
# coordinates for working field
x, y, t = Symbol("x"), Symbol("y"), Symbol("t")
time_range = {t : (0, 1)}
panel_origin = (-.02, -.055) # (m, m) left bottom corner
panel_dim = (0.04, 0.11) # (m, m)
circle1_origin = (-.02, -.01) # (m, m)
circle2_origin = (.02, .01) # (m, m)
panel = Rectangle(panel_origin, (panel_origin[0] + panel_dim[0], panel_origin[1] + panel_dim[1]))
thckness = 3e-3 #m

# final geometry after remove those circle form the rectangular panel
geo = panel # - c1 - c2

# geometry 2 in the paper:
rec_origin = (0, 0) # (m, m)
u0 = 1e-3 # m

# bounds on the geometry
bounds_x = (panel_origin[0], panel_origin[0] + panel_dim[0])
bounds_y = (panel_origin[1], panel_origin[1] + panel_dim[1])

# bound on time scale
bounds_t = (0, 1)

# DOMAIN
domain = Domain()

# Boundary conditions and Initial Conditions
# left wall: no force / tractions, free surface
panel_left_1 = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"traction_x": 0.0, "traction_y": 0.0},
batch_size = cfg.batch_size.panel_left,
criteria = Eq(x, panel_origin[0]), # equations to be need to defined here of constarints
parameterization={t: bounds_t, x:panel_origin[0]}, # t is a sympy symbol!!
batch_per_epoch=500,
)

# right wall
panel_right_1 = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"traction_x": 0.0, "traction_y": 0.0},
batch_size = cfg.batch_size.panel_right,
criteria = Eq(x, panel_origin[0] + panel_dim[0]), # defining the equations (constraints)
parameterization={t: bounds_t, x:(panel_origin[0] + panel_dim[0])}, # parametric terms
batch_per_epoch=500,
)

# bottom wall
panel_bottom = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar = {"v": 0.0, "d": 0.0},
batch_size = cfg.batch_size.panel_bottom,
criteria = Eq(y, panel_origin[1]),
parameterization={t: bounds_t, y: panel_origin[1]}, # t is a sympy symbol!!
batch_per_epoch=500,
)

# top wall
# u_top = constant*t {the pseudo time constarint introduced here}
panel_top = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"v": u_max*t},
batch_size = cfg.batch_size.panel_top,
criteria = Eq(y, panel_origin[1] + panel_dim[1]),
parameterization={t: bounds_t, y:(panel_origin[1] + panel_dim[1])},
batch_per_epoch=500,

)

# inferenceer if required

# constarints to be hold: interior
LE = PointwiseInteriorConstraint(
nodes=nodes,
geometry=geo,
outvar={
"equilibrium_x": 0.0,
"equilibrium_y": 0.0,
"DM":0.0,
},
batch_size=cfg.batch_size.lr_interior,
bounds={x: bounds_x, y: bounds_y, t: bounds_t},
lambda_weighting={
"equilibrium_x": Symbol("sdf"),
"equilibrium_y": Symbol("sdf"),
"DM": Symbol("sdf"),
},
parameterization={t: bounds_t, x: bounds_x, y: bounds_y},
batch_per_epoch=500,
)

nx = 400
ny = 1100

for i, specific_time in enumerate(np.linspace(0, bounds_t[1], 10)):
vtk_obj = VTKUniformGrid(
# give the geometric bounds below
bounds=[bounds_x, bounds_y],
# set the required number of points
npoints=[nx, ny],
# export map for the same
export_map={"U": ["u", "v"], "D": ["d"], "F":["f"]},
)
grid_inference = PointVTKInferencer(
vtk_obj=vtk_obj,
nodes=nodes,
input_vtk_map={"x": "x", "y": "y"}, # geometry
output_names=["u", "v", "d", "f"],
# create the time varible at instant that would be constant
invar={"t": np.full([nx*ny, 1], specific_time)},
batch_size=1000,
)
domain.add_inferencer(grid_inference, name="time_slice_" + str(i).zfill(3) + "_&_actual_time = " + str(specific_time))

# initial condition @t = 0 => d = 0
IC = PointwiseInteriorConstraint(
nodes = nodes,
geometry = geo,
outvar={"d": 0.0},
batch_size=cfg.batch_size.lr_interior,
# criteria = Eq(t, bounds_t[0]),
# bounds={x: bounds_x, y: bounds_y},
lambda_weighting={"d": 1.0},
parameterization={t : 0.0},
batch_per_epoch=500,
)

# make solver
slv = Solver(cfg, domain)

# start solver
slv.solve()

if __name__ == "__main__":
run()

``````
1. The GE.py file; support script
``````
import numpy as np
from sympy import Symbol, Function, Number
from modulus.sym.eq.pde import PDE
from modulus.sym.node import Node

class Model2D(PDE):
# 2D stress case if assumed here; Linear Elasticity
def __init__(self, E=None, nu=None, lambda_=None, mu=None, rho=1, beta=0, flag=0, r1=1000, r2=0):

# spacial coordinates
x, y = Symbol("x"), Symbol("y")
normal_x, normal_y = Symbol("normal_x"), Symbol("normal_y")

# time coordinates
t = Symbol("t")

# material properties; in case not mentioned thus adding (E, /nu) to outputs of the neural network
if lambda_ is None:
if isinstance(nu, str):
nu = Function(nu)(*input_variables)
elif isinstance(nu, (float, int)):
nu = Number(nu)
if isinstance(E, str):
E = Function(E)(*input_variables)
elif isinstance(E, (float, int)):
E = Number(E)
lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
else:
if isinstance(lambda_, str):
lambda_ = Function(lambda_)(*input_variables)
elif isinstance(lambda_, (float, int)):
lambda_ = Number(lambda_)
if isinstance(mu, str):
mu = Function(mu)(*input_variables)
elif isinstance(mu, (float, int)):
mu = Number(mu)
if isinstance(rho, str):
rho = Function(rho)(*input_variables)
elif isinstance(rho, (float, int)):
rho = Number(rho)

# make input variables
input_variables = {"x": x, "y": y, "t": t}

# Symbols
u, v = Symbol("u"), Symbol("v")
sigma_xx, sigma_xy, sigma_yy = Symbol("sigma_xx"), Symbol("sigma_xy"), Symbol("sigma_yy")

# displacement componets
u = Function("u")(*input_variables)
v = Function("v")(*input_variables)

# Linear Elasticity: Energy equation
# strains
e_xx = u.diff(x)
e_yy = v.diff(y)
e_xy = 0.5*(v.diff(x) + u.diff(y))

# stress components predicted
sigma_xx = lambda_ * (e_xx + e_yy + (-lambda_ / (lambda_ + 2 * mu) * (u.diff(x) + v.diff(y)))) + 2 * mu * e_xx
sigma_yy = lambda_ * (e_xx + e_yy + (-lambda_ / (lambda_ + 2 * mu) * (u.diff(x) + v.diff(y)))) + 2 * mu * e_yy
sigma_xy = mu * e_xy
w_z = -lambda_ / (lambda_ + 2 * mu) * (u.diff(x) + v.diff(y))

# set equations
self.equations = {}

# Equations of equilibrium: COM
self.equations["equilibrium_x"] = rho * ((u.diff(t)).diff(t)) - (sigma_xx.diff(x) + sigma_xy.diff(y))
self.equations["equilibrium_y"] = rho * ((v.diff(t)).diff(t)) - (sigma_xy.diff(x) + sigma_yy.diff(y))

# Traction equations
self.equations["traction_x"] = normal_x * sigma_xx + normal_y * sigma_xy
self.equations["traction_y"] = normal_x * sigma_xy + normal_y * sigma_yy

# phi_{0}
self.equations["phi_0"] = 0.5*(sigma_xx*e_xx + sigma_yy*e_yy + 2*sigma_xy*e_xy)

# damage variable
d = Function("d")(*input_variables)
self.equations["d"] = 0.5*(d + (d**2)**(0.5)) # ensuring always positive

if flag == 0:
# exponential
# f = 1 - (np.exp(2*self.equations["d"]) - 1)/(np.exp(2*self.equations["d"]) + 1) # f = 1 - tanh(d)
f = 1 - np.exp(-1*self.equations["d"])
else:
# linear
f = 1 - d

self.equations["f"] = f # damage function

# /phi
self.equations["phi"] = f*self.equations["phi_0"] + 0.5*beta*((f.diff(x))**2 + (f.diff(y))**2)

# /delta^{diss}
self.equations["diss"] = r1*(f.diff(t)) + 0.5*r2*(f.diff(t))**2

# p
self.equations["p"] = -1*self.equations["phi_0"] + beta*((f.diff(x)).diff(x) + (f.diff(y)).diff(y))

# p_hat
self.equations["p_hat"] = self.equations["p"] - r1

# p_plus
self.equations["p_plus"] = 0.5*(self.equations["p_hat"] + (self.equations["p_hat"]**2)**0.5)

# DM
self.equations["DM"] = r2*f.diff(t) + self.equations["p_plus"]

``````
1. I am trying to capture the transient extension in solid as: Please check if this correct method to do so?
``````   # top wall
# u_top = constant*t {the pseudo time constarint introduced here}
panel_top = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"v": u_max*t},
batch_size = cfg.batch_size.panel_top,
criteria = Eq(y, panel_origin[1] + panel_dim[1]),
parameterization={t: bounds_t, y:(panel_origin[1] + panel_dim[1])},
batch_per_epoch=500,

)
``````
1. As the results are at t = t_max :

but U_y should have been max and the top region should have been red (max) as per the point 4 but its not, which forces me to question whether my model is correct or not.

Thank you,

If a particular boundary or region is not converging fully (assuming your PDE and problem set are all good), consider tuning the `lambda_weighting` of your constraints. For example, if your top region is incorrect try increasing the weighting of the top boundary constraint.

If things still look incorrect, then perhaps its a set up issue. But try experimenting with the weighting first, this is looking like you may be in the hyper-parameter search/tuning phase.