A modified case of 3D collapse model

ZenanH commented 6 months ago

Hi,

I modified the configuration of the case DPmaterial.py in the example folder, below is the new file:

from geotaichi import *

init()

mpm = MPM()

mpm.set_configuration(domain=ti.Vector([0.0750, 0.5625, 0.1250]), 
                      background_damping=0.0, 
                      alphaPIC=0.00, 
                      mapping="MUSL", 
                      stabilize=None,
                      shape_function="GIMP",
                      gauss_number=0)

mpm.set_solver(solver={
                           "Timestep":                   6.495190528383289e-5,
                           "SimulationTime":             1.0,
                           "SaveInterval":               1#0.01
                      })

mpm.memory_allocate(memory={
                                "max_material_number":    1,
                                "max_particle_number":    5.12e5,
                                "max_constraint_number":  {
                                                               "max_velocity_constraint":   2004300,#80935,
                                                               #"max_friction_constraint":   53703
                                                          }
                            })

mpm.add_material(model="DruckerPrager",
                 material={
                               "MaterialID":           1,
                               "Density":              2700.,
                               "YoungModulus":         1e6,
                               "PossionRatio":         0.,
                               "Cohesion":             0.,
                               "Friction":             19.8,
                               "Dilation":             0.,
                               "Tensile":              0.
                 })

mpm.add_element(element={
                             "ElementType":               "R8N3D",
                             "ElementSize":               ti.Vector([0.0025, 0.0025, 0.0025]),
                             "Contact":    {
                                                "ContactDetection":                False
                                           }
                        })

mpm.add_region(region={
                            "Name": "region1",
                            "Type": "Rectangle",
                            "BoundingBoxPoint": ti.Vector([0.005, 0.005, 0.005]),
                            "BoundingBoxSize": ti.Vector([0.05, 0.2, 0.1]),
                            "zdirection": ti.Vector([0., 0., 1.])
                      })

mpm.add_body(body={
                       "Template": {
                                       "RegionName":         "region1",
                                       "nParticlesPerCell":  2,
                                       "BodyID":             0,
                                       "MaterialID":         1,
                                       "ParticleStress": {
                                                              "GravityField":     True,
                                                              "InternalStress":   ti.Vector([-0., -0., -0., 0., 0., 0.]),
                                                              "Traction":         {}
                                                         },
                                       "InitialVelocity":ti.Vector([0, 0, 0]),
                                       "FixVelocity":    ["Free", "Free", "Free"]    

                                   }
                   })

mpm.add_boundary_condition(boundary=[

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., 0., None],
                                             "StartPoint":     [0., 0.005, 0.],
                                             "EndPoint":       [0.0750, 0.005, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., None, None],
                                             "StartPoint":     [0.005, 0., 0.],
                                             "EndPoint":       [0.005, 0.5625, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., None, None],
                                             "StartPoint":     [0.055, 0., 0.],
                                             "EndPoint":       [0.055, 0.5625, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., 0., 0.],
                                             "StartPoint":     [0., 0., 0.005],
                                             "EndPoint":       [0.0750, 0.5625, 0.005]
                                        },
                                    ])

mpm.select_save_data()

mpm.run()

mpm.postprocessing()

I found that:

it seems like the maximum run-out distance (fig1) is a little bit larger than my mpm code (fig2) and the results from other paper (fig3, fig4), do you have any comments about this point?
I'm doing a benchmarks across different code. Upper code on RTX3090 took 219.16s without saving data during the simulation. Is it reasonable? Or maybe there are some performance tricks I can try?

ZenanH commented 6 months ago

Fig1. Result from geotaichi, the farthest particle is around 0.55

Fig2. Result from my code, around 0.45

Fig3. Results from other paper (less than 0.5):

Huang, P., Li, S., Guo, H., and Hao, Z.: Large deformation failure analysis of the soil slope based on the material point method, Comput Geosci, 19, 951–963, https://doi.org/10.1007/s10596-015-9512-9, 2015.

Wyser, E.: Elasto-plastic deformations within a material point framework on modern GPU architectures, 2021.

ZenanH commented 6 months ago

day2024_05_14.log

