[Solver] Two-Phase One-Point Explicit #680

tianchiTJ · 2020-07-29T08:05:28Z

Describe the PR
This PR is to implement the two-phase one-layer explicit solver for fully saturated soil simulation. One particle includes both solid phase and liquid phase, which are defined by the corresponding volume fractions (n_s and n_w). The theory of porous media is adopted to set up the governing equations (momentum balance and mass balance). The pore pressure is updated based on the strain of the solid skeleton, and a weakly compressible fluid assumption is assumed to solve the fluid and mixture balance equations. All implementations work coherently with the current MPI and OpenMP parallelization features.

New implementation
TwoPhase Particle

New class ParticleTwoPhase derived from Particle, which include all necessary functions for two-phase computation and new variables for two-phase (etc. pore pressure, porosity)
New hdf5 structure for ParticleTwoPhase (also includes refactoring of base hdf5 class)

MPMExplicitTwoPhase

New class MPMExplicitTwoPhase derived from MPMBase, which includes the main loop for two-phase computation.
New input functions in MPMBase to read the new types of input files for two-phase.

TwoPhase functions for nodes

New functions for the two-phase are implemented in the current Node class directly (see node_multiphase.tcc)

Free surface detection functions in mesh

New functions to detect the changing surface are implemented in the Mesh class (see mesh_multiphase.tcc)

Unit tests and benchmarks
All new implementations are validated by corresponding test functions, and the unit-test for two-phase computation is also included. Related 2D and 3D one-dimensional consolidation tests can be found here: cb-geo/mpm-benchmarks#19.

…xplicit

…/particles

kks32

Partial review: Could we move the fluid function implementations in mesh to a separate "mesh_fluid.tcc" file so it's easier to manage?

kks32 · 2020-11-03T16:21:22Z

include/cell.h

+
+  //! Map cell volume to the nodes
+  //! \param[in] phase to map volume
+  void map_cell_volume_to_nodes(unsigned phase);


Would this affect GIMP and other non-local algorithms?

Hmm... That's a good point, but I don't think the free-surface detection can be used for non-local MPM yet in general, so I am expecting this function won't work well too for non-local MPM.

include/cell.tcc

include/node.tcc

include/node_base.h

bodhinandach · 2020-11-04T17:06:21Z

@kks32 I take some of the functions out to mesh_multiphase.tcc as you suggested. Will move the other functions related to two-phase out in the semi-implicit and 2P2P too

…/mpm into solver/two-phase-explicit

kks32

Partial review: Review contains comments for mesh amd mesh_mutliphase files.

CMakeLists.txt

kks32 · 2020-11-09T23:27:53Z

include/cell.tcc

