From e51430a2038be6ba40bc59cc8a91d6f90ff0ca49 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 09:39:48 -0500 Subject: [PATCH] Add time dependent BC with linear ramp --- examples/crack_branching.cpp | 8 +- examples/read_file.cpp | 9 +- src/CabanaPD_Boundary.hpp | 182 ++++++++++++++++++++++++++--------- src/CabanaPD_Solver.hpp | 9 +- 4 files changed, 152 insertions(+), 56 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index d71812a4..cc045cc3 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -106,9 +106,11 @@ int main( int argc, char* argv[] ) CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy, high_y + dy, low_z, high_z ); std::vector planes = { plane1, plane2 }; - auto bc = - createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, - exec_space{}, *particles, planes, b0 ); + int bc_dim = 1; + double center = 0.0; + auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, + exec_space{}, *particles, planes, b0, + bc_dim, center ); auto init_functor = KOKKOS_LAMBDA( const int pid ) { diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 39181a30..8d7fe8bb 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -123,7 +123,8 @@ int main( int argc, char* argv[] ) using memory_space = typename exec_space::memory_space; // Time - double t_final = 1.0; + double t_final = 1e-1; + double t_ramp = 1e-3; double dt = 1e-7; double output_frequency = 100; @@ -184,9 +185,9 @@ int main( int argc, char* argv[] ) int bc_dim = 1; double center = particles->local_mesh_ext[bc_dim] / 2.0 + particles->local_mesh_lo[bc_dim]; - auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, - exec_space{}, *particles, planes, - 2e4, 0.0, bc_dim, center ); + auto bc = createBoundaryCondition( + CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, + 2e4, 0.0, bc_dim, center, 0.0, t_ramp ); auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, bc, prenotch ); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 975190ce..0b5392b5 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -151,10 +151,19 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; + double _time_start; + double _time_ramp; + double _time_end; - BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + BoundaryCondition( const double value, BCIndexSpace bc_index_space, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) : _value( value ) , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) { } @@ -166,16 +175,28 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& ) + void apply( ExecSpace, FieldType& f, PositionType&, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - for ( int d = 0; d < 3; d++ ) - f( pid, d ) = value; + double current = value; + // Initial linear ramp. + if ( t < ramp ) + { + current = value * ( t - start ) / ( ramp - start ); + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + for ( int d = 0; d < 3; d++ ) + f( pid, d ) = current; + } } ); } }; @@ -186,11 +207,23 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; + double _time_start; + double _time_ramp; + double _time_end; - BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + BoundaryCondition( const double value, BCIndexSpace bc_index_space, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) : _value( value ) , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); } template @@ -201,16 +234,28 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& ) + void apply( ExecSpace, FieldType& f, PositionType&, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - for ( int d = 0; d < 3; d++ ) - f( pid, d ) += value; + double current = value; + // Initial linear ramp. + if ( t < ramp ) + { + current = value * ( t - start ) / ( ramp - start ); + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + for ( int d = 0; d < 3; d++ ) + f( pid, d ) += current; + } } ); } }; @@ -220,31 +265,56 @@ struct BoundaryCondition template struct BoundaryCondition { + BCIndexSpace _index_space; double _value_top; double _value_bottom; - BCIndexSpace _index_space; + double _time_start; + double _time_ramp; + double _time_end; + int _dim; double _center; - BoundaryCondition( const double value, BCIndexSpace bc_index_space, - const int dim = 1, const double center = 0.0 ) - : _value_top( value ) + BoundaryCondition( BCIndexSpace bc_index_space, const double value, + const int dim, const double center, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) + : _index_space( bc_index_space ) + , _value_top( value ) , _value_bottom( value ) - , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) , _dim( dim ) , _center( center ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); + assert( _dim >= 0 ); + assert( _dim < 3 ); } - BoundaryCondition( const double value_top, const double value_bottom, - BCIndexSpace bc_index_space, const int dim = 1, - const double center = 0.0 ) - : _value_top( value_top ) + BoundaryCondition( BCIndexSpace bc_index_space, const double value_top, + const double value_bottom, const int dim, + const double center, const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) + : _index_space( bc_index_space ) + , _value_top( value_top ) , _value_bottom( value_bottom ) - , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) , _dim( dim ) , _center( center ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); + assert( _dim >= 0 ); + assert( _dim < 3 ); } template @@ -255,7 +325,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& x ) + void apply( ExecSpace, FieldType& f, PositionType& x, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); @@ -263,13 +333,28 @@ struct BoundaryCondition auto center = _center; auto value_top = _value_top; auto value_bottom = _value_bottom; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - if ( x( pid, dim ) > center ) - f( pid, dim ) += value_top; - else - f( pid, dim ) += value_bottom; + auto current_top = value_top; + auto current_bottom = value_bottom; + // Initial linear ramp. + if ( t < ramp ) + { + auto t_factor = ( t - start ) / ( ramp - start ); + current_top = value_top * t_factor; + current_bottom = value_bottom * t_factor; + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + if ( x( pid, dim ) > center ) + f( pid, dim ) += current_top; + else + f( pid, dim ) += current_bottom; + } } ); } }; @@ -277,50 +362,57 @@ struct BoundaryCondition // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles, - std::vector planes, - const double value, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + BCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); - return BoundaryCondition( value, bc_indices ); + return BoundaryCondition( bc_indices, value, start, + ramp, end ); } // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, - Particles particles, - std::vector planes, - const double value, const int dim, - const double center, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value, const int dim, + const double center, const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( - value, bc_indices, dim, center ); + bc_indices, value, dim, center, start, ramp, end ); } // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, - Particles particles, - std::vector planes, - const double value_top, const double value_bottom, - const int dim, const double center, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value_top, + const double value_bottom, const int dim, const double center, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( - value_top, value_bottom, bc_indices, dim, center ); + bc_indices, value_top, value_bottom, dim, center, start, ramp, end ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 2157db8c..800e78d3 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -398,12 +398,12 @@ class SolverFracture if ( force_bc ) { auto f = particles->sliceForce(); - boundary_condition.apply( exec_space(), f, x ); + boundary_condition.apply( exec_space(), f, x, 0.0 ); } else { auto u = particles->sliceDisplacement(); - boundary_condition.apply( exec_space(), u, x ); + boundary_condition.apply( exec_space(), u, x, 0.0 ); } particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); @@ -450,15 +450,16 @@ class SolverFracture // Add boundary condition. auto x = particles->sliceReferencePosition(); // FIXME: this needs to be generalized to any field. + auto time = step * inputs->timestep; if ( force_bc ) { auto f = particles->sliceForce(); - boundary_condition.apply( exec_space(), f, x ); + boundary_condition.apply( exec_space(), f, x, time ); } else { auto u = particles->sliceDisplacement(); - boundary_condition.apply( exec_space(), u, x ); + boundary_condition.apply( exec_space(), u, x, time ); } // Integrate - velocity Verlet second half.