Skip to content

Commit

Permalink
Add time dependent BC with linear ramp
Browse files Browse the repository at this point in the history
  • Loading branch information
@streeve
streeve committed Nov 10, 2023
1 parent d47622e commit ab6f3e5
Show file tree
Hide file tree
Showing 4 changed files with 154 additions and 58 deletions.
8 changes: 5 additions & 3 deletions examples/crack_branching.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<CabanaPD::RegionBoundary> 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 )
{
Expand Down
9 changes: 5 additions & 4 deletions examples/read_file.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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<memory_space>(
inputs, particles, force_model, bc, prenotch );
Expand Down
186 changes: 139 additions & 47 deletions src/CabanaPD_Boundary.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,10 +151,19 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
{
double _value;
BCIndexSpace _index_space;

BoundaryCondition( const double value, BCIndexSpace bc_index_space )
double _time_start;
double _time_ramp;
double _time_end;

BoundaryCondition( BCIndexSpace bc_index_space, const double value,
const double start = 0.0,
const double ramp = std::numeric_limits<double>::max(),
const double end = std::numeric_limits<double>::max() )
: _value( value )
, _index_space( bc_index_space )
, _time_start( start )
, _time_ramp( ramp )
, _time_end( end )
{
}

Expand All @@ -166,16 +175,28 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
}

template <class ExecSpace, class FieldType, class PositionType>
void apply( ExecSpace, FieldType& f, PositionType& )
void apply( ExecSpace, FieldType& f, PositionType&, const double t )
{
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> 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;
}
} );
}
};
Expand All @@ -186,11 +207,23 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
{
double _value;
BCIndexSpace _index_space;

BoundaryCondition( const double value, BCIndexSpace bc_index_space )
double _time_start;
double _time_ramp;
double _time_end;

BoundaryCondition( BCIndexSpace bc_index_space, const double value,
const double start = 0.0,
const double ramp = std::numeric_limits<double>::max(),
const double end = std::numeric_limits<double>::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 <class ExecSpace, class Particles>
Expand All @@ -201,16 +234,28 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
}

template <class ExecSpace, class FieldType, class PositionType>
void apply( ExecSpace, FieldType& f, PositionType& )
void apply( ExecSpace, FieldType& f, PositionType&, const double t )
{
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> 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;
}
} );
}
};
Expand All @@ -220,31 +265,56 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ForceSymmetric1dBCTag>
{
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<double>::max(),
const double end = std::numeric_limits<double>::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<double>::max(),
const double end = std::numeric_limits<double>::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 <class ExecSpace, class Particles>
Expand All @@ -255,72 +325,94 @@ struct BoundaryCondition<BCIndexSpace, ForceSymmetric1dBCTag>
}

template <class ExecSpace, class FieldType, class PositionType>
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<ExecSpace> policy( 0, index_space.size() );
auto dim = _dim;
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;
}
} );
}
};

// FIXME: relatively large initial guess for allocation.
template <class BoundaryType, class BCTag, class ExecSpace, class Particles,
class... Args>
auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
std::vector<BoundaryType> planes,
const double value,
const double initial_guess = 0.5 )
auto createBoundaryCondition(
BCTag, ExecSpace exec_space, Particles particles,
std::vector<BoundaryType> planes, const double value,
const double start = 0.0,
const double ramp = std::numeric_limits<double>::max(),
const double end = std::numeric_limits<double>::max(),
const double initial_guess = 0.5 )
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
return BoundaryCondition<bc_index_type, BCTag>( bc_indices, value, start,
ramp, end );
}

// FIXME: relatively large initial guess for allocation.
template <class BoundaryType, class ExecSpace, class Particles>
auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space,
Particles particles,
std::vector<BoundaryType> 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<BoundaryType> planes, const double value, const int dim,
const double center, const double start = 0.0,
const double ramp = std::numeric_limits<double>::max(),
const double end = std::numeric_limits<double>::max(),
const double initial_guess = 0.5 )
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, ForceSymmetric1dBCTag>(
value, bc_indices, dim, center );
bc_indices, value, dim, center, start, ramp, end );
}

// FIXME: relatively large initial guess for allocation.
template <class BoundaryType, class ExecSpace, class Particles>
auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space,
Particles particles,
std::vector<BoundaryType> 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<BoundaryType> 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<double>::max(),
const double end = std::numeric_limits<double>::max(),
const double initial_guess = 0.5 )
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, ForceSymmetric1dBCTag>(
value_top, value_bottom, bc_indices, dim, center );
bc_indices, value_top, value_bottom, dim, center, start, ramp, end );
}
} // namespace CabanaPD

Expand Down
9 changes: 5 additions & 4 deletions src/CabanaPD_Solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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.
Expand Down

0 comments on commit ab6f3e5

Please sign in to comment.