+//! Map cell volume to nodes
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::map_cell_volume_to_nodes(unsigned phase) {
+  if (this->status()) {


Why are we doing this only for cells with particles? The function looks generic so adding this for only cells with particles looks strange. Could we instead pass it as a function argument?

We only do this for cells with particles as otherwise, the nodes_->update_volume() will be called as many numbers as the MPI thread. Therefore, without checking the status(), we won't have the free_surface detection working in parallel.

include/mesh.tcc

include/mesh_multiphase.tcc

kks32 · 2020-11-09T23:56:55Z

include/mesh_multiphase.tcc

+
+        //! Check internal cell
+        for (const auto neighbour_cell_id : (*citr)->neighbours()) {
+#if USE_MPI


Why not use the same functions? solving_status?

In that case, we need to assign solving_status for a single MPI threat too, which is just a copy of status

kks32 · 2020-11-09T23:57:55Z

include/mesh_multiphase.tcc

+                // Assign free surface nodes
+                if (!common_node_id.empty()) {
+                  for (const auto common_id : common_node_id) {
+                    map_nodes_[common_id]->assign_free_surface(true);


This needs to be reduced across MPI tasks

Good point! I think, since we did the solving_status all_reduce prior to this call, we don't have to do so.

include/mesh_multiphase.tcc

kks32 · 2020-11-09T23:59:47Z

include/mesh_multiphase.tcc

+  bool status = true;
+
+  try {
+    if (!particles_.size())


Use assert to check this condition

Actually, we just followed the format of other assign_particles_xxx functions. If we change this one, will that be better to change all other functions in mesh.tcc?

kks32 · 2020-11-10T00:01:04Z

include/mesh_multiphase.tcc

+                                             mpm::ParticlePhase::Liquid);
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());


Are we OK if this fails and we can continue the calculation with default values?

I think that would be better if we stop the simulation.

kks32

Could you make sure all return functions are marked const

Partial review

kks32 · 2020-11-14T12:34:25Z

include/node.h

@@ -245,6 +260,21 @@ class Node : public NodeBase<Tdim> {
  //! Set ghost id
  void ghost_id(Index gid) override { ghost_id_ = gid; }

+  //! Return real density at a given node for a given phase


What's a real density? is it the interpolated density of the particles?

Please make sure the return statements are all const

Yes, it is the interpolated particle density, but not from mass_density. It is computed by the division of mass_/volume_ inside nodes. Those two variables are the interpolated variables from particles. I think I will change the detail to "approximated" density instead of "real" density.

kks32 · 2020-11-14T12:40:44Z

include/node_twophase.tcc

+
+    // Apply velocity constraints, which also sets acceleration to 0,
+    // when velocity is set.
+    this->apply_velocity_constraints();


Does this apply velocity constraints to both solid and fluid phases?

Yes, according to this:

mpm/include/node.tcc

Lines 362 to 367 in 689d9fb

// Direction value in the constraint (0, Dim * Nphases)

const unsigned dir = constraint.first;

// Direction: dir % Tdim (modulus)

const auto direction = static_cast<unsigned>(dir % Tdim);

// Phase: Integer value of division (dir / Tdim)

const auto phase = static_cast<unsigned>(dir / Tdim);

include/particles/particle.h

kks32 · 2020-11-14T12:45:49Z

include/particles/particle.h

-      const mpm::dense_map& state_vars,
-      const std::shared_ptr<mpm::Material<Tdim>>& material,
-      unsigned phase = mpm::ParticlePhase::Solid) override;
+      PODParticle& particle,


could we pass it as a const&?

I checked this one, it is not possible as we reinterpret_cast it inside the derived function:

mpm/include/particles/particle_twophase.tcc

Lines 167 to 170 in 689d9fb

bool mpm::TwoPhaseParticle<Tdim>::initialise_particle(PODParticle& particle) {

// Initialise solid phase

bool status = mpm::Particle<Tdim>::initialise_particle(particle);

auto twophase_particle = reinterpret_cast<PODParticleTwoPhase*>(&particle);

mpm/include/particles/particle_twophase.tcc

Lines 199 to 204 in 689d9fb

template <unsigned Tdim>

bool mpm::TwoPhaseParticle<Tdim>::initialise_particle(

PODParticle& particle,

const std::vector<std::shared_ptr<mpm::Material<Tdim>>>& materials) {

auto twophase_particle = reinterpret_cast<PODParticleTwoPhase*>(&particle);

bool status = this->initialise_particle(*twophase_particle);

include/particles/particle.h

kks32 · 2020-11-14T12:47:15Z

include/particles/particle.h

+  //! Free surface
+  bool free_surface_{false};
+  //! Free surface
+  Eigen::Matrix<double, Tdim, 1> normal_;


Specifically a free_surface_normal_?

Do I need to change the name of the assign and return functions too?

kks32 · 2020-11-14T12:52:49Z

include/particles/particle.tcc

@@ -284,6 +290,16 @@ bool mpm::Particle<Tdim>::assign_material_state_vars(
  return status;
 }

+//! Assign a state variable
+template <unsigned Tdim>
+void mpm::Particle<Tdim>::assign_state_variable(const std::string& var,


Why do we need this function to assign a specific state variable?

We should check this in assert, and go ahead without throwing an error. No reason to throw, if we can't recover.

I think we need this to assign_pressure or pore_pressure, just to make it generalize as when we call pressure(), we actually calling state_variable("pressure").

I agree with you to make it in assert instead of throwing an error.

kks32 · 2020-11-14T13:01:29Z

include/particles/particle_base.h

+  //! Update porosity
+  //! \param[in] dt Analysis time step
+  virtual void update_porosity(double dt) {
+    throw std::runtime_error(


Why do we have to do this? I would keep this pure virtual and throw error elsewhere?

That means we need to populate the particle.h with all these functions. Are you okay with that?

kks32 · 2020-11-14T13:04:20Z

include/solvers/mpm_base.tcc

+        // Pore pressure constraint phase indice
+        unsigned constraint_phase = constraints["phase_id"];
+
+        if (constraint_phase >= Tnphases)


We can't check phases in MPM class, we should check it in the nodes

Super cool. I removed this check. It will be automatically called as the constraints_->assign_nodal_pressure_constraints() function is called. In short, they will reach this assertion.

mpm/include/node.tcc

Lines 178 to 186 in 689d9fb

//! Assign pressure constraint

template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>

bool mpm::Node<Tdim, Tdof, Tnphases>::assign_pressure_constraint(

const unsigned phase, const double pressure,

const std::shared_ptr<FunctionBase>& function) {

bool status = true;

try {

// Constrain directions can take values between 0 and Tnphases

assert(phase < Tnphases * 2);

src/particle.cc

bodhinandach · 2020-11-17T12:01:09Z

@kks32 Thanks for your time reviewing this PR. It looks much better now. I partially modified the PR according to your comments. Can you spend some time to check it again?

bodhinandach · 2020-12-03T13:25:49Z

ping @kks32, any further issues with this branch. Other than dynamic load balancing and resume feature in MPI, I am good.

bodhinandach · 2020-12-08T16:34:38Z

@kks32 Some missing features which are not available in MPI while running 2Phase:

compute_free_surface_by_geometry: Please use compute_free_surface_by_density instead. You can change the setting below from the input json:

"free_surface_detection": {
	"type": "density",
	"volume_tolerance": 0.25
}

dynamic_load_balancing and resume functions are not working properly in 3D (It was checked to be OK for 2D). Problem was observed to be occurring in the halo exchange function for shared nodes between MPI rank.
Error in halo exchange:

NOTE: For instance, at node 1050 (a shared node), solid and liquid mass are expected to be 0.01456 and 0.0024, respectively (they are not properly summed). Thus, producing this:

Error in dynamic load balancing at step 10:

tianchiTJ and others added 30 commits June 28, 2020 22:08

Add twophase particle functions

bbf7c51

💻 Refactor particle to handle scalar and vector properties

f3f6fb2

🎯 Add testing to particle refactor of scalar properties

ad09dfa

Merge remote-tracking branch 'origin/develop' into solver/two-phase-e…

ac84455

…xplicit

Add explicit twophase solver

26e5ed9

💻 Compute mass free function

e68ab89

Merge branch 'develop' of https://github.com/cb-geo/mpm into refactor…

d40b6b0

…/particles

🔧 modify mass, mass_density_volume

67303d8

🔧 modify compute_mass function

0a6ac16

🔧 refactor particle velocity and displacement

adfd93f

🔧 refactor nodal mass and volume

2b92b3f

🔧 refactor nodal velocity and acceleration

336cfe6

🔧 refactor nodal momentum, external, and internal forces

91d5c97

🔧 refactor nodal pressure

8ed2eb8

🔧 Add scalar and vector properties enum

02911bf

🔧 aesthetic clean up order of initialise()

0feaf6d

🔧 fix particle vector_properties_ size order

bcfb878

🎯 add testings for update and return functions

99b9eea

Merge branch 'develop' of https://github.com/cb-geo/mpm into refactor…

315cfae

…/particles

🔧 add overload functions of mapping

2f54b90

🚧 🎯 test nightly benchmarks with refactored particles and nodes

0facda0

🚧 🎯 refactor map_mass_momentum and tests

4136654

Merge branch 'develop' of https://github.com/cb-geo/mpm into refactor…

9da224f

…/particles

🚧 🎯 refactor map_body_force and update_volume

cb71c90

🚧 🎯 change assign_velocity to void and cleanup

52eb2b3

🎯 revert mass change in particle test

a20a2ab

🚧 🎯 add return cell pointer and refactor compute_volume

0a0cf7b

🔧 add assert in return and update scalar and vector properties

fe0f9bf

🚧 🎯 return set traction and refactor map traction

b693709

🚧 🎯 refactor pressure smoothing

5fc100e

tianchiTJ and others added 4 commits October 23, 2020 13:26

Add test file for read particle scalar properties

c2b48f2

Add test function for read_pressure_constraints

9f3cc4a

Add test functions for pressure constraint

f6022df

🚧 🎯 modify write mesh particle

1b34a27

kks32 requested changes Nov 3, 2020

View reviewed changes

tianchiTJ and others added 6 commits November 3, 2020 11:29

Change the density tolerance

ad851bf

🚧 correct some suggestions

ae2b6fe

🚧 move multiphase functions outside mesh.tcc

b8a7d89

🚧 add groups in cell and node functions

ae92ec0

🎯 modify test to fix error in assertion

62e065c

🚧 🎯 change nPhase to NPhase in nodes

32f1c5a

tianchiTJ added 2 commits November 4, 2020 09:48

Correct the names of nodal phases

2b195ee

Merge branch 'solver/two-phase-explicit' of https://github.com/cb-geo…

689d9fb

…/mpm into solver/two-phase-explicit

kks32 requested changes Nov 10, 2020

View reviewed changes

kks32 requested changes Nov 14, 2020

View reviewed changes

kks32 reviewed Nov 16, 2020

View reviewed changes

src/particle.cc Outdated Show resolved Hide resolved

bodhinandach added 4 commits November 17, 2020 15:29

🚧 cleanup as suggested

c24b27e

🚧 use iterate_over function with lambda function input

4f550b1

🚧 modify particle and node functions as suggested

2ddbf23

🚧 additional fixes on assertion and some info

fe725f6

🎯 🚧 remove throw test

eae0fbf

bodhinandach mentioned this pull request Nov 18, 2020

2Phase Consolidation in 2D and 3D cb-geo/mpm-benchmarks#19

Open

bodhinandach added 2 commits November 23, 2020 17:27

🚧 grouping node multiphase

72dbc24

🎯 🚧 fix test error due to initialise twophase

2ad9bc2

bodhinandach mentioned this pull request Nov 27, 2020

Add math function from file #692 #701

Closed

bodhinandach approved these changes Dec 3, 2020

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[Solver] Two-Phase One-Point Explicit #680

[Solver] Two-Phase One-Point Explicit #680

tianchiTJ commented Jul 29, 2020 •

edited by bodhinandach

Loading

kks32 left a comment

kks32 Nov 3, 2020

bodhinandach Nov 4, 2020

bodhinandach commented Nov 4, 2020

kks32 left a comment

kks32 Nov 9, 2020

bodhinandach Nov 17, 2020

kks32 Nov 9, 2020

bodhinandach Nov 17, 2020

kks32 Nov 9, 2020

bodhinandach Nov 17, 2020 •

edited

Loading

kks32 Nov 9, 2020

bodhinandach Nov 17, 2020

kks32 Nov 10, 2020

bodhinandach Nov 17, 2020

kks32 left a comment

kks32 Nov 14, 2020

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020 •

edited

Loading

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020 •

edited

Loading

kks32 Nov 14, 2020

bodhinandach Nov 17, 2020

bodhinandach commented Nov 17, 2020

bodhinandach commented Dec 3, 2020

bodhinandach commented Dec 8, 2020

	// Direction value in the constraint (0, Dim * Nphases)
	const unsigned dir = constraint.first;
	// Direction: dir % Tdim (modulus)
	const auto direction = static_cast<unsigned>(dir % Tdim);
	// Phase: Integer value of division (dir / Tdim)
	const auto phase = static_cast<unsigned>(dir / Tdim);

	bool mpm::TwoPhaseParticle<Tdim>::initialise_particle(PODParticle& particle) {
	// Initialise solid phase
	bool status = mpm::Particle<Tdim>::initialise_particle(particle);
	auto twophase_particle = reinterpret_cast<PODParticleTwoPhase*>(&particle);

	template <unsigned Tdim>
	bool mpm::TwoPhaseParticle<Tdim>::initialise_particle(
	PODParticle& particle,
	const std::vector<std::shared_ptr<mpm::Material<Tdim>>>& materials) {
	auto twophase_particle = reinterpret_cast<PODParticleTwoPhase*>(&particle);
	bool status = this->initialise_particle(*twophase_particle);

	//! Assign pressure constraint
	template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
	bool mpm::Node<Tdim, Tdof, Tnphases>::assign_pressure_constraint(
	const unsigned phase, const double pressure,
	const std::shared_ptr<FunctionBase>& function) {
	bool status = true;
	try {
	// Constrain directions can take values between 0 and Tnphases
	assert(phase < Tnphases * 2);

[Solver] Two-Phase One-Point Explicit #680

Are you sure you want to change the base?

[Solver] Two-Phase One-Point Explicit #680

Conversation

tianchiTJ commented Jul 29, 2020 • edited by bodhinandach Loading

kks32 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bodhinandach commented Nov 4, 2020

kks32 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bodhinandach Nov 17, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kks32 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bodhinandach Nov 17, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bodhinandach Nov 17, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bodhinandach commented Nov 17, 2020

bodhinandach commented Dec 3, 2020

bodhinandach commented Dec 8, 2020

tianchiTJ commented Jul 29, 2020 •

edited by bodhinandach

Loading

bodhinandach Nov 17, 2020 •

edited

Loading

bodhinandach Nov 17, 2020 •

edited

Loading

bodhinandach Nov 17, 2020 •

edited

Loading