Merge pull request #132 from streeve/fixup_bc_order

fixup: BC created before ghost particles are incorrect
ORNL · Sep 24, 2024 · d53312a · d53312a
2 parents 67109b9 + 0343d8d
commit d53312a
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 30 deletions.
diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp
@@ -78,12 +78,9 @@ void crackBranchingExample( const std::string filename )
         exec_space{}, inputs );
 
     // ====================================================
-    //                Boundary conditions
+    //                Boundary conditions planes
     // ====================================================
-    double sigma0 = inputs["traction"];
     double dy = particles->dx[1];
-    double b0 = sigma0 / dy;
-
     CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane1(
         low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy,
         low_corner[2], high_corner[2] );
@@ -92,21 +89,6 @@ void crackBranchingExample( const std::string filename )
         low_corner[2], high_corner[2] );
     std::vector<CabanaPD::RegionBoundary<CabanaPD::RectangularPrism>> planes = {
         plane1, plane2 };
-    auto particles_f = particles->getForce();
-    auto particles_x = particles->getReferencePosition();
-    // Create a symmetric force BC in the y-direction.
-    auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
-    {
-        // Get a modifiable copy of force.
-        auto p_f = particles_f.getParticleView( pid );
-        // Get a copy of the position.
-        auto p_x = particles_x.getParticle( pid );
-        auto ypos = Cabana::get( p_x, CabanaPD::Field::ReferencePosition(), 1 );
-        auto sign = std::abs( ypos ) / ypos;
-        Cabana::get( p_f, CabanaPD::Field::Force(), 1 ) += b0 * sign;
-    };
-    auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes,
-                                       true );
 
     // ====================================================
     //            Custom particle initialization
@@ -134,6 +116,24 @@ void crackBranchingExample( const std::string filename )
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
         inputs, particles, force_model, prenotch );
 
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+    // Create BC last to ensure ghost particles are included.
+    double sigma0 = inputs["traction"];
+    double b0 = sigma0 / dy;
+    f = particles->sliceForce();
+    x = particles->sliceReferencePosition();
+    // Create a symmetric force BC in the y-direction.
+    auto bc_op = KOKKOS_LAMBDA( const int pid, const double )
+    {
+        auto ypos = x( pid, 1 );
+        auto sign = std::abs( ypos ) / ypos;
+        f( pid, 1 ) += b0 * sign;
+    };
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes,
+                                       true );
+
     // ====================================================
     //                   Simulation run
     // ====================================================

diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp
@@ -91,17 +91,6 @@ void kalthoffWinklerExample( const std::string filename )
     auto particles = CabanaPD::createParticles<memory_space, model_type>(
         exec_space(), low_corner, high_corner, num_cells, halo_width );
 
-    // ====================================================
-    //                Boundary conditions
-    // ====================================================
-    double dx = particles->dx[0];
-    double x_bc = -0.5 * height;
-    CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
-        x_bc - dx, x_bc + dx, y_prenotch1 - 0.25 * dx, y_prenotch2 + 0.25 * dx,
-        -thickness, thickness );
-    auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0,
-                                       exec_space{}, *particles, plane );
-
     // ====================================================
     //            Custom particle initialization
     // ====================================================
@@ -110,6 +99,7 @@ void kalthoffWinklerExample( const std::string filename )
     auto v = particles->sliceVelocity();
     auto f = particles->sliceForce();
 
+    double dx = particles->dx[0];
     double v0 = inputs["impactor_velocity"];
     auto init_functor = KOKKOS_LAMBDA( const int pid )
     {
@@ -128,6 +118,17 @@ void kalthoffWinklerExample( const std::string filename )
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
         inputs, particles, force_model, prenotch );
 
+    // ====================================================
+    //                Boundary conditions
+    // ====================================================
+    // Create BC last to ensure ghost particles are included.
+    double x_bc = -0.5 * height;
+    CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
+        x_bc - dx, x_bc + dx, y_prenotch1 - 0.25 * dx, y_prenotch2 + 0.25 * dx,
+        -thickness, thickness );
+    auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0,
+                                       exec_space{}, *particles, plane );
+
     // ====================================================
     //                   Simulation run
     // ====================================================