From 3a0923ea5bbc58e4b8627f3edc33bfb942b7ef68 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 28 May 2024 10:33:24 -0400 Subject: [PATCH] Rearrange boundary conditions after solver init; remove zero boundary --- examples/mechanics/crack_branching.cpp | 6 +- examples/mechanics/elastic_wave.cpp | 10 +- examples/mechanics/kalthoff_winkler.cpp | 6 +- .../thermomechanics/thermal_deformation.cpp | 38 +++-- src/CabanaPD_Boundary.hpp | 25 ---- src/CabanaPD_Solver.hpp | 140 ++++++++++++------ 6 files changed, 122 insertions(+), 103 deletions(-) diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index fd29321c..d9ea11e4 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -137,9 +137,9 @@ void crackBranchingExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model, prenotch ); + cabana_pd->init(); + cabana_pd->run( bc ); } // Initialize MPI+Kokkos. diff --git a/examples/mechanics/elastic_wave.cpp b/examples/mechanics/elastic_wave.cpp index d70904a0..27fe55a9 100644 --- a/examples/mechanics/elastic_wave.cpp +++ b/examples/mechanics/elastic_wave.cpp @@ -68,12 +68,6 @@ void elasticWaveExample( const std::string filename ) typename model_type::thermal_type>>( exec_space(), low_corner, high_corner, num_cells, halo_width ); - // ==================================================== - // Boundary conditions - // ==================================================== - auto bc = CabanaPD::createBoundaryCondition( - CabanaPD::ZeroBCTag{} ); - // ==================================================== // Custom particle initialization // ==================================================== @@ -111,8 +105,8 @@ void elasticWaveExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverElastic( - inputs, particles, force_model, bc ); - cabana_pd->init_force(); + inputs, particles, force_model ); + cabana_pd->init(); cabana_pd->run(); // ==================================================== diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index a3b2b904..efb00ec8 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -129,9 +129,9 @@ void kalthoffWinklerExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model, prenotch ); + cabana_pd->init(); + cabana_pd->run( bc ); } // Initialize MPI+Kokkos. diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index e25687d6..64554b72 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -71,18 +71,6 @@ void thermalDeformationExample( const std::string filename ) CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); - // ==================================================== - // Imposed field - // ==================================================== - auto x = particles->sliceReferencePosition(); - auto temp = particles->sliceTemperature(); - const double low_corner_y = low_corner[1]; - auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) - { - temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; - }; - auto body_term = CabanaPD::createBodyTerm( temp_func ); - // ==================================================== // Custom particle initialization // ==================================================== @@ -98,12 +86,30 @@ void thermalDeformationExample( const std::string filename ) *particles, delta, K, alpha, temp0 ); // ==================================================== - // Simulation run + // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverElastic( - inputs, particles, force_model, body_term ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model ); + cabana_pd->init(); + + // ==================================================== + // Imposed field + // ==================================================== + auto x = particles->sliceReferencePosition(); + auto temp = particles->sliceTemperature(); + const double low_corner_y = low_corner[1]; + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + }; + auto body_term = CabanaPD::createBodyTerm( temp_func ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->run( body_term ); } // Initialize MPI+Kokkos. diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index d4ac4627..925edc05 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -131,10 +131,6 @@ struct ForceUpdateBCTag { }; -struct ZeroBCTag -{ -}; - template struct BoundaryCondition { @@ -168,20 +164,6 @@ struct BoundaryCondition } }; -template -struct BoundaryCondition -{ - template - void update( ExecSpace, Particles, RegionBoundary ) - { - } - - template - void apply( ExecSpace, ParticleType, double ) - { - } -}; - template struct BoundaryCondition { @@ -282,13 +264,6 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space, user_functor ); } -template -auto createBoundaryCondition( ZeroBCTag ) -{ - using bc_index_type = BoundaryIndexSpace; - return BoundaryCondition(); -} - } // namespace CabanaPD #endif diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 81b3864a..50952ee6 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -90,7 +90,7 @@ class SolverBase }; template + class ForceModel> class SolverElastic { public: @@ -108,14 +108,12 @@ class SolverElastic Cabana::VerletLayout2D, Cabana::TeamOpTag>; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = InputType; - using bc_type = BoundaryCondition; SolverElastic( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc ) + force_model_type force_model ) : inputs( _inputs ) , particles( _particles ) - , boundary_condition( bc ) { neighbor_time = 0; force_time = 0; @@ -199,7 +197,7 @@ class SolverElastic init_time += init_timer.seconds(); } - void init_force() + void init() { init_timer.reset(); @@ -226,7 +224,8 @@ class SolverElastic init_time += init_timer.seconds(); } - void run() + template + void run( BoundaryType boundary_condition ) { init_output(); @@ -270,26 +269,79 @@ class SolverElastic integrator->finalHalfStep( *particles ); integrate_time += integrate_timer.seconds(); - // Print output. - if ( step % output_frequency == 0 ) - { - energy_timer.reset(); - auto W = computeEnergy( *force, *particles, *neighbors, - neigh_iter_tag() ); - energy_time += energy_timer.seconds(); + output( step ); + } - output_timer.reset(); - step_output( step, W ); - particles->output( step / output_frequency, step * dt, - output_reference ); - output_time += output_timer.seconds(); - } + // Final output and timings. + final_output(); + } + + void run() + { + init_output(); + + // Main timestep loop. + for ( int step = 1; step <= num_steps; step++ ) + { + // Integrate - velocity Verlet first half. + integrate_timer.reset(); + integrator->initialHalfStep( *particles ); + integrate_time += integrate_timer.seconds(); + + // Update ghost particles. + comm_timer.reset(); + comm->gatherDisplacement(); + comm_time += comm_timer.seconds(); + + // Do not need to recompute LPS weighted volume here without damage. + // Compute/communicate LPS dilatation (does nothing for PMB). + force_timer.reset(); + force->computeDilatation( *particles, *neighbors, + neigh_iter_tag{} ); + force_time += force_timer.seconds(); + comm_timer.reset(); + comm->gatherDilatation(); + comm_time += comm_timer.seconds(); + + // Compute internal forces. + force_timer.reset(); + computeForce( *force, *particles, *neighbors, neigh_iter_tag{} ); + force_time += force_timer.seconds(); + + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + + // Integrate - velocity Verlet second half. + integrate_timer.reset(); + integrator->finalHalfStep( *particles ); + integrate_time += integrate_timer.seconds(); + + output( step ); } // Final output and timings. final_output(); } + void output( const int step ) + { + // Print output. + if ( step % output_frequency == 0 ) + { + energy_timer.reset(); + auto W = computeEnergy( *force, *particles, *neighbors, + neigh_iter_tag() ); + energy_time += energy_timer.seconds(); + + output_timer.reset(); + step_output( step, W ); + particles->output( step / output_frequency, step * dt, + output_reference ); + output_time += output_timer.seconds(); + } + } + void init_output() { // Output after construction and initial forces. @@ -360,7 +412,6 @@ class SolverElastic std::shared_ptr integrator; std::shared_ptr force; std::shared_ptr neighbors; - bc_type boundary_condition; std::string output_file; std::string error_file; @@ -386,14 +437,13 @@ class SolverElastic }; template + class ForceModel, class PrenotchType> class SolverFracture - : public SolverElastic + : public SolverElastic { public: - using base_type = SolverElastic; + using base_type = + SolverElastic; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; @@ -404,15 +454,13 @@ class SolverFracture using force_model_type = ForceModel; using force_type = typename base_type::force_type; using neigh_iter_tag = Cabana::SerialOpTag; - using bc_type = BoundaryCondition; using prenotch_type = PrenotchType; using input_type = typename base_type::input_type; SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc, - prenotch_type prenotch ) - : base_type( _inputs, _particles, force_model, bc ) + force_model_type force_model, prenotch_type prenotch ) + : base_type( _inputs, _particles, force_model ) { init_timer.reset(); @@ -426,8 +474,8 @@ class SolverFracture SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc ) - : base_type( _inputs, _particles, force_model, bc ) + force_model_type force_model ) + : base_type( _inputs, _particles, force_model ) { init_timer.reset(); @@ -447,7 +495,7 @@ class SolverFracture Kokkos::deep_copy( mu, 1 ); } - void init_force() + void init() { init_timer.reset(); // Compute/communicate weighted volume for LPS (does nothing for PMB). @@ -461,14 +509,12 @@ class SolverFracture computeForce( *force, *particles, *neighbors, mu, neigh_iter_tag{} ); computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() ); - // Add boundary condition. - boundary_condition.apply( exec_space(), *particles, 0 ); - particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); } - void run() + template + void run( BoundaryType boundary_condition ) { this->init_output(); @@ -540,7 +586,6 @@ class SolverFracture using base_type::output_reference; protected: - using base_type::boundary_condition; using base_type::comm; using base_type::force; using base_type::inputs; @@ -571,26 +616,25 @@ class SolverFracture }; template + class ForceModel> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model, BCType bc ) + ForceModel model ) { - return std::make_shared>( - inputs, particles, model, bc ); + return std::make_shared< + SolverElastic>( + inputs, particles, model ); } template + class ForceModel, class PrenotchType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, BCType bc, PrenotchType prenotch ) + ForceModel model, PrenotchType prenotch ) { - return std::make_shared< - SolverFracture>( inputs, particles, model, bc, - prenotch ); + return std::make_shared>( + inputs, particles, model, prenotch ); } } // namespace CabanaPD