Merge pull request #44 from streeve/input_file

Input file

CMakeLists.txt

@@ -28,6 +28,17 @@ endmacro()
 CabanaPD_check_optional( OPTION HDF5 )
 CabanaPD_check_optional( OPTION SILO )
 
+if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.24")
+ cmake_policy(SET CMP0135 NEW)
+endif()
+find_package(nlohmann_json 3.10.0 QUIET)
+if(NOT NLOHMANN_JSON_FOUND)
+ include(FetchContent)
+ # Using most recent release here
+ FetchContent_Declare(json URL https://github.com/nlohmann/json/releases/download/v3.11.2/json.tar.xz)
+ FetchContent_MakeAvailable(json)
+endif()
+
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cmake/CabanaPD_Config.cmakein
  ${CMAKE_CURRENT_BINARY_DIR}/CabanaPD_Config.cmake @ONLY)
 

examples/crack_branching.cpp

@@ -29,47 +29,34 @@ int main( int argc, char* argv[] )
  using exec_space = Kokkos::DefaultExecutionSpace;
  using memory_space = typename exec_space::memory_space;
 
- // Plate dimensions (m)
- double height = 0.1;
- double width = 0.04;
- double thickness = 0.002;
-
- // Domain
- std::array<int, 3> num_cell = { 400, 160, 8 }; // 400 x 160 x 8
- double low_x = -0.5 * height;
- double low_y = -0.5 * width;
- double low_z = -0.5 * thickness;
- double high_x = 0.5 * height;
- double high_y = 0.5 * width;
- double high_z = 0.5 * thickness;
- std::array<double, 3> low_corner = { low_x, low_y, low_z };
- std::array<double, 3> high_corner = { high_x, high_y, high_z };
-
- // Time
- double t_final = 43e-6;
- double dt = 5e-8;
- double output_frequency = 5;
- bool output_reference = true;
+ CabanaPD::Inputs inputs( argv[1] );
 
  // Material constants
- double E = 72e+9; // [Pa]
- double nu = 0.25; // unitless
- double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa]
- double rho0 = 2440; // [kg/m^3]
- double G0 = 3.8; // [J/m^2]
+ double E = inputs["elastic_modulus"];
+ double nu = 0.25;
+ double K = E / ( 3 * ( 1 - 2 * nu ) );
+ double rho0 = inputs["density"];
+ double G0 = inputs["fracture_energy"];
 
  // PD horizon
- double delta = 0.001 + 1e-10;
+ double delta = inputs["horizon"];
+ delta += 1e-10;
 
  // FIXME: set halo width based on delta
+ std::array<double, 3> low_corner = inputs["low_corner"];
+ std::array<double, 3> high_corner = inputs["high_corner"];
+ std::array<int, 3> num_cells = inputs["num_cells"];
  int m = std::floor(
- delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) );
- int halo_width = m + 1;
+ delta / ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) );
+ int halo_width = m + 1; // Just to be safe.
 
  // Prenotch
+ double height = inputs["system_size"][0];
+ double thickness = inputs["system_size"][2];
  double L_prenotch = height / 2.0;
  double y_prenotch1 = 0.0;
- Kokkos::Array<double, 3> p01 = { low_x, y_prenotch1, low_z };
+ Kokkos::Array<double, 3> p01 = { low_corner[0], y_prenotch1,
+ low_corner[2] };
  Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
  Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
  Kokkos::Array<Kokkos::Array<double, 3>, 1> notch_positions = { p01 };
@@ -79,16 +66,12 @@ int main( int argc, char* argv[] )
  using model_type =
  CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
  model_type force_model( delta, K, G0 );
- 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<memory_space, typename model_type::base_model>>(
- 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();
@@ -97,14 +80,16 @@ int main( int argc, char* argv[] )
  auto rho = particles->sliceDensity();
  auto nofail = particles->sliceNoFail();
 
- // Relying on uniform grid here.
  double dy = particles->dx[1];
- double b0 = 2e6 / dy; // Pa/m
-
- CabanaPD::RegionBoundary plane1( low_x, high_x, low_y - dy, low_y + dy,
- low_z, high_z );
- CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy,
- high_y + dy, low_z, high_z );
+ double sigma0 = inputs["traction"];
+ double b0 = sigma0 / dy;
+
+ CabanaPD::RegionBoundary plane1( low_corner[0], high_corner[0],
+ low_corner[1] - dy, low_corner[1] + dy,
+ low_corner[2], high_corner[2] );
+ CabanaPD::RegionBoundary plane2(
+ low_corner[0], high_corner[0], high_corner[1] - dy,
+ high_corner[1] + dy, low_corner[2], high_corner[2] );
  std::vector<CabanaPD::RegionBoundary> planes = { plane1, plane2 };
  auto bc =
  createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{},