Yihao-Shi commented 6 months ago

Thanks for the questions,

Appendix C has detailed discussed the reason why there shows insignificant difference on run-out distance. Besides, I noticed that the local damping is set zero in this script, thus it undoubtedly results in a larger run-out distance.
I think it is quite normal, cause there are 500,000 particles in the simulation. I am also curious about the computational efficiency about your own code. Is it more faster?

Yihao-Shi commented 6 months ago

You can try this Python script, which may be better:

from geotaichi import *

init(kernel_profiler=True)

mpm = MPM()

mpm.set_configuration(domain=ti.Vector([0.0750, 0.5625, 0.1250]), 
                      background_damping=0.025, 
                      alphaPIC=0.00, 
                      mapping="USF", 
                      stabilize=None,
                      shape_function="GIMP",
                      gauss_number=0)

mpm.set_solver(solver={
                           "Timestep":                   6.495190528383289e-5,
                           "SimulationTime":             1.0,
                           "SaveInterval":               1#0.01
                      })

mpm.memory_allocate(memory={
                                "max_material_number":    1,
                                "max_particle_number":    5.12e5,
                                "max_constraint_number":  {
                                                               "max_velocity_constraint":   2004300,#80935,
                                                               #"max_friction_constraint":   53703
                                                          }
                            })

mpm.add_material(model="DruckerPrager",
                 material={
                               "MaterialID":           1,
                               "Density":              2700.,
                               "YoungModulus":         8.4e5,
                               "PossionRatio":         0.,
                               "Cohesion":             0.,
                               "Friction":             19.8,
                               "Dilation":             0.,
                               "Tensile":              0.
                 })

mpm.add_element(element={
                             "ElementType":               "R8N3D",
                             "ElementSize":               ti.Vector([0.0025, 0.0025, 0.0025]),
                             "Contact":    {
                                                "ContactDetection":                False
                                           }
                        })

mpm.add_region(region={
                            "Name": "region1",
                            "Type": "Rectangle",
                            "BoundingBoxPoint": ti.Vector([0.005, 0.005, 0.005]),
                            "BoundingBoxSize": ti.Vector([0.05, 0.2, 0.1]),
                            "zdirection": ti.Vector([0., 0., 1.])
                      })

mpm.add_body(body={
                       "Template": {
                                       "RegionName":         "region1",
                                       "nParticlesPerCell":  2,
                                       "BodyID":             0,
                                       "MaterialID":         1,
                                       "ParticleStress": {
                                                              "GravityField":     True,
                                                              "InternalStress":   ti.Vector([-0., -0., -0., 0., 0., 0.]),
                                                              "Traction":         {}
                                                         },
                                       "InitialVelocity":ti.Vector([0, 0, 0]),
                                       "FixVelocity":    ["Free", "Free", "Free"]    

                                   }
                   })

mpm.add_boundary_condition(boundary=[

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., 0., None],
                                             "StartPoint":     [0., 0.005, 0.],
                                             "EndPoint":       [0.0750, 0.005, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., None, None],
                                             "StartPoint":     [0.005, 0., 0.],
                                             "EndPoint":       [0.005, 0.5625, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., None, None],
                                             "StartPoint":     [0.055, 0., 0.],
                                             "EndPoint":       [0.055, 0.5625, 0.1250]
                                        },

                                        {    
                                             "BoundaryType":   "VelocityConstraint",
                                             "Velocity":       [0., 0., 0.],
                                             "StartPoint":     [0., 0., 0.005],
                                             "EndPoint":       [0.0750, 0.5625, 0.005]
                                        },
                                    ])

mpm.select_save_data()

mpm.run()

mpm.postprocessing()

ti.profiler.print_kernel_profiler_info()

ZenanH commented 6 months ago

I agree with your point and the reference you mentioned.

I think there is no local damping for all the figures I showed, and also without locking-free method, and I adopted the same situation in geotaichi (no locking-free, no damping).

Anyway, this model has a higher resolution, maybe already changed the accuracy, its not a big issue😅

Thanks for the code, I will try it~

Yihao-Shi / GeoTaichi

A modified case of 3D collapse model #7