From 788ef8de3d26b1b538a158d9c96822da4720d7a1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 20 Oct 2023 11:43:08 -0400 Subject: [PATCH 01/15] Add temperature field --- src/CabanaPD_Particles.hpp | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 34a70c28..23114a78 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -100,9 +100,9 @@ class Particles using scalar_type = Cabana::MemberTypes; // no-fail. using int_type = Cabana::MemberTypes; - // type, W, v, rho, damage. + // v, W, rho, damage, temperature, type. using other_types = - Cabana::MemberTypes; + Cabana::MemberTypes; // Potentially needed later: body force (b), ID. // FIXME: add vector length. @@ -252,6 +252,7 @@ class Particles auto y = sliceCurrentPosition(); auto vol = sliceVolume(); auto nofail = sliceNoFail(); + auto temp = sliceTemperature(); // Initialize particles. auto create_functor = @@ -280,6 +281,7 @@ class Particles type( pid ) = 0; nofail( pid ) = 0; rho( pid ) = 1.0; + temp( pid ) = 0.0; return create; }; @@ -345,8 +347,8 @@ class Particles { return Cabana::slice<0>( _aosoa_vol, "volume" ); } - auto sliceType() { return Cabana::slice<0>( _aosoa_other, "type" ); } - auto sliceType() const { return Cabana::slice<0>( _aosoa_other, "type" ); } + auto sliceType() { return Cabana::slice<5>( _aosoa_other, "type" ); } + auto sliceType() const { return Cabana::slice<5>( _aosoa_other, "type" ); } auto sliceStrainEnergy() { return Cabana::slice<1>( _aosoa_other, "strain_energy" ); @@ -357,21 +359,29 @@ class Particles } auto sliceVelocity() { - return Cabana::slice<2>( _aosoa_other, "velocities" ); + return Cabana::slice<0>( _aosoa_other, "velocities" ); } auto sliceVelocity() const { - return Cabana::slice<2>( _aosoa_other, "velocities" ); + return Cabana::slice<0>( _aosoa_other, "velocities" ); } - auto sliceDensity() { return Cabana::slice<3>( _aosoa_other, "density" ); } + auto sliceDensity() { return Cabana::slice<2>( _aosoa_other, "density" ); } auto sliceDensity() const { - return Cabana::slice<3>( _aosoa_other, "density" ); + return Cabana::slice<2>( _aosoa_other, "density" ); } - auto sliceDamage() { return Cabana::slice<4>( _aosoa_other, "damage" ); } + auto sliceDamage() { return Cabana::slice<3>( _aosoa_other, "damage" ); } auto sliceDamage() const { - return Cabana::slice<4>( _aosoa_other, "damage" ); + return Cabana::slice<3>( _aosoa_other, "damage" ); + } + auto sliceTemperature() + { + return Cabana::slice<4>( _aosoa_other, "temperature" ); + } + auto sliceTemperature() const + { + return Cabana::slice<4>( _aosoa_other, "temperature" ); } auto sliceNoFail() { From 430a8d44ec0f7919bf5edb2b707d4887b0a39330 Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Mon, 20 Nov 2023 17:16:54 -0500 Subject: [PATCH 02/15] Initial thermo-mechanics version --- examples/CMakeLists.txt | 5 +- examples/thermal_deformation.cpp | 159 +++++++++++++++++++++++++++++++ src/CabanaPD_Force.hpp | 36 +++++-- 3 files changed, 191 insertions(+), 9 deletions(-) create mode 100644 examples/thermal_deformation.cpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index df1b8e15..f45e6dad 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -7,4 +7,7 @@ target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD) add_executable(CrackBranching crack_branching.cpp) target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(ThermalDeformation thermal_deformation.cpp) +target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching ThermalDeformation DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp new file mode 100644 index 00000000..e53cd452 --- /dev/null +++ b/examples/thermal_deformation.cpp @@ -0,0 +1,159 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + + { + Kokkos::ScopeGuard scope_guard( argc, argv ); + + // FIXME: change backend at compile time for now. + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + std::array num_cell = { 100, 30, 3 }; + std::array low_corner = { -0.5, 0.0, -0.015 }; + std::array high_corner = { 0.5, 0.3, 0.015 }; + double t_final = 0.0093; + double dt = 7.5E-7; + int output_frequency = 5; + bool output_reference = true; + + double rho0 = 3980; // [kg/m^3] + double E = 370e+9; // [Pa] + double nu = 0.25; // unitless + double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] + double alpha = 7.5e-6; // [1/oC] + double delta = 0.03; + + // Reference temperature + double temp0 = 0.0; + + int m = std::floor( + delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); + int halo_width = m + 1; // Just to be safe. + + // Choose force model type. + using model_type = + CabanaPD::ForceModel; + // model_type force_model( delta, K ); + model_type force_model( delta, K, alpha ); + // using model_type = + // CabanaPD::ForceModel; + // model_type force_model( delta, K, G ); + + CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final, + dt, output_frequency, output_reference ); + inputs.read_args( argc, argv ); + + // Create particles from mesh. + // Does not set displacements, velocities, etc. + auto particles = std::make_shared< + CabanaPD::Particles>( + exec_space(), inputs.low_corner, inputs.high_corner, + inputs.num_cells, halo_width ); + + // Define particle initialization. + auto x = particles->sliceReferencePosition(); + auto u = particles->sliceDisplacement(); + auto v = particles->sliceVelocity(); + auto rho = particles->sliceDensity(); + + /* + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + double a = 0.001; + double r0 = 0.25; + double l = 0.07; + double norm = std::sqrt( x( pid, 0 ) * x( pid, 0 ) + + x( pid, 1 ) * x( pid, 1 ) + + x( pid, 2 ) * x( pid, 2 ) ); + double diff = norm - r0; + double arg = diff * diff / l / l; + for ( int d = 0; d < 3; d++ ) + { + double comp = 0.0; + if ( norm > 0.0 ) + comp = x( pid, d ) / norm; + u( pid, d ) = a * std::exp( -arg ) * comp; + v( pid, d ) = 0.0; + } + rho( pid ) = rho0; + }; + + */ + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + rho( pid ) = rho0; + }; + particles->updateParticles( exec_space{}, init_functor ); + + auto cabana_pd = CabanaPD::createSolverElastic( + inputs, particles, force_model ); + cabana_pd->init_force(); + cabana_pd->run(); + + x = particles->sliceReferencePosition(); + u = particles->sliceDisplacement(); + double num_cell_x = inputs.num_cells[0]; + auto profile = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "displacement_profile" ), + num_cell_x ); + int mpi_rank; + MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); + Kokkos::View count( "c", 1 ); + // double dx = particles->dx[0]; + + /* + auto measure_profile = KOKKOS_LAMBDA( const int pid ) + { + if ( x( pid, 1 ) < dx / 2.0 && x( pid, 1 ) > -dx / 2.0 && + x( pid, 2 ) < dx / 2.0 && x( pid, 2 ) > -dx / 2.0 ) + { + auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); + profile( c, 0 ) = x( pid, 0 ); + profile( c, 1 ) = u( pid, 0 ); + } + }; + */ + + /* + Kokkos::RangePolicy policy( 0, x.size() ); + Kokkos::parallel_for( "displacement_profile", policy, measure_profile + ); auto count_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace{}, count ); auto profile_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, profile + ); + */ + /* + std::fstream fout; + std::string file_name = "displacement_profile.txt"; + fout.open( file_name, std::ios::app ); + for ( int p = 0; p < count_host( 0 ); p++ ) + { + fout << mpi_rank << " " << profile_host( p, 0 ) << " " + << profile_host( p, 1 ) << std::endl; + } + */ + } + + MPI_Finalize(); +} diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 13f52459..ade17bf2 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -223,12 +223,13 @@ struct ForceModel : public BaseForceModel double c; double K; + double alpha; ForceModel(){}; - ForceModel( const double delta, const double K ) + ForceModel( const double delta, const double K, const double alpha ) : base_type( delta ) { - set_param( delta, K ); + set_param( delta, K, alpha ); } ForceModel( const ForceModel& model ) @@ -236,12 +237,14 @@ struct ForceModel : public BaseForceModel { c = model.c; K = model.K; + alpha = model.alpha; } - void set_param( const double _delta, const double _K ) + void set_param( const double _delta, const double _K, const double _alpha ) { delta = _delta; K = _K; + alpha = _alpha; c = 18.0 * K / ( 3.141592653589793 * delta * delta * delta * delta ); } }; @@ -253,6 +256,7 @@ struct ForceModel : public ForceModel using base_model = PMB; using fracture_type = Fracture; + using base_type::alpha; using base_type::c; using base_type::delta; using base_type::K; @@ -261,10 +265,11 @@ struct ForceModel : public ForceModel double bond_break_coeff; ForceModel() {} - ForceModel( const double delta, const double K, const double G0 ) - : base_type( delta, K ) + ForceModel( const double delta, const double K, const double alpha, + const double G0 ) + : base_type( delta, K, alpha ) { - set_param( delta, K, G0 ); + set_param( delta, K, alpha, G0 ); } ForceModel( const ForceModel& model ) @@ -275,9 +280,10 @@ struct ForceModel : public ForceModel bond_break_coeff = model.bond_break_coeff; } - void set_param( const double _delta, const double _K, const double _G0 ) + void set_param( const double _delta, const double _K, const double _alpha, + const double _G0 ) { - base_type::set_param( _delta, _K ); + base_type::set_param( _delta, _K, _alpha ); G0 = _G0; s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); @@ -293,6 +299,7 @@ struct ForceModel : public ForceModel using base_type::base_type; + using base_type::alpha; using base_type::c; using base_type::delta; using base_type::K; @@ -307,6 +314,7 @@ struct ForceModel : public ForceModel using base_type::base_type; + using base_type::alpha; using base_type::c; using base_type::delta; using base_type::K; @@ -925,7 +933,9 @@ class Force> ParallelType& neigh_op_tag ) const { auto c = _model.c; + auto alpha = _model.alpha; const auto vol = particles.sliceVolume(); + const auto temp = particles.sliceTemperature(); auto force_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -936,6 +946,16 @@ class Force> double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Compute average temperature + // It should be: + // double T_av = 0.5*( temp( i ) + temp( j ) ) - temp0; + // assume temp0 = 0 for now. + double T_av = 0.5 * ( temp( i ) + temp( j ) ); + + // Add thermal stretch + s -= alpha * T_av; + const double coeff = c * s * vol( j ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; From c8ee37e393347df1578f7a617402073ab886cb6f Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Wed, 22 Nov 2023 16:37:07 -0500 Subject: [PATCH 03/15] Initial thermo-mechanics version: constant T profile --- examples/crack_branching.cpp | 3 +++ examples/kalthoff_winkler.cpp | 7 ++----- examples/thermal_deformation.cpp | 6 ++++-- src/CabanaPD_Force.hpp | 3 ++- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index 206bf8fc..cc20169d 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -63,6 +63,9 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions ); // Choose force model type. + // using model_type = + // CabanaPD::ForceModel; + // model_type force_model( delta, K, G0 ); using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); diff --git a/examples/kalthoff_winkler.cpp b/examples/kalthoff_winkler.cpp index f5d8e320..a783a91a 100644 --- a/examples/kalthoff_winkler.cpp +++ b/examples/kalthoff_winkler.cpp @@ -73,11 +73,8 @@ int main( int argc, char* argv[] ) // Choose force model type. using model_type = - CabanaPD::ForceModel; - model_type force_model( delta, K, G0 ); - // using model_type = - // CabanaPD::ForceModel; - // model_type force_model( delta, K, G, G0 ); + CabanaPD::ForceModel; + model_type force_model( delta, K, G, G0 ); // Create particles from mesh. // Does not set displacements, velocities, etc. diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index e53cd452..1e0c7f8b 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -34,7 +34,7 @@ int main( int argc, char* argv[] ) std::array high_corner = { 0.5, 0.3, 0.015 }; double t_final = 0.0093; double dt = 7.5E-7; - int output_frequency = 5; + int output_frequency = 200; bool output_reference = true; double rho0 = 3980; // [kg/m^3] @@ -45,7 +45,7 @@ int main( int argc, char* argv[] ) double delta = 0.03; // Reference temperature - double temp0 = 0.0; + // double temp0 = 0.0; int m = std::floor( delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); @@ -76,6 +76,7 @@ int main( int argc, char* argv[] ) auto u = particles->sliceDisplacement(); auto v = particles->sliceVelocity(); auto rho = particles->sliceDensity(); + auto temp = particles->sliceTemperature(); /* auto init_functor = KOKKOS_LAMBDA( const int pid ) @@ -102,6 +103,7 @@ int main( int argc, char* argv[] ) */ auto init_functor = KOKKOS_LAMBDA( const int pid ) { + temp( pid ) = 5000 * x( pid, 1 ); rho( pid ) = rho0; }; particles->updateParticles( exec_space{}, init_functor ); diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index ade17bf2..4b90970e 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -950,7 +950,8 @@ class Force> // Compute average temperature // It should be: // double T_av = 0.5*( temp( i ) + temp( j ) ) - temp0; - // assume temp0 = 0 for now. + + // assume temp0 = 0 for now double T_av = 0.5 * ( temp( i ) + temp( j ) ); // Add thermal stretch From 8848f957bafbc64783aaf459f38144488398dfd6 Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Thu, 7 Dec 2023 16:01:49 -0500 Subject: [PATCH 04/15] Changes for thermal deformation example --- src/CabanaPD_Boundary.hpp | 7 ++++--- src/CabanaPD_Force.hpp | 15 ++++++++++++--- src/CabanaPD_Solver.hpp | 5 +++-- 3 files changed, 19 insertions(+), 8 deletions(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 6a792f2a..9a7fc4bf 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -155,15 +155,16 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& ) + void apply( ExecSpace, ParticleType& particles, double t ) { - auto user = _user_functor; + auto temp = particles.sliceTemperature(); + auto x = particles.sliceReferencePosition(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); - user( pid ); + temp( pid ) = 5000 * x( pid, 1 ) * t; } ); } }; diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 4b90970e..c8246090 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -1035,8 +1035,11 @@ class Force> const int n_local, ParallelType& ) const { auto c = _model.c; - auto break_coeff = _model.bond_break_coeff; + auto alpha = _model.alpha; + // auto break_coeff = _model.bond_break_coeff; + auto break_coeff = _model.s0; const auto vol = particles.sliceVolume(); + const auto temp = particles.sliceTemperature(); const auto nofail = particles.sliceNoFail(); auto force_full = KOKKOS_LAMBDA( const int i ) @@ -1059,9 +1062,15 @@ class Force> double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + // assume temp0 = 0 for now + double T_av = 0.5 * ( temp( i ) + temp( j ) ); + // Break if beyond critical stretch unless in no-fail zone. - if ( r * r >= break_coeff * xi * xi && !nofail( i ) && - !nofail( j ) ) + if ( r * r >= ( 1.0 + s0 + alpha * T_av ) * + ( 1.0 + s0 + alpha * T_av ) * xi * xi && + !nofail( i ) && !nofail( j ) ) + // if ( r * r >= break_coeff * xi * xi && !nofail( i ) && + // !nofail( j ) ) { mu( i, n ) = 0; } diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index c1e48135..c71c7ef1 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -448,7 +448,7 @@ class SolverFracture computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, 0 ); particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); @@ -493,7 +493,8 @@ class SolverFracture force_time += force_timer.seconds(); // Add boundary condition. - boundary_condition.apply( exec_space{}, *particles ); + boundary_condition.apply( exec_space{}, *particles, + step * inputs->timestep ); // Integrate - velocity Verlet second half. integrate_timer.reset(); From 152ea83fde372bc28b4eb0f4de47703891b6a144 Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Tue, 12 Dec 2023 16:16:18 -0500 Subject: [PATCH 05/15] Changes for time-dependent thermal deformation --- examples/inputs/kalthoff_winkler_m3.json | 13 +++++ examples/inputs/kalthoff_winkler_m4.json | 13 +++++ examples/inputs/kalthoff_winkler_m6.json | 13 +++++ examples/thermal_deformation.cpp | 62 ++---------------------- src/CabanaPD_Boundary.hpp | 2 +- src/CabanaPD_Force.hpp | 2 +- 6 files changed, 44 insertions(+), 61 deletions(-) create mode 100644 examples/inputs/kalthoff_winkler_m3.json create mode 100644 examples/inputs/kalthoff_winkler_m4.json create mode 100644 examples/inputs/kalthoff_winkler_m6.json diff --git a/examples/inputs/kalthoff_winkler_m3.json b/examples/inputs/kalthoff_winkler_m3.json new file mode 100644 index 00000000..c79728d4 --- /dev/null +++ b/examples/inputs/kalthoff_winkler_m3.json @@ -0,0 +1,13 @@ +{ + "num_cells" : {"value": [151, 301, 14]}, + "system_size" : {"value": [0.1, 0.2, 0.009], "unit": "m"}, + "density" : {"value": 8000, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 191e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 42408, "unit": "J/m^2"}, + "horizon" : {"value": 0.002, "unit": "m"}, + "impactor_velocity" : {"value": 16, "unit": "m/sec"}, + "final_time" : {"value": 70e-6, "unit": "s"}, + "timestep" : {"value": 0.133e-6, "unit": "s"}, + "output_frequency" : {"value": 500}, + "output_reference" : {"value": true} +} diff --git a/examples/inputs/kalthoff_winkler_m4.json b/examples/inputs/kalthoff_winkler_m4.json new file mode 100644 index 00000000..791625cf --- /dev/null +++ b/examples/inputs/kalthoff_winkler_m4.json @@ -0,0 +1,13 @@ +{ + "num_cells" : {"value": [200, 400, 18]}, + "system_size" : {"value": [0.1, 0.2, 0.009], "unit": "m"}, + "density" : {"value": 8000, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 191e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 42408, "unit": "J/m^2"}, + "horizon" : {"value": 0.002, "unit": "m"}, + "impactor_velocity" : {"value": 16, "unit": "m/sec"}, + "final_time" : {"value": 70e-6, "unit": "s"}, + "timestep" : {"value": 0.133e-6, "unit": "s"}, + "output_frequency" : {"value": 500}, + "output_reference" : {"value": true} +} diff --git a/examples/inputs/kalthoff_winkler_m6.json b/examples/inputs/kalthoff_winkler_m6.json new file mode 100644 index 00000000..182c797c --- /dev/null +++ b/examples/inputs/kalthoff_winkler_m6.json @@ -0,0 +1,13 @@ +{ + "num_cells" : {"value": [300, 600, 27]}, + "system_size" : {"value": [0.1, 0.2, 0.009], "unit": "m"}, + "density" : {"value": 8000, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 191e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 42408, "unit": "J/m^2"}, + "horizon" : {"value": 0.002, "unit": "m"}, + "impactor_velocity" : {"value": 16, "unit": "m/sec"}, + "final_time" : {"value": 70e-6, "unit": "s"}, + "timestep" : {"value": 0.133e-6, "unit": "s"}, + "output_frequency" : {"value": 500}, + "output_reference" : {"value": true} +} diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index 1e0c7f8b..efaf58b3 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -76,34 +76,11 @@ int main( int argc, char* argv[] ) auto u = particles->sliceDisplacement(); auto v = particles->sliceVelocity(); auto rho = particles->sliceDensity(); - auto temp = particles->sliceTemperature(); - - /* - auto init_functor = KOKKOS_LAMBDA( const int pid ) - { - double a = 0.001; - double r0 = 0.25; - double l = 0.07; - double norm = std::sqrt( x( pid, 0 ) * x( pid, 0 ) + - x( pid, 1 ) * x( pid, 1 ) + - x( pid, 2 ) * x( pid, 2 ) ); - double diff = norm - r0; - double arg = diff * diff / l / l; - for ( int d = 0; d < 3; d++ ) - { - double comp = 0.0; - if ( norm > 0.0 ) - comp = x( pid, d ) / norm; - u( pid, d ) = a * std::exp( -arg ) * comp; - v( pid, d ) = 0.0; - } - rho( pid ) = rho0; - }; - - */ + // auto temp = particles->sliceTemperature(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) { - temp( pid ) = 5000 * x( pid, 1 ); + // temp( pid ) = 5000 * x( pid, 1 ); rho( pid ) = rho0; }; particles->updateParticles( exec_space{}, init_functor ); @@ -122,39 +99,6 @@ int main( int argc, char* argv[] ) int mpi_rank; MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); Kokkos::View count( "c", 1 ); - // double dx = particles->dx[0]; - - /* - auto measure_profile = KOKKOS_LAMBDA( const int pid ) - { - if ( x( pid, 1 ) < dx / 2.0 && x( pid, 1 ) > -dx / 2.0 && - x( pid, 2 ) < dx / 2.0 && x( pid, 2 ) > -dx / 2.0 ) - { - auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); - profile( c, 0 ) = x( pid, 0 ); - profile( c, 1 ) = u( pid, 0 ); - } - }; - */ - - /* - Kokkos::RangePolicy policy( 0, x.size() ); - Kokkos::parallel_for( "displacement_profile", policy, measure_profile - ); auto count_host = Kokkos::create_mirror_view_and_copy( - Kokkos::HostSpace{}, count ); auto profile_host = - Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, profile - ); - */ - /* - std::fstream fout; - std::string file_name = "displacement_profile.txt"; - fout.open( file_name, std::ios::app ); - for ( int p = 0; p < count_host( 0 ); p++ ) - { - fout << mpi_rank << " " << profile_host( p, 0 ) << " " - << profile_host( p, 1 ) << std::endl; - } - */ } MPI_Finalize(); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 9a7fc4bf..e0b45a7b 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -203,7 +203,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles ) + void apply( ExecSpace, ParticleType& particles, double t ) { auto f = particles.sliceForce(); auto index_space = _index_space._view; diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index c8246090..ae305db4 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -1037,7 +1037,7 @@ class Force> auto c = _model.c; auto alpha = _model.alpha; // auto break_coeff = _model.bond_break_coeff; - auto break_coeff = _model.s0; + auto s0 = _model.s0; const auto vol = particles.sliceVolume(); const auto temp = particles.sliceTemperature(); const auto nofail = particles.sliceNoFail(); From d9adf22d286935f8c66a37f5da5f81e659a9cf19 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 12 Dec 2023 14:30:32 -0700 Subject: [PATCH 06/15] fixup: hardcoded temp BC --- src/CabanaPD_Boundary.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index e0b45a7b..87f88bc0 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -238,7 +238,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles ) + void apply( ExecSpace, ParticleType& particles, double ) { auto f = particles.sliceForce(); auto index_space = _index_space._view; From 67410b48f131d1f85177d49bd9f0ae3c9bf5725a Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Tue, 12 Dec 2023 18:37:56 -0500 Subject: [PATCH 07/15] Changes for time-dependent thermal deformation --- examples/thermal_deformation.cpp | 16 +++++++++++++--- src/CabanaPD_Boundary.hpp | 7 +++++-- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index efaf58b3..b008ad48 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -74,13 +74,25 @@ int main( int argc, char* argv[] ) // Define particle initialization. auto x = particles->sliceReferencePosition(); auto u = particles->sliceDisplacement(); + auto f = particles->sliceForce(); auto v = particles->sliceVelocity(); auto rho = particles->sliceDensity(); // auto temp = particles->sliceTemperature(); + // Domain to apply b.c. + CabanaPD::RegionBoundary domain1( low_corner[0], high_corner[0], + low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + std::vector domain = { domain1 }; + + double b0 = 0.0; + auto bc = + createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, + exec_space{}, *particles, domain, b0 ); + auto init_functor = KOKKOS_LAMBDA( const int pid ) { - // temp( pid ) = 5000 * x( pid, 1 ); + // temp( pid ) = 5000 * x( pid, 1 ) * t_final; rho( pid ) = rho0; }; particles->updateParticles( exec_space{}, init_functor ); @@ -90,8 +102,6 @@ int main( int argc, char* argv[] ) cabana_pd->init_force(); cabana_pd->run(); - x = particles->sliceReferencePosition(); - u = particles->sliceDisplacement(); double num_cell_x = inputs.num_cells[0]; auto profile = Kokkos::View( Kokkos::ViewAllocateWithoutInitializing( "displacement_profile" ), diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 87f88bc0..35bbc2e9 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -164,7 +164,10 @@ struct BoundaryCondition Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); + // This is specifically for the thermal deformation problem temp( pid ) = 5000 * x( pid, 1 ) * t; + std::cout << "1: temp(" << pid << ")=" << temp( pid ) + << std::endl; } ); } }; @@ -205,7 +208,7 @@ struct BoundaryCondition template void apply( ExecSpace, ParticleType& particles, double t ) { - auto f = particles.sliceForce(); + auto temp = particles.sliceTemperature(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; @@ -238,7 +241,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles, double ) + void apply( ExecSpace, ParticleType& particles, double t ) { auto f = particles.sliceForce(); auto index_space = _index_space._view; From cc106dd697f1f07947b4e08f761b3ae507c5897c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 18 Dec 2023 10:48:30 -0500 Subject: [PATCH 08/15] fix thermal problem --- examples/thermal_deformation.cpp | 2 +- src/CabanaPD_Boundary.hpp | 5 ++--- src/CabanaPD_Particles.hpp | 14 +++++++------- src/CabanaPD_Solver.hpp | 5 +++-- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index b008ad48..91b57163 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -98,7 +98,7 @@ int main( int argc, char* argv[] ) particles->updateParticles( exec_space{}, init_functor ); auto cabana_pd = CabanaPD::createSolverElastic( - inputs, particles, force_model ); + inputs, particles, force_model, bc ); cabana_pd->init_force(); cabana_pd->run(); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 35bbc2e9..7d75d88e 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -166,8 +166,6 @@ struct BoundaryCondition auto pid = index_space( b ); // This is specifically for the thermal deformation problem temp( pid ) = 5000 * x( pid, 1 ) * t; - std::cout << "1: temp(" << pid << ")=" << temp( pid ) - << std::endl; } ); } }; @@ -261,7 +259,8 @@ template auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space, Particles particles, std::vector planes, - const double initial_guess = 0.5 ) + const double value, + const double initial_guess = 1.1 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 23114a78..9571f8c9 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -442,15 +442,15 @@ class Particles Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, n_local, getPosition( use_reference ), sliceStrainEnergy(), - sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() ); + sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), + sliceTemperature() ); #else #ifdef Cabana_ENABLE_SILO - Cabana::Grid::Experimental::SiloParticleOutput:: - writePartialRangeTimeStep( - "particles", local_grid->globalGrid(), output_step, output_time, - 0, n_local, getPosition( use_reference ), sliceStrainEnergy(), - sliceForce(), sliceDisplacement(), sliceVelocity(), - sliceDamage() ); + Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep( + "particles", local_grid->globalGrid(), output_step, output_time, 0, + n_local, getPosition( use_reference ), sliceStrainEnergy(), + sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), + sliceTemperature() ); #else log( std::cout, "No particle output enabled." ); #endif diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index c71c7ef1..cc44e5f5 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -210,7 +210,7 @@ class SolverElastic computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, 0.0 ); particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); @@ -249,7 +249,8 @@ class SolverElastic force_time += force_timer.seconds(); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, + step * inputs->timestep ); // Integrate - velocity Verlet second half. integrate_timer.reset(); From 8bb6693d30e627650535f17aac190d448114a576 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 13 Feb 2024 11:36:36 -0500 Subject: [PATCH 09/15] fixup: rebase mistakes --- src/CabanaPD_Boundary.hpp | 41 ++++++++++++++++++++++++++++++++++++--- src/CabanaPD_Solver.hpp | 6 ++---- 2 files changed, 40 insertions(+), 7 deletions(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 7d75d88e..9c96dc70 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -131,6 +131,10 @@ struct ForceUpdateBCTag { }; +struct TempBCTag +{ +}; + struct ZeroBCTag { }; @@ -154,6 +158,37 @@ struct BoundaryCondition _index_space.update( exec_space, particles, plane ); } + template + void apply( ExecSpace, ParticleType& particles, double ) + { + auto user = _user_functor; + auto index_space = _index_space._view; + Kokkos::RangePolicy policy( 0, index_space.size() ); + Kokkos::parallel_for( + "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { + auto pid = index_space( b ); + user( pid ); + } ); + } +}; + +template +struct BoundaryCondition +{ + BCIndexSpace _index_space; + + BoundaryCondition( BCIndexSpace bc_index_space, UserFunctor user ) + : _index_space( bc_index_space ) + { + } + + template + void update( ExecSpace exec_space, Particles particles, + RegionBoundary plane ) + { + _index_space.update( exec_space, particles, plane ); + } + template void apply( ExecSpace, ParticleType& particles, double t ) { @@ -179,7 +214,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType ) + void apply( ExecSpace, ParticleType, double ) { } }; @@ -204,7 +239,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles, double t ) + void apply( ExecSpace, ParticleType& particles, double ) { auto temp = particles.sliceTemperature(); auto index_space = _index_space._view; @@ -239,7 +274,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles, double t ) + void apply( ExecSpace, ParticleType& particles, double ) { auto f = particles.sliceForce(); auto index_space = _index_space._view; diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index cc44e5f5..dd0f5153 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -249,8 +249,7 @@ class SolverElastic force_time += force_timer.seconds(); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles, - step * inputs->timestep ); + boundary_condition.apply( exec_space(), *particles, step * dt ); // Integrate - velocity Verlet second half. integrate_timer.reset(); @@ -494,8 +493,7 @@ class SolverFracture force_time += force_timer.seconds(); // Add boundary condition. - boundary_condition.apply( exec_space{}, *particles, - step * inputs->timestep ); + boundary_condition.apply( exec_space{}, *particles, step * dt ); // Integrate - velocity Verlet second half. integrate_timer.reset(); From 7e46e167f89a6904dd47ef507b05c92057da2d4c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 16 Jan 2024 21:15:35 -0500 Subject: [PATCH 10/15] Convert temperature example to json --- examples/inputs/thermal.json | 12 +++++++++ examples/thermal_deformation.cpp | 46 ++++++++++++++++---------------- 2 files changed, 35 insertions(+), 23 deletions(-) create mode 100644 examples/inputs/thermal.json diff --git a/examples/inputs/thermal.json b/examples/inputs/thermal.json new file mode 100644 index 00000000..93cf1d42 --- /dev/null +++ b/examples/inputs/thermal.json @@ -0,0 +1,12 @@ +{ + "num_cells" : {"value": [100, 30, 3]}, + "system_size" : {"value": [1.0, 0.3, 0.03], "unit": ""}, + "density" : {"value": 3980.0, "unit": ""}, + "elastic_modulus" : {"value": 370e+9, "unit": ""}, + "shear_modulus" : {"value": 0.5, "unit": ""}, + "horizon" : {"value": 0.0375, "unit": ""}, + "final_time" : {"value": 0.0093, "unit": ""}, + "timestep" : {"value": 7.5E-7, "unit": ""}, + "output_frequency": {"value": 200}, + "output_reference": {"value": true} +} diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index 91b57163..e363bf36 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -29,26 +29,22 @@ int main( int argc, char* argv[] ) using exec_space = Kokkos::DefaultExecutionSpace; using memory_space = typename exec_space::memory_space; - std::array num_cell = { 100, 30, 3 }; - std::array low_corner = { -0.5, 0.0, -0.015 }; - std::array high_corner = { 0.5, 0.3, 0.015 }; - double t_final = 0.0093; - double dt = 7.5E-7; - int output_frequency = 200; - bool output_reference = true; - - double rho0 = 3980; // [kg/m^3] - double E = 370e+9; // [Pa] + CabanaPD::Inputs inputs( argv[1] ); + double E = inputs["elastic_modulus"]; + double rho0 = inputs["density"]; double nu = 0.25; // unitless double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] - double alpha = 7.5e-6; // [1/oC] - double delta = 0.03; + double delta = inputs["horizon"]; + double alpha = 7.5e-6; // [1/oC] // Reference temperature // double temp0 = 0.0; + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array num_cells = inputs["num_cells"]; int m = std::floor( - delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); + delta / ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) ); int halo_width = m + 1; // Just to be safe. // Choose force model type. @@ -60,16 +56,11 @@ int main( int argc, char* argv[] ) // CabanaPD::ForceModel; // model_type force_model( delta, K, G ); - CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final, - dt, output_frequency, output_reference ); - inputs.read_args( argc, argv ); - // Create particles from mesh. // Does not set displacements, velocities, etc. auto particles = std::make_shared< CabanaPD::Particles>( - exec_space(), inputs.low_corner, inputs.high_corner, - inputs.num_cells, halo_width ); + exec_space(), low_corner, high_corner, num_cells, halo_width ); // Define particle initialization. auto x = particles->sliceReferencePosition(); @@ -85,10 +76,19 @@ int main( int argc, char* argv[] ) low_corner[2], high_corner[2] ); std::vector domain = { domain1 }; - double b0 = 0.0; + // This really shouldn't be applied as a BC. + auto bc_op = KOKKOS_LAMBDA( const int pid ) + { + // Get a modifiable copy of temperature. + auto p_t = particles_t.getParticleView( pid ); + // Get a copy of the position. + auto p_x = particles_x.getParticle( pid ); + auto yref = + Cabana::get( p_x, CabanaPD::Field::ReferencePosition(), 1 ); + temp( pid ) += 5000 * yref * t; + }; auto bc = - createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, - exec_space{}, *particles, domain, b0 ); + createBoundaryCondition( bc_op, exec_space{}, *particles, domain ); auto init_functor = KOKKOS_LAMBDA( const int pid ) { @@ -102,7 +102,7 @@ int main( int argc, char* argv[] ) cabana_pd->init_force(); cabana_pd->run(); - double num_cell_x = inputs.num_cells[0]; + double num_cell_x = num_cells[0]; auto profile = Kokkos::View( Kokkos::ViewAllocateWithoutInitializing( "displacement_profile" ), num_cell_x ); From adee4593f456462430ff54b812881a2dfd8ce280 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 13 Feb 2024 11:45:41 -0500 Subject: [PATCH 11/15] fixup: from rebase and temperature json --- examples/crack_branching.cpp | 2 +- examples/kalthoff_winkler.cpp | 6 +++--- examples/thermal_deformation.cpp | 15 ++------------- src/CabanaPD_Boundary.hpp | 18 ++++++++++-------- 4 files changed, 16 insertions(+), 25 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index cc20169d..9f7b7ebb 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -68,7 +68,7 @@ int main( int argc, char* argv[] ) // model_type force_model( delta, K, G0 ); using model_type = CabanaPD::ForceModel; - model_type force_model( delta, K, G0 ); + model_type force_model( delta, K, 0.0, G0 ); // Create particles from mesh. // Does not set displacements, velocities, etc. diff --git a/examples/kalthoff_winkler.cpp b/examples/kalthoff_winkler.cpp index a783a91a..8e87f4ed 100644 --- a/examples/kalthoff_winkler.cpp +++ b/examples/kalthoff_winkler.cpp @@ -36,7 +36,7 @@ int main( int argc, char* argv[] ) double K = E / ( 3.0 * ( 1.0 - 2.0 * nu ) ); double rho0 = inputs["density"]; double G0 = inputs["fracture_energy"]; - // double G = E / ( 2.0 * ( 1.0 + nu ) ); // Only for LPS. + double G = E / ( 2.0 * ( 1.0 + nu ) ); // Only for LPS. // PD horizon double delta = inputs["horizon"]; @@ -73,8 +73,8 @@ int main( int argc, char* argv[] ) // Choose force model type. using model_type = - CabanaPD::ForceModel; - model_type force_model( delta, K, G, G0 ); + CabanaPD::ForceModel; + model_type force_model( delta, K, G, G0 ); // Create particles from mesh. // Does not set displacements, velocities, etc. diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index e363bf36..4486f92d 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -76,19 +76,8 @@ int main( int argc, char* argv[] ) low_corner[2], high_corner[2] ); std::vector domain = { domain1 }; - // This really shouldn't be applied as a BC. - auto bc_op = KOKKOS_LAMBDA( const int pid ) - { - // Get a modifiable copy of temperature. - auto p_t = particles_t.getParticleView( pid ); - // Get a copy of the position. - auto p_x = particles_x.getParticle( pid ); - auto yref = - Cabana::get( p_x, CabanaPD::Field::ReferencePosition(), 1 ); - temp( pid ) += 5000 * yref * t; - }; - auto bc = - createBoundaryCondition( bc_op, exec_space{}, *particles, domain ); + auto bc = createBoundaryCondition( CabanaPD::TempBCTag{}, 5000.0, + exec_space{}, *particles, domain ); auto init_functor = KOKKOS_LAMBDA( const int pid ) { diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 9c96dc70..a0d33046 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -159,7 +159,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles, double ) + void apply( ExecSpace, ParticleType&, double ) { auto user = _user_functor; auto index_space = _index_space._view; @@ -172,13 +172,15 @@ struct BoundaryCondition } }; -template -struct BoundaryCondition +template +struct BoundaryCondition { + double _value; BCIndexSpace _index_space; - BoundaryCondition( BCIndexSpace bc_index_space, UserFunctor user ) - : _index_space( bc_index_space ) + BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + : _value( value ) + , _index_space( bc_index_space ) { } @@ -194,13 +196,14 @@ struct BoundaryCondition { auto temp = particles.sliceTemperature(); auto x = particles.sliceReferencePosition(); + auto value = _value; auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); // This is specifically for the thermal deformation problem - temp( pid ) = 5000 * x( pid, 1 ) * t; + temp( pid ) = value * x( pid, 1 ) * t; } ); } }; @@ -241,7 +244,7 @@ struct BoundaryCondition template void apply( ExecSpace, ParticleType& particles, double ) { - auto temp = particles.sliceTemperature(); + auto f = particles.sliceForce(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; @@ -294,7 +297,6 @@ template auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space, Particles particles, std::vector planes, - const double value, const double initial_guess = 1.1 ) { using memory_space = typename Particles::memory_space; From a3c51343878bab442521275b853545343657b538 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 18 Jan 2024 11:55:00 -0500 Subject: [PATCH 12/15] Custom geometry from Pablo --- examples/inputs/thermal.json | 14 ++++++++------ examples/thermal_deformation.cpp | 17 +++++++++++++++-- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/examples/inputs/thermal.json b/examples/inputs/thermal.json index 93cf1d42..2db5784c 100644 --- a/examples/inputs/thermal.json +++ b/examples/inputs/thermal.json @@ -1,12 +1,14 @@ { - "num_cells" : {"value": [100, 30, 3]}, - "system_size" : {"value": [1.0, 0.3, 0.03], "unit": ""}, - "density" : {"value": 3980.0, "unit": ""}, - "elastic_modulus" : {"value": 370e+9, "unit": ""}, + "num_cells" : {"value": [100, 100, 43]}, + "system_size" : {"value": [0.028, 0.028, 0.012], "unit": "m"}, + "radius" : {"value": 0.0085, "unit": "m"}, + "density" : {"value": 19.25e+3, "unit": ""}, + "elastic_modulus" : {"value": 340e+9, "unit": ""}, "shear_modulus" : {"value": 0.5, "unit": ""}, - "horizon" : {"value": 0.0375, "unit": ""}, + "thermal_coeff" : {"value": 4.2E-6, "unit": ""}, + "horizon" : {"value": 8.4E-4, "unit": ""}, "final_time" : {"value": 0.0093, "unit": ""}, "timestep" : {"value": 7.5E-7, "unit": ""}, - "output_frequency": {"value": 200}, + "output_frequency": {"value": 10}, "output_reference": {"value": true} } diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index 4486f92d..04b6c866 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -36,7 +36,7 @@ int main( int argc, char* argv[] ) double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] double delta = inputs["horizon"]; - double alpha = 7.5e-6; // [1/oC] + double alpha = inputs["thermal_coeff"]; // [1/oC] // Reference temperature // double temp0 = 0.0; @@ -60,7 +60,20 @@ int main( int argc, char* argv[] ) // Does not set displacements, velocities, etc. auto particles = std::make_shared< CabanaPD::Particles>( - exec_space(), low_corner, high_corner, num_cells, halo_width ); + low_corner, high_corner, num_cells, halo_width ); + // Do not create particles in the center. + double x_center = 0.0; + double y_center = -0.0005; + double radius = inputs["radius"]; + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + if ( ( ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ) ) < + radius * radius ) + return false; + return true; + }; + particles->createParticles( exec_space(), init_op ); // Define particle initialization. auto x = particles->sliceReferencePosition(); From 390dcd0a7d137862c4b62d3ae99f60d4e8cc8a89 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 13 Feb 2024 11:50:19 -0500 Subject: [PATCH 13/15] fixup: warning --- src/CabanaPD_Particles.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 9571f8c9..341dffd0 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -446,11 +446,12 @@ class Particles sliceTemperature() ); #else #ifdef Cabana_ENABLE_SILO - Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep( - "particles", local_grid->globalGrid(), output_step, output_time, 0, - n_local, getPosition( use_reference ), sliceStrainEnergy(), - sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), - sliceTemperature() ); + Cabana::Grid::Experimental::SiloParticleOutput:: + writePartialRangeTimeStep( + "particles", local_grid->globalGrid(), output_step, output_time, + 0, n_local, getPosition( use_reference ), sliceStrainEnergy(), + sliceForce(), sliceDisplacement(), sliceVelocity(), + sliceDamage(), sliceTemperature() ); #else log( std::cout, "No particle output enabled." ); #endif From 9b68c264531747229e234168affd56b498b1930e Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 13 Feb 2024 11:55:50 -0500 Subject: [PATCH 14/15] fixup: constructor change --- examples/thermal_deformation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index 04b6c866..68405d61 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -60,7 +60,7 @@ int main( int argc, char* argv[] ) // Does not set displacements, velocities, etc. auto particles = std::make_shared< CabanaPD::Particles>( - low_corner, high_corner, num_cells, halo_width ); + exec_space(), low_corner, high_corner, num_cells, halo_width ); // Do not create particles in the center. double x_center = 0.0; double y_center = -0.0005; From 100050aaad79c5c1f963c53450f564b79412d5f9 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 14 Feb 2024 15:09:12 -0500 Subject: [PATCH 15/15] fixup: set corner for temperature BC --- examples/thermal_deformation.cpp | 1 + src/CabanaPD_Boundary.hpp | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/thermal_deformation.cpp b/examples/thermal_deformation.cpp index 68405d61..3673d8a8 100644 --- a/examples/thermal_deformation.cpp +++ b/examples/thermal_deformation.cpp @@ -91,6 +91,7 @@ int main( int argc, char* argv[] ) auto bc = createBoundaryCondition( CabanaPD::TempBCTag{}, 5000.0, exec_space{}, *particles, domain ); + bc.setCorner( low_corner[1] ); auto init_functor = KOKKOS_LAMBDA( const int pid ) { diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index a0d33046..8902acaa 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -177,6 +177,7 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; + double _low_corner = 0.0; BoundaryCondition( const double value, BCIndexSpace bc_index_space ) : _value( value ) @@ -184,6 +185,8 @@ struct BoundaryCondition { } + void setCorner( const double low ) { _low_corner = low; } + template void update( ExecSpace exec_space, Particles particles, RegionBoundary plane ) @@ -197,13 +200,14 @@ struct BoundaryCondition auto temp = particles.sliceTemperature(); auto x = particles.sliceReferencePosition(); auto value = _value; + auto low_corner = _low_corner; auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); // This is specifically for the thermal deformation problem - temp( pid ) = value * x( pid, 1 ) * t; + temp( pid ) = value * ( x( pid, 1 ) - low_corner ) * t; } ); } };