@@ -120,7 +105,6 @@ int main( int argc, char* argv[] )
  };
  particles->updateParticles( exec_space{}, init_functor );
 
- // FIXME: use createSolver to switch backend at runtime.
  auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
  inputs, particles, force_model, bc, prenotch );
  cabana_pd->init_force();

examples/elastic_wave.cpp

@@ -25,22 +25,26 @@ int main( int argc, char* 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<int, 3> num_cell = { 41, 41, 41 };
- std::array<double, 3> low_corner = { -0.5, -0.5, -0.5 };
- std::array<double, 3> high_corner = { 0.5, 0.5, 0.5 };
- double t_final = 0.6;
- double dt = 0.01;
- int output_frequency = 5;
- bool output_reference = true;
- double K = 1.0;
- double G = 0.5;
- double delta = 0.075;
+ CabanaPD::Inputs inputs( argv[1] );
+
+ // Material constants
+ auto K = inputs["bulk_modulus"];
+ double G = inputs["shear_modulus"];
+ double rho0 = inputs["density"];
+
+ // PD horizon
+ double delta = inputs["horizon"];
+ delta += 1e-10;
+
+ // FIXME: set halo width based on delta
+ std::array<double, 3> low_corner = inputs["low_corner"];
+ std::array<double, 3> high_corner = inputs["high_corner"];
+ std::array<int, 3> 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.
@@ -51,16 +55,11 @@ int main( int argc, char* argv[] )
  CabanaPD::ForceModel<CabanaPD::LinearLPS, CabanaPD::Elastic>;
  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<memory_space, typename model_type::base_model>>(
- 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();
@@ -86,7 +85,7 @@ int main( int argc, char* argv[] )
  u( pid, d ) = a * std::exp( -arg ) * comp;
  v( pid, d ) = 0.0;
  }
- rho( pid ) = 100.0;
+ rho( pid ) = rho0;
  };
  particles->updateParticles( exec_space{}, init_functor );
 
@@ -97,7 +96,7 @@ int main( int argc, char* argv[] )
 
  x = particles->sliceReferencePosition();
  u = particles->sliceDisplacement();
- double num_cell_x = inputs.num_cells[0];
+ int num_cell_x = num_cells[0];
  auto profile = Kokkos::View<double* [2], memory_space>(
  Kokkos::ViewAllocateWithoutInitializing( "displacement_profile" ),
  num_cell_x );

examples/inputs/crack_branching.json

@@ -0,0 +1,13 @@
+{
+ "num_cells" : {"value": [400, 160, 8]},
+ "system_size" : {"value": [0.1, 0.04, 0.002], "unit": "m"},
+ "density" : {"value": 2440, "unit": "kg/m^3"},
+ "elastic_modulus" : {"value": 72e+9, "unit": "Pa"},
+ "fracture_energy" : {"value": 3.8, "unit": "J/m^2"},
+ "horizon" : {"value": 0.001, "unit": "m"}, 
+ "traction" : {"value": 2e6, "unit": "Pa"},
+ "final_time" : {"value": 43e-6, "unit": "s"},
+ "timestep" : {"value": 5e-8, "unit": "s"},
+ "output_frequency": {"value": 5},
+ "output_reference": {"value": true}
+}

examples/inputs/elastic_wave.json

@@ -0,0 +1,12 @@
+{
+ "num_cells" : {"value": [41, 41, 41]},
+ "system_size" : {"value": [1.0, 1.0, 1.0], "unit": ""},
+ "density" : {"value": 100.0, "unit": ""},
+ "bulk_modulus" : {"value": 1.0, "unit": ""},
+ "shear_modulus" : {"value": 0.5, "unit": ""},
+ "horizon" : {"value": 0.075, "unit": ""}, 
+ "final_time" : {"value": 0.6, "unit": ""},
+ "timestep" : {"value": 0.01, "unit": ""},
+ "output_frequency": {"value": 5},
+ "output_reference": {"value": true}
+}

examples/inputs/kalthoff_winkler.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": 10},
+ "output_reference" : {"value": true}
+}

examples/kalthoff_winkler.cpp

@@ -25,73 +25,65 @@ int main( int argc, char* 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;
 
- // Plate dimension)
- double height = 0.1; // [m] (100 mm)
- double width = 0.2; // [m] (200 mm)
- double thickness = 0.009; // [m] ( 9 mm)
-
- // Domain
- // This is a relatively large example for CPU - reduce the number of
- // cells and increase delta if needed. Note this is also a relatively
- // small example for GPU.
- std::array<int, 3> num_cell = { 151, 301, 14 };
- std::array<double, 3> low_corner = { -0.5 * height, -0.5 * width,
- -0.5 * thickness };
- std::array<double, 3> high_corner = { 0.5 * height, 0.5 * width,
- 0.5 * thickness };
- double t_final = 70e-6;
- double dt = 0.133e-6;
- int output_frequency = 10;
- bool output_reference = true;
+ CabanaPD::Inputs inputs( argv[1] );
 
  // Material constants
- double E = 191e+9; // [Pa]
- double nu = 1.0 / 3.0; // unitless
- double K = E / ( 3.0 * ( 1.0 - 2.0 * nu ) ); // [Pa]
- double rho0 = 8000; // [kg/m^3]
- double G0 = 42408; // [J/m^2]
+ double E = inputs["elastic_modulus"];
+ double nu = 1.0 / 3.0;
+ 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 v0 = 16; // [m/sec] (Half impactor's velocity)
- double L_prenotch = 0.05; // [m] (50 mm)
- double y_prenotch1 = -0.025; // [m] (-25 mm)
- double y_prenotch2 = 0.025; // [m] ( 25 mm)
- Kokkos::Array<double, 3> p01 = { low_corner[0], y_prenotch1,
- low_corner[2] };
- Kokkos::Array<double, 3> p02 = { low_corner[0], y_prenotch2,
- low_corner[2] };
+ // PD horizon
+ double delta = inputs["horizon"];
+ delta += 1e-10;
+
+ // Impactor velocity
+ double v0 = inputs["impactor_velocity"];
+
+ // FIXME: set halo width based on delta
+ std::array<double, 3> low_corner = inputs["low_corner"];
+ std::array<double, 3> high_corner = inputs["high_corner"];
+ std::array<int, 3> num_cells = inputs["num_cells"];
+ int m = std::floor(
+ delta / ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) );
+ int halo_width = m + 1; // Just to be safe.
+
+ // Prenotches
+ std::array<double, 3> system_size = inputs["system_size"];
+ double height = system_size[0];
+ double width = system_size[1];
+ double thickness = system_size[2];
+ double L_prenotch = height / 2.0;
+ double y_prenotch1 = -width / 8.0;
+ double y_prenotch2 = width / 8.0;
+ double low_x = low_corner[0];
+ double low_z = low_corner[2];
+ Kokkos::Array<double, 3> p01 = { low_x, y_prenotch1, low_z };
+ Kokkos::Array<double, 3> p02 = { low_x, y_prenotch2, low_z };
  Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
  Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
  Kokkos::Array<Kokkos::Array<double, 3>, 2> notch_positions = { p01,
  p02 };
  CabanaPD::Prenotch<2> prenotch( v1, v2, notch_positions );
 
- double delta = 0.0020000001;
- 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<CabanaPD::PMB, CabanaPD::Fracture>;
  model_type force_model( delta, K, G0 );
  // using model_type =
  // CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
  // model_type force_model( delta, K, G, G0 );
- 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<memory_space, typename model_type::base_model>>(
- 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();
@@ -100,7 +92,6 @@ int main( int argc, char* argv[] )
  auto rho = particles->sliceDensity();
 
  double dx = particles->dx[0];
-
  double x_bc = -0.5 * height;
  CabanaPD::RegionBoundary plane(
  x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25,

src/CMakeLists.txt

@@ -11,6 +11,6 @@ target_include_directories(CabanaPD INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}
  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
  $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>)
 
-target_link_libraries(CabanaPD INTERFACE Cabana::cabanacore Cabana::Cajita)
+target_link_libraries(CabanaPD INTERFACE Cabana::cabanacore Cabana::Cajita nlohmann_json::nlohmann_json)
 
 install(TARGETS CabanaPD DESTINATION ${CMAKE_INSTALL_LIBDIR})