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.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

defaults :
  - modulus_default
  - arch:
      - fully_connected
  - scheduler: tf_exponential_lr
  - optimizer: adam
  - 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)
    circle_radius = 0.005 # (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]))
    c1 = Circle(circle1_origin, circle_radius)
    c2 = Circle(circle2_origin, circle_radius)
    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,
        )
    domain.add_constraint(panel_left_1, "panel_left_1")

    # 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,
        )
    domain.add_constraint(panel_right_1, "panel_right_1")

    # 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,
        )
    domain.add_constraint(panel_bottom, "panel_bottom")

    # 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,
            
        )
    domain.add_constraint(panel_top, "panel_top")

    # 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,
    )
    domain.add_constraint(LE, "interior")

    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"],
            # requires_grad=False,
            # 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,
    )
    domain.add_constraint(IC, "IC")

    # 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,
            
        )
    domain.add_constraint(panel_top, "panel_top")
  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,

Hi @nihalpushkar11

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.