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})