diff --git a/single-node-refactor/CMakeLists.txt b/single-node-refactor/CMakeLists.txt index 78d9f53e7..c8f008e1f 100755 --- a/single-node-refactor/CMakeLists.txt +++ b/single-node-refactor/CMakeLists.txt @@ -54,9 +54,6 @@ include_directories(src/input) add_subdirectory(src/input) -include_directories(src/Solvers/SGH_solver) -add_subdirectory(src/Solvers/SGH_solver) - add_subdirectory(src) diff --git a/single-node-refactor/input-impact-erosion.yaml b/single-node-refactor/input-impact-erosion.yaml index ebb56847c..ff05f5b20 100644 --- a/single-node-refactor/input-impact-erosion.yaml +++ b/single-node-refactor/input-impact-erosion.yaml @@ -85,20 +85,13 @@ materials: - 1.0E-4 - 1.0 erosion_model: basic - void_mat_id: 2 erode_tension_val: 2.0e-7 erode_density_val: 0.01 - - # eroded material - - material: - id: 2 - eos_model_type: decoupled - eos_model: void # background air - material: - id: 3 + id: 2 eos_model_type: decoupled eos_model: gamma_law_gas # strength_model: none @@ -115,7 +108,7 @@ regions: # air - fill_volume: type: global - material_id: 3 + material_id: 2 den: 0.010 sie: 3.0e-4 velocity: cartesian diff --git a/single-node-refactor/input-origin-erosion.yaml b/single-node-refactor/input-origin-erosion.yaml index 3805d6dba..24f3d3051 100644 --- a/single-node-refactor/input-origin-erosion.yaml +++ b/single-node-refactor/input-origin-erosion.yaml @@ -5,7 +5,7 @@ dynamic_options: dt_min: 1.e-8 dt_max: 1.e-2 dt_start: 1.e-5 - cycle_stop: 300000 + cycle_stop: 30000 # mesh_options: @@ -79,14 +79,23 @@ materials: - 1.0E-14 - 1.0 erosion_model: basic - void_mat_id: 1 erode_tension_val: 0.00003 erode_density_val: 0.02 + - material: id: 1 eos_model_type: decoupled - eos_model: void + eos_model: gamma_law_gas + # strength_model: none + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 1.666666666666667 + - 1.0E-14 + - 1.0 regions: - fill_volume: diff --git a/single-node-refactor/input.yaml b/single-node-refactor/input.yaml index 4802a3c53..a657ecbf2 100644 --- a/single-node-refactor/input.yaml +++ b/single-node-refactor/input.yaml @@ -42,9 +42,6 @@ boundary_conditions: geometry: x_plane direction: x_dir value: 0.0 - u: 0.0 - v: 0.0 - w: 1.0 type: reflected_velocity @@ -63,6 +60,8 @@ boundary_conditions: direction: z_dir value: 0.0 type: reflected_velocity + + materials: - material: diff --git a/single-node-refactor/src/CMakeLists.txt b/single-node-refactor/src/CMakeLists.txt index 859a059e9..71b573cf8 100755 --- a/single-node-refactor/src/CMakeLists.txt +++ b/single-node-refactor/src/CMakeLists.txt @@ -75,5 +75,9 @@ endif() include_directories(common) -add_executable(Fierro main.cpp solver.cpp ) -target_link_libraries(Fierro PRIVATE matar parse_yaml sgh_solver Kokkos::kokkos) +# Add SGH Solver +include_directories(Solvers/SGH_solver/include) +add_subdirectory(Solvers/SGH_solver) + +add_executable(Fierro main.cpp solver.cpp ${YAML_SRC_Files} ${SGH_SRC_Files} ) +target_link_libraries(Fierro PRIVATE matar Kokkos::kokkos) #sgh_solver parse_yaml diff --git a/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt b/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt index 76ac6da4e..9dc2c4025 100755 --- a/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt +++ b/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt @@ -1,6 +1,5 @@ cmake_minimum_required(VERSION 3.1.3) - find_package(Matar REQUIRED) find_package(Kokkos REQUIRED) @@ -19,29 +18,18 @@ endif() include_directories(include) include_directories(src) - -set(SRC_Files -src/boundary.cpp -src/energy_sgh.cpp -src/force_sgh.cpp -src/position.cpp -src/momentum.cpp -src/properties.cpp -src/sgh_solve.cpp -src/time_integration.cpp) - -# set(SGH_Solver_SRC src/sgh_solver.cpp ) - -message("\n ****** ADDING SGH LIBRARY ******** \n ") - -add_library(sgh_solver ${SRC_Files} ) # ${SGH_Solver_SRC} -target_include_directories(sgh_solver PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include) - - -target_link_libraries(sgh_solver matar Kokkos::kokkos) - - - - - - +message("\n ****** ADDING SGH SOURCE FILES ******** \n ") + +set(SGH_SRC_Files +${CMAKE_CURRENT_SOURCE_DIR}/src/boundary.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/energy_sgh.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/force_sgh.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/position.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/momentum.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/properties.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/sgh_solve.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/time_integration.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/sgh_setup.cpp +${CMAKE_CURRENT_SOURCE_DIR}/include/sgh_solver.h +PARENT_SCOPE +) \ No newline at end of file diff --git a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h index 4c93c4ff7..a8dec6470 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h +++ b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h @@ -35,13 +35,16 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #ifndef SGH_SOLVER_H #define SGH_SOLVER_H +#include "solver.h" #include "geometry_new.h" #include "matar.h" #include "simulation_parameters.h" #include "boundary_conditions.h" #include "material.h" -#include "solver.h" - +#include "mesh.h" +#include "state.h" +#include "io_utils.h" +#include "dynamic_options.h" using namespace mtr; // matar namespace @@ -70,27 +73,25 @@ class SGH : public Solver // Initialize data specific to the SGH solver void initialize(SimulationParameters_t& SimulationParamaters, Material_t& Materials, - BoundaryCondition_t& Boundary) const override + mesh_t& mesh, + BoundaryCondition_t& Boundary, + State_t& State) const override { + // stuff goes here } ///////////////////////////////////////////////////////////////////////////// /// /// \fn setup /// - /// \brief Calls setup_sgh, which initializes state, and material data + /// \brief Calls setup_sgh, which initializes state and material data /// ///////////////////////////////////////////////////////////////////////////// void setup(SimulationParameters_t& SimulationParamaters, Material_t& Materials, - BoundaryCondition_t& Boundary, mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) const override - { - } + BoundaryCondition_t& Boundary, + State_t& State) override; ///////////////////////////////////////////////////////////////////////////// /// @@ -104,10 +105,7 @@ class SGH : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) override; + State_t& State) override; ///////////////////////////////////////////////////////////////////////////// /// @@ -132,6 +130,28 @@ class SGH : public Solver // Any finalize goes here, remove allocated memory, etc } + // **** Functions defined in sgh_setup.cpp **** // + void fill_regions_sgh( + const Material_t& Materials, + const mesh_t& mesh, + const DCArrayKokkos & node_coords, + DCArrayKokkos & node_vel, + DCArrayKokkos & GaussPoint_den, + DCArrayKokkos & GaussPoint_sie, + DCArrayKokkos & elem_mat_id, + DCArrayKokkos & voxel_elem_mat_id, + const CArrayKokkos & region_fills, + const CArray & region_fills_host, + const size_t num_fills, + const size_t num_elems, + const size_t num_nodes, + const size_t rk_num_bins) const; + + void init_corner_node_masses_zero( + const mesh_t& mesh, + const DCArrayKokkos& node_mass, + const DCArrayKokkos& corner_mass) const; + // **** Functions defined in boundary.cpp **** // void boundary_velocity( const mesh_t& mesh, @@ -147,19 +167,26 @@ class SGH : public Solver // **** Functions defined in energy_sgh.cpp **** // void update_energy( - double rk_alpha, - double dt, + const double rk_alpha, + const double dt, const mesh_t& mesh, const DCArrayKokkos& node_vel, const DCArrayKokkos& node_coords, - DCArrayKokkos& MaterialPoints_sie, + const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& corner_force) const; + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t corners_in_mat_elem, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems) const; // **** Functions defined in force_sgh.cpp **** // void get_force( const Material_t& Materials, const mesh_t& mesh, + const DCArrayKokkos& GaussPoints_vol, + const DCArrayKokkos& GaussPoints_div, + const DCArrayKokkos& GaussPoints_eroded, + const DCArrayKokkos& corner_force, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, const DCArrayKokkos& MaterialPoints_den, @@ -167,20 +194,23 @@ class SGH : public Solver const DCArrayKokkos& MaterialPoints_pres, const DCArrayKokkos& MaterialPoints_stress, const DCArrayKokkos& MaterialPoints_sspd, - const DCArrayKokkos& GaussPoints_vol, - const DCArrayKokkos& GaussPoints_div, - const DCArrayKokkos& GaussPoints_mat_id, - const DCArrayKokkos& MaterialPoints_eroded, - const DCArrayKokkos& corner_force, + const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems, + const size_t mat_id, const double fuzz, const double small, - const DCArrayKokkos& MaterialPoints_statev, const double dt, const double rk_alpha) const; void get_force_2D( const Material_t& Materials, const mesh_t& mesh, + const DCArrayKokkos& GaussPoints_vol, + const DCArrayKokkos& GaussPoints_div, + const DCArrayKokkos& corner_force, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, const DCArrayKokkos& MaterialPoints_den, @@ -188,13 +218,14 @@ class SGH : public Solver const DCArrayKokkos& MaterialPoints_pres, const DCArrayKokkos& MaterialPoints_stress, const DCArrayKokkos& MaterialPoints_sspd, - const DCArrayKokkos& GaussPoints_vol, - const DCArrayKokkos& GaussPoints_div, - const DCArrayKokkos& GaussPoints_mat_id, - DCArrayKokkos& corner_force, + const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems, + const size_t mat_id, const double fuzz, const double small, - const DCArrayKokkos& MaterialPoints_statev, const double dt, const double rk_alpha) const; @@ -261,35 +292,40 @@ class SGH : public Solver const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_den, - DCArrayKokkos& MaterialPoints_pres, - DCArrayKokkos& MaterialPoints_stress, - DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& GaussPoints_vol, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& GaussPoints_mat_id, const DCArrayKokkos& MaterialPoints_statev, - const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& GaussPoints_eroded, + const DCArrayKokkos& MaterialToMeshMaps_elem, const double dt, - const double rk_alpha) const; + const double rk_alpha, + const size_t num_material_elems, + const size_t mat_id) const; void update_state2D( const Material_t& Materials, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_den, - DCArrayKokkos& MaterialPoints_pres, - DCArrayKokkos& MaterialPoints_stress, - DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& GaussPoints_vol, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& GaussPoints_mat_id, const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& GaussPoints_eroded, + const DCArrayKokkos& MaterialToMeshMaps_elem, const double dt, - const double rk_alpha) const; + const double rk_alpha, + const size_t num_material_elems, + const size_t mat_id) const; // **** Functions defined in time_integration.cpp **** // // NOTE: Consider pulling up @@ -300,14 +336,18 @@ class SGH : public Solver DCArrayKokkos& MaterialPoints_stress, const size_t num_dims, const size_t num_elems, - const size_t num_nodes) const; + const size_t num_nodes, + const size_t num_mat_points) const; void get_timestep( mesh_t& mesh, DCArrayKokkos& node_coords, DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_sspd, DCArrayKokkos& GaussPoints_vol, + DCArrayKokkos& MaterialPoints_sspd, + DCArrayKokkos& MaterialPoints_eroded, + DCArrayKokkos& MaterialToMeshMaps_elem, + size_t num_mat_elems, double time_value, const double graphics_time, const double time_final, @@ -321,8 +361,11 @@ class SGH : public Solver mesh_t& mesh, DCArrayKokkos& node_coords, DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_sspd, DCArrayKokkos& GaussPoints_vol, + DCArrayKokkos& MaterialPoints_sspd, + DCArrayKokkos& MaterialPoints_eroded, + DCArrayKokkos& MaterialToMeshMaps_elem, + size_t num_mat_elems, double time_value, const double graphics_time, const double time_final, @@ -364,4 +407,35 @@ class SGH : public Solver const double rk_alpha); }; +void calc_extensive_node_mass(const CArrayKokkos& node_extensive_mass, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass, + const double num_dims, + const double num_nodes); + +void calc_node_areal_mass(const mesh_t& mesh, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass, + CArrayKokkos node_extensive_mass, + double tiny); + +double sum_domain_internal_energy(const DCArrayKokkos& MaterialPoints_mass, + const DCArrayKokkos& MaterialPoints_sie, + const size_t num_mat_points); + +double sum_domain_kinetic_energy(const mesh_t& mesh, + const DCArrayKokkos& node_vel, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass); + +double sum_domain_material_mass(const DCArrayKokkos& MaterialPoints_mass, + const size_t num_mat_points); + +double sum_domain_node_mass(const mesh_t& mesh, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass); + +void set_corner_force_zero(const mesh_t& mesh, + const DCArrayKokkos& corner_force); + #endif // end HEADER_H diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp index bb28f30e6..0063e366f 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/energy_sgh.cpp @@ -50,30 +50,50 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. /// \param A view into the corner force data /// ///////////////////////////////////////////////////////////////////////////// -void SGH::update_energy(double rk_alpha, - double dt, +void SGH::update_energy(const double rk_alpha, + const double dt, const mesh_t& mesh, const DCArrayKokkos& node_vel, const DCArrayKokkos& node_coords, - DCArrayKokkos& MaterialPoints_sie, + const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& corner_force) const + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t corners_in_mat_elem, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems + ) const { + // loop over all the elements in the mesh - FOR_ALL(elem_gid, 0, mesh.num_elems, { + FOR_ALL(mat_elem_lid, 0, num_mat_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + // the material point index = the material elem index for a 1-point element + size_t mat_point_lid = mat_elem_lid; + + double MaterialPoints_power = 0.0; + // --- tally the contribution from each corner to the element --- // Loop over the nodes in the element for (size_t node_lid = 0; node_lid < mesh.num_nodes_in_elem; node_lid++) { + + // corner lid and node lid size_t corner_lid = node_lid; // Get node global id for the local node id size_t node_gid = mesh.nodes_in_elem(elem_gid, node_lid); // Get the corner global id for the local corner id - size_t corner_gid = mesh.corners_in_elem(elem_gid, corner_lid); + //size_t corner_gid = mesh.corners_in_elem(elem_gid, corner_lid); + + // Get the material corner lid + size_t mat_corner_lid = corners_in_mat_elem(mat_elem_lid, corner_lid); + double node_radius = 1; if (mesh.num_dims == 2) { @@ -82,14 +102,18 @@ void SGH::update_energy(double rk_alpha, // calculate the Power=F dot V for this corner for (size_t dim = 0; dim < mesh.num_dims; dim++) { + double half_vel = (node_vel(1, node_gid, dim) + node_vel(0, node_gid, dim)) * 0.5; - MaterialPoints_power += corner_force(corner_gid, dim) * node_radius * half_vel; + MaterialPoints_power += MaterialCorners_force(mat_corner_lid, dim) * node_radius * half_vel; + } // end for dim } // end for node_lid // update the specific energy - MaterialPoints_sie(1, elem_gid) = MaterialPoints_sie(0, elem_gid) - - rk_alpha * dt / MaterialPoints_mass(elem_gid) * MaterialPoints_power; + MaterialPoints_sie(1, mat_point_lid) = MaterialPoints_sie(0, mat_point_lid) - + rk_alpha * dt / MaterialPoints_mass(mat_point_lid) * MaterialPoints_power; + + }); // end parallel loop over the elements return; diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp index 6faed136e..dac1858f0 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp @@ -60,29 +60,43 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. /// ///////////////////////////////////////////////////////////////////////////// void SGH::get_force(const Material_t& Materials, - const mesh_t& mesh, - const DCArrayKokkos& node_coords, - const DCArrayKokkos& node_vel, - const DCArrayKokkos& MaterialPoints_den, - const DCArrayKokkos& MaterialPoints_sie, - const DCArrayKokkos& MaterialPoints_pres, - const DCArrayKokkos& MaterialPoints_stress, - const DCArrayKokkos& MaterialPoints_sspd, - const DCArrayKokkos& GaussPoints_vol, - const DCArrayKokkos& GaussPoints_div, - const DCArrayKokkos& GaussPoints_mat_id, - const DCArrayKokkos& GaussPoints_eroded, - const DCArrayKokkos& corner_force, - const double fuzz, - const double small, - const DCArrayKokkos& MaterialPoints_statev, - const double dt, - const double rk_alpha) const + const mesh_t& mesh, + const DCArrayKokkos& GaussPoints_vol, + const DCArrayKokkos& GaussPoints_div, + const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& corner_force, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_vel, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_sie, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t corners_in_mat_elem, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems, + const size_t mat_id, + const double fuzz, + const double small, + const double dt, + const double rk_alpha) const { + const size_t num_dims = 3; + const size_t num_nodes_in_elem = 8; + + + // --- calculate the forces acting on the nodes from the element --- - FOR_ALL(elem_gid, 0, mesh.num_elems, { - const size_t num_dims = 3; - const size_t num_nodes_in_elem = 8; + FOR_ALL(mat_elem_lid, 0, num_mat_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + // the material point index = the material elem index for a 1-point element + size_t mat_point_lid = mat_elem_lid; + // total Cauchy stress double tau_array[9]; @@ -118,11 +132,9 @@ void SGH::get_force(const Material_t& Materials, // element volume double vol = GaussPoints_vol(elem_gid); - // the id for the material in this element - size_t mat_id = GaussPoints_mat_id(elem_gid); // create a view of the stress_matrix - ViewCArrayKokkos stress(&MaterialPoints_stress(1, elem_gid, 0, 0), 3, 3); + ViewCArrayKokkos stress(&MaterialPoints_stress(1, mat_point_lid, 0, 0), 3, 3); // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 8); @@ -172,7 +184,7 @@ void SGH::get_force(const Material_t& Materials, // add the pressure if a decoupled model is used if (Materials.MaterialEnums(mat_id).EOSType == model::decoupledEOSType) { for (int i = 0; i < num_dims; i++) { - tau(i, i) -= MaterialPoints_pres(elem_gid); + tau(i, i) -= MaterialPoints_pres(mat_point_lid); } // end for } @@ -231,8 +243,8 @@ void SGH::get_force(const Material_t& Materials, // if there is no velocity change, then use the surface area // normal as the shock direction mag = sqrt(area_normal(node_lid, 0) * area_normal(node_lid, 0) - + area_normal(node_lid, 1) * area_normal(node_lid, 1) - + area_normal(node_lid, 2) * area_normal(node_lid, 2) ); + + area_normal(node_lid, 1) * area_normal(node_lid, 1) + + area_normal(node_lid, 2) * area_normal(node_lid, 2) ); // estimate of the shock direction for (int dim = 0; dim < num_dims; dim++) { @@ -242,13 +254,13 @@ void SGH::get_force(const Material_t& Materials, // cell divergence indicates compression or expansions if (div < 0) { // element in compression - muc(node_lid) = MaterialPoints_den(elem_gid) * - (Materials.MaterialFunctions(mat_id).q1 * MaterialPoints_sspd(elem_gid) + + muc(node_lid) = MaterialPoints_den(mat_point_lid) * + (Materials.MaterialFunctions(mat_id).q1 * MaterialPoints_sspd(mat_point_lid) + Materials.MaterialFunctions(mat_id).q2 * mag_vel); } else{ // element in expansion - muc(node_lid) = MaterialPoints_den(elem_gid) * - (Materials.MaterialFunctions(mat_id).q1ex * MaterialPoints_sspd(elem_gid) + + muc(node_lid) = MaterialPoints_den(mat_point_lid) * + (Materials.MaterialFunctions(mat_id).q1ex * MaterialPoints_sspd(mat_point_lid) + Materials.MaterialFunctions(mat_id).q2ex * mag_vel); } // end if on divergence sign @@ -262,15 +274,15 @@ void SGH::get_force(const Material_t& Materials, // direction mu_term = muc(node_lid) * fabs(shock_dir(0) * area_normal(node_lid, 0) - + shock_dir(1) * area_normal(node_lid, 1) - + shock_dir(2) * area_normal(node_lid, 2) ); + + shock_dir(1) * area_normal(node_lid, 1) + + shock_dir(2) * area_normal(node_lid, 2) ); } else{ // Using a full tensoral Riemann jump relation mu_term = muc(node_lid) * sqrt(area_normal(node_lid, 0) * area_normal(node_lid, 0) - + area_normal(node_lid, 1) * area_normal(node_lid, 1) - + area_normal(node_lid, 2) * area_normal(node_lid, 2) ); + + area_normal(node_lid, 1) * area_normal(node_lid, 1) + + area_normal(node_lid, 2) * area_normal(node_lid, 2) ); } sum(0) += mu_term * vel(0); @@ -336,7 +348,7 @@ void SGH::get_force(const Material_t& Materials, double omega = 20.0; // 20.0; // weighting factor on Mach number double third = 1.0 / 3.0; double c_length = pow(vol, third); // characteristic length - double alpha = fmin(1.0, omega * (c_length * fabs(div)) / (MaterialPoints_sspd(elem_gid) + fuzz) ); + double alpha = fmin(1.0, omega * (c_length * fabs(div)) / (MaterialPoints_sspd(mat_point_lid) + fuzz) ); // use Mach based detector with standard shock detector @@ -353,6 +365,8 @@ void SGH::get_force(const Material_t& Materials, // loop over the each node in the elem for (size_t node_lid = 0; node_lid < num_nodes_in_elem; node_lid++) { + + // the local corner id is the local node id size_t corner_lid = node_lid; // Get corner gid @@ -361,22 +375,32 @@ void SGH::get_force(const Material_t& Materials, // Get node gid size_t node_gid = mesh.nodes_in_elem(elem_gid, node_lid); - // loop over dimension + // Get the material corner lid + size_t mat_corner_lid = corners_in_mat_elem(mat_elem_lid, corner_lid); + //printf("corner difference = %zu \n", mat_corner_lid-corner_gid); - if (GaussPoints_eroded(elem_gid) == true) { // material(mat_id).blank_mat_id) + // loop over dimensions and calc corner forces + if (MaterialPoints_eroded(mat_point_lid) == true) { // material(mat_id).blank_mat_id) for (int dim = 0; dim < num_dims; dim++) { corner_force(corner_gid, dim) = 0.0; + MaterialCorners_force(mat_corner_lid, dim) = 0.0; } } + else{ for (int dim = 0; dim < num_dims; dim++) { - corner_force(corner_gid, dim) = + + double force_component = area_normal(node_lid, 0) * tau(0, dim) + area_normal(node_lid, 1) * tau(1, dim) + area_normal(node_lid, 2) * tau(2, dim) + phi * muc(node_lid) * (vel_star(dim) - node_vel(1, node_gid, dim)); + + MaterialCorners_force(mat_corner_lid, dim) = force_component; + corner_force(corner_gid, dim) += force_component; // tally all forces to the corner } // end loop over dimension - } + } // end if + } // end for loop over nodes in elem // --- Update Stress --- @@ -389,12 +413,12 @@ void SGH::get_force(const Material_t& Materials, // --- call strength model --- Materials.MaterialFunctions(mat_id).calc_stress(MaterialPoints_pres, MaterialPoints_stress, - elem_gid, + mat_point_lid, mat_id, MaterialPoints_statev, MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(elem_gid), + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1,mat_point_lid), vel_grad, elem_node_gids, node_coords, @@ -434,28 +458,49 @@ void SGH::get_force(const Material_t& Materials, /// ///////////////////////////////////////////////////////////////////////////// void SGH::get_force_2D(const Material_t& Materials, - const mesh_t& mesh, - const DCArrayKokkos& node_coords, - const DCArrayKokkos& node_vel, - const DCArrayKokkos& MaterialPoints_den, - const DCArrayKokkos& MaterialPoints_sie, - const DCArrayKokkos& MaterialPoints_pres, - const DCArrayKokkos& MaterialPoints_stress, - const DCArrayKokkos& MaterialPoints_sspd, - const DCArrayKokkos& GaussPoints_vol, - const DCArrayKokkos& GaussPoints_div, - const DCArrayKokkos& GaussPoints_mat_id, - DCArrayKokkos& corner_force, - const double fuzz, - const double small, - const DCArrayKokkos& MaterialPoints_statev, - const double dt, - const double rk_alpha) const + const mesh_t& mesh, + const DCArrayKokkos& GaussPoints_vol, + const DCArrayKokkos& GaussPoints_div, + const DCArrayKokkos& corner_force, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_vel, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_sie, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& MaterialCorners_force, + const corners_in_mat_t corners_in_mat_elem, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems, + const size_t mat_id, + const double fuzz, + const double small, + const double dt, + const double rk_alpha) const { + const size_t num_dims = 2; + const size_t num_nodes_in_elem = 4; + + // set corner force to zero + FOR_ALL(corner_gid, 0, mesh.num_corners, { + for (int dim = 0; dim < num_dims; dim++) { + corner_force(corner_gid, dim) = 0.0; + } + }); // end parallel for corners + + // --- calculate the forces acting on the nodes from the element --- - FOR_ALL_CLASS(elem_gid, 0, mesh.num_elems, { - const size_t num_dims = 2; - const size_t num_nodes_in_elem = 4; + FOR_ALL(mat_elem_lid, 0, num_mat_elems, { + + // get mesh elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + //size_t guass_gid = elem_gid; // 1 gauss point per element + + // the material point index = the material elem index for a 1-point element + size_t mat_point_lid = mat_elem_lid; // total Cauchy stress double tau_array[9]; @@ -489,7 +534,7 @@ void SGH::get_force_2D(const Material_t& Materials, ViewCArrayKokkos vel_grad(vel_grad_array, 3, 3); // create a view of the stress_matrix - ViewCArrayKokkos stress(&MaterialPoints_stress(1, elem_gid, 0, 0), 3, 3); + ViewCArrayKokkos stress(&MaterialPoints_stress(1, mat_point_lid, 0, 0), 3, 3); // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 4); @@ -546,7 +591,7 @@ void SGH::get_force_2D(const Material_t& Materials, // add the pressure for (int i = 0; i < 3; i++) { - tau(i, i) -= MaterialPoints_pres(elem_gid); + tau(i, i) -= MaterialPoints_pres(mat_point_lid); } // end for // ---- Multidirectional Approximate Riemann solver (MARS) ---- @@ -611,15 +656,14 @@ void SGH::get_force_2D(const Material_t& Materials, } // end if mag_vel // cell divergence indicates compression or expansions - size_t mat_id = GaussPoints_mat_id(elem_gid); if (div < 0) { // element in compression - muc(node_lid) = MaterialPoints_den(elem_gid) * - (Materials.MaterialFunctions(mat_id).q1 * MaterialPoints_sspd(elem_gid) + + muc(node_lid) = MaterialPoints_den(mat_point_lid) * + (Materials.MaterialFunctions(mat_id).q1 * MaterialPoints_sspd(mat_point_lid) + Materials.MaterialFunctions(mat_id).q2 * mag_vel); } else{ // element in expansion - muc(node_lid) = MaterialPoints_den(elem_gid) * - (Materials.MaterialFunctions(mat_id).q1ex * MaterialPoints_sspd(elem_gid) + + muc(node_lid) = MaterialPoints_den(mat_point_lid) * + (Materials.MaterialFunctions(mat_id).q1ex * MaterialPoints_sspd(mat_point_lid) + Materials.MaterialFunctions(mat_id).q2ex * mag_vel); } // end if on divergence sign @@ -754,7 +798,6 @@ void SGH::get_force_2D(const Material_t& Materials, // --- Update Stress --- // calculate the new stress at the next rk level, if it is a increment_based model - size_t mat_id = GaussPoints_mat_id(elem_gid); // increment_based elastic plastic model if (Materials.MaterialEnums(mat_id).StrengthType == model::incrementBased) { @@ -769,7 +812,7 @@ void SGH::get_force_2D(const Material_t& Materials, // MaterialPoints_statev, // MaterialPoints_sspd, // MaterialPoints_den(elem_gid), - // MaterialPoints_sie(elem_gid), + // MaterialPoints_sie(1,elem_gid), // vel_grad, // elem_node_gids, // node_coords, diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp index 1a748ab08..1724ec12c 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/properties.cpp @@ -57,44 +57,133 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. /// \param The current Runge Kutta integration alpha value /// ///////////////////////////////////////////////////////////////////////////// -void SGH::update_state(const Material_t& Materials, +void SGH::update_state( + const Material_t& Materials, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_den, - DCArrayKokkos& MaterialPoints_pres, - DCArrayKokkos& MaterialPoints_stress, - DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& GaussPoints_vol, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& GaussPoints_mat_id, const DCArrayKokkos& MaterialPoints_statev, - const DCArrayKokkos& GaussPoints_eroded, + const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& MaterialToMeshMaps_elem, const double dt, - const double rk_alpha) const + const double rk_alpha, + const size_t num_material_elems, + const size_t mat_id) const { - // loop over all the elements in the mesh - FOR_ALL(elem_gid, 0, mesh.num_elems, { - const size_t num_dims = mesh.num_dims; - const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + + const size_t num_dims = mesh.num_dims; + + + // --- pressure --- + if (Materials.MaterialEnums.host(mat_id).EOSType == model::decoupledEOSType) { + + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; - // cut out the node_gids for this element - ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); + // --- Density --- + MaterialPoints_den(mat_point_lid) = MaterialPoints_mass(mat_point_lid) / GaussPoints_vol(gauss_gid); + + // --- Pressure --- + Materials.MaterialFunctions(mat_id).calc_pressure( + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1, mat_point_lid), + Materials.eos_global_vars); + + // --- Sound Speed --- + Materials.MaterialFunctions(mat_id).calc_sound_speed( + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1, mat_point_lid), + Materials.eos_global_vars); + + }); // end parallel for over mat elem lid + + } // if decoupled EOS + else { + // only calculate density as pressure and sound speed come from the coupled strength model // --- Density --- - MaterialPoints_den(elem_gid) = MaterialPoints_mass(elem_gid) / GaussPoints_vol(elem_gid); + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; + + + // --- Density --- + MaterialPoints_den(mat_point_lid) = MaterialPoints_mass(mat_point_lid) / GaussPoints_vol(gauss_gid); + + }); // end parallel for over mat elem lid + Kokkos::fence(); + } // end if + + + // --- Stress --- + + // state_based elastic plastic model + if (Materials.MaterialEnums.host(mat_id).StrengthType == model::stateBased) { + + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; - size_t mat_id = GaussPoints_mat_id(elem_gid); - // --- Stress --- - // state_based elastic plastic model - if (Materials.MaterialEnums(mat_id).StrengthType == model::stateBased) { // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); + // --- Density --- - MaterialPoints_den(elem_gid) = MaterialPoints_mass(elem_gid) / GaussPoints_vol(elem_gid); + MaterialPoints_den(mat_point_lid) = MaterialPoints_mass(mat_point_lid) / GaussPoints_vol(gauss_gid); + // corner area normals double area_array[24]; @@ -115,81 +204,80 @@ void SGH::update_state(const Material_t& Materials, GaussPoints_vol(elem_gid), elem_gid); + // --- call strength model --- Materials.MaterialFunctions(mat_id).calc_stress( - MaterialPoints_pres, - MaterialPoints_stress, - elem_gid, - mat_id, - MaterialPoints_statev, - MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(elem_gid), - vel_grad, - elem_node_gids, - node_coords, - node_vel, - GaussPoints_vol(elem_gid), - dt, - rk_alpha); - } // end logical on state_based strength model - - // apply the element erosion model - //if (material(mat_id).erosion_type == model::erosion) { - // // starting simple, but in the future call an erosion model - // if (MaterialPoints_pres(elem_gid) <= material(mat_id).erode_tension_val - // || MaterialPoints_den(elem_gid) <= material(mat_id).erode_density_val) { - // GaussPoints_mat_id(elem_gid) = material(mat_id).void_mat_id; - // - // GaussPoints_eroded(elem_gid) = true; - // } // end if - //} // end if - if (Materials.MaterialFunctions(mat_id).erode != NULL) { - + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1,mat_point_lid), + vel_grad, + elem_node_gids, + node_coords, + node_vel, + GaussPoints_vol(gauss_gid), + dt, + rk_alpha); + + + }); // end parallel for over mat elem lid + + } // end if state_based strength model + + + // --- mat point erosion --- + if (Materials.MaterialEnums.host(mat_id).ErosionModels != model::noErosion) { + + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; // --- Element erosion model --- Materials.MaterialFunctions(mat_id).erode( - MaterialPoints_pres, + MaterialPoints_eroded, MaterialPoints_stress, - GaussPoints_eroded, - GaussPoints_mat_id, - elem_gid, - Materials.MaterialFunctions(mat_id).void_mat_id, + MaterialPoints_pres(mat_point_lid), + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1, mat_point_lid), + MaterialPoints_sspd(mat_point_lid), Materials.MaterialFunctions(mat_id).erode_tension_val, Materials.MaterialFunctions(mat_id).erode_density_val, - MaterialPoints_sspd, - MaterialPoints_den, - MaterialPoints_sie(1, elem_gid)); + mat_point_lid); - } // end if - if (Materials.MaterialEnums(mat_id).EOSType == model::decoupledEOSType) { - // --- Pressure --- - Materials.MaterialFunctions(mat_id).calc_pressure( - MaterialPoints_pres, - MaterialPoints_stress, - elem_gid, - GaussPoints_mat_id(elem_gid), - MaterialPoints_statev, - MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(1, elem_gid), - Materials.eos_global_vars); - // --- Sound Speed --- - Materials.MaterialFunctions(mat_id).calc_sound_speed( - MaterialPoints_pres, - MaterialPoints_stress, - elem_gid, - GaussPoints_mat_id(elem_gid), - MaterialPoints_statev, - MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(1, elem_gid), - Materials.eos_global_vars); - } - }); // end parallel for - Kokkos::fence(); + // apply a void eos if mat_point is eroded + if(MaterialPoints_eroded(mat_point_lid)){ + + MaterialPoints_pres(mat_point_lid) = 0.0; + MaterialPoints_sspd(mat_point_lid) = 1.0e-32; + + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < 3; j++) { + MaterialPoints_stress(1, mat_point_lid, i, j) = 0.0; + } + } // end for i,j + + } // end if on eroded + + + + }); // end parallel for + } // end if elem errosion return; } // end method to update state @@ -217,127 +305,213 @@ void SGH::update_state(const Material_t& Materials, /// \param The current Runge Kutta integration alpha value /// ///////////////////////////////////////////////////////////////////////////// -void SGH::update_state2D(const Material_t& Materials, +void SGH::update_state2D( + const Material_t& Materials, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_den, - DCArrayKokkos& MaterialPoints_pres, - DCArrayKokkos& MaterialPoints_stress, - DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, const DCArrayKokkos& MaterialPoints_sie, const DCArrayKokkos& GaussPoints_vol, const DCArrayKokkos& MaterialPoints_mass, - const DCArrayKokkos& GaussPoints_mat_id, const DCArrayKokkos& MaterialPoints_statev, + const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& MaterialToMeshMaps_elem, const double dt, - const double rk_alpha) const + const double rk_alpha, + const size_t num_material_elems, + const size_t mat_id) const { - // loop over all the elements in the mesh - FOR_ALL(elem_gid, 0, mesh.num_elems, { - const size_t num_dims = mesh.num_dims; - const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + + const size_t num_dims = mesh.num_dims; + + // --- Density --- + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; - // cut out the node_gids for this element - ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); // --- Density --- - MaterialPoints_den(elem_gid) = MaterialPoints_mass(elem_gid) / GaussPoints_vol(elem_gid); + MaterialPoints_den(mat_point_lid) = MaterialPoints_mass(mat_point_lid) / GaussPoints_vol(gauss_gid); - size_t mat_id = GaussPoints_mat_id(elem_gid); + }); // end parallel for over mat elem lid + Kokkos::fence(); - // --- Stress --- - // state_based elastic plastic model - if (Materials.MaterialEnums(mat_id).StrengthType == model::stateBased) { - // cut out the node_gids for this element - ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); - // --- Density --- - MaterialPoints_den(elem_gid) = MaterialPoints_mass(elem_gid) / GaussPoints_vol(elem_gid); + // --- pressure --- + if (Materials.MaterialEnums(mat_id).EOSType == model::decoupledEOSType) { - // corner area normals - double area_array[8]; - ViewCArrayKokkos area(area_array, num_nodes_in_elem, num_dims); + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { - // velocity gradient - double vel_grad_array[4]; - ViewCArrayKokkos vel_grad(vel_grad_array, num_dims, num_dims); + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); - // get the B matrix which are the OUTWARD corner area normals - geometry::get_bmatrix(area, elem_gid, node_coords, elem_node_gids); - // --- Calculate the velocity gradient --- - get_velgrad(vel_grad, - elem_node_gids, - node_vel, - area, - GaussPoints_vol(elem_gid), - elem_gid); + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; - // --- call strength model --- - // Material.MaterialFunctions(mat_id).strength_model(MaterialPoints_pres, - // MaterialPoints_stress, - // elem_gid, - // mat_id, - // MaterialPoints_statev, - // MaterialPoints_sspd, - // MaterialPoints_den(elem_gid), - // MaterialPoints_sie(elem_gid), - // vel_grad, - // elem_node_gids, - // node_coords, - // node_vel, - // GaussPoints_vol(elem_gid), - // dt, - // rk_alpha); - } // end logical on state_based strength model - - // --- Erosion --- - // apply the element erosion model - //if (Materials.MaterialFunctions(mat_id).erode != NULL) { - // - // // --- Element erosion model --- - // material.MaterialFunctions(mat_id).erode(MaterialPoints_pres, - // MaterialPoints_stress, - // GaussPoints_eroded, - // GaussPoints_mat_id, - // elem_gid, - // material(mat_id).void_mat_id, - // material(mat_id).erode_tension_val, - // material(mat_id).erode_density_val, - // MaterialPoints_sspd, - // MaterialPoints_den, - // MaterialPoints_sie(1, elem_gid)); - //} // end if - - // --- Pressure --- - if (Materials.MaterialEnums(mat_id).EOSType == model::decoupledEOSType) { + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; // --- Pressure --- Materials.MaterialFunctions(mat_id).calc_pressure( - MaterialPoints_pres, - MaterialPoints_stress, - elem_gid, - GaussPoints_mat_id(elem_gid), - MaterialPoints_statev, - MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(1, elem_gid), - Materials.eos_global_vars); + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(0, mat_point_lid), + Materials.eos_global_vars); // --- Sound Speed --- Materials.MaterialFunctions(mat_id).calc_sound_speed( - MaterialPoints_pres, - MaterialPoints_stress, - elem_gid, - GaussPoints_mat_id(elem_gid), - MaterialPoints_statev, - MaterialPoints_sspd, - MaterialPoints_den(elem_gid), - MaterialPoints_sie(1, elem_gid), - Materials.eos_global_vars); - } - }); // end parallel for - Kokkos::fence(); + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(0, mat_point_lid), + Materials.eos_global_vars); + + }); // end parallel for over mat elem lid + + } // if decoupled EOS + + + // --- Stress --- + + // state_based elastic plastic model + if (Materials.MaterialEnums.host(mat_id).StrengthType == model::stateBased) { + + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; + + + // cut out the node_gids for this element + ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); + + + // --- Density --- + MaterialPoints_den(mat_point_lid) = MaterialPoints_mass(mat_point_lid) / GaussPoints_vol(gauss_gid); + + + // corner area normals + double area_array[24]; + ViewCArrayKokkos area(area_array, num_nodes_in_elem, num_dims); + + // velocity gradient + double vel_grad_array[9]; + ViewCArrayKokkos vel_grad(vel_grad_array, num_dims, num_dims); + + // get the B matrix which are the OUTWARD corner area normals + geometry::get_bmatrix(area, elem_gid, node_coords, elem_node_gids); + + // --- Calculate the velocity gradient --- + get_velgrad(vel_grad, + elem_node_gids, + node_vel, + area, + GaussPoints_vol(elem_gid), + elem_gid); + + // --- call strength model --- + Materials.MaterialFunctions(mat_id).calc_stress( + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1,mat_point_lid), + vel_grad, + elem_node_gids, + node_coords, + node_vel, + GaussPoints_vol(gauss_gid), + dt, + rk_alpha); + + + }); // end parallel for over mat elem lid + + } // end if state_based strength model + + + // --- mat point erosion --- + if (Materials.MaterialEnums.host(mat_id).ErosionModels != model::noErosion) { + + // loop over all the elements the material lives in + FOR_ALL(mat_elem_lid, 0, num_material_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + + // get the material points for this material + // Note, with the SGH method, they are equal + size_t mat_point_lid = mat_elem_lid; + + // for this method, gauss point is equal to elem_gid + size_t gauss_gid = elem_gid; + + // --- Element erosion model --- + Materials.MaterialFunctions(mat_id).erode( + MaterialPoints_eroded, + MaterialPoints_stress, + MaterialPoints_pres(mat_point_lid), + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(1, mat_point_lid), + MaterialPoints_sspd(mat_point_lid), + Materials.MaterialFunctions(mat_id).erode_tension_val, + Materials.MaterialFunctions(mat_id).erode_density_val, + mat_point_lid); + + // apply a void eos if mat_point is eroded + double phi_fail = 1.0 - (double)MaterialPoints_eroded(mat_point_lid); + MaterialPoints_pres(mat_point_lid) *= phi_fail; // phi_fail = 1 or 0 + MaterialPoints_sspd(mat_point_lid) *= phi_fail; + MaterialPoints_sspd(mat_point_lid) = fmax(MaterialPoints_sspd(mat_point_lid), 1e-32); + + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < 3; j++) { + MaterialPoints_stress(1, mat_point_lid, i, j) *= phi_fail; + } + } // end for i,j + + }); // end parallel for + } // end if elem errosion return; -} // end method to update state +} // end method to update state2D \ No newline at end of file diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/sgh_setup.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/sgh_setup.cpp new file mode 100644 index 000000000..643b64e4f --- /dev/null +++ b/single-node-refactor/src/Solvers/SGH_solver/src/sgh_setup.cpp @@ -0,0 +1,561 @@ +/********************************************************************************************** +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ + +#include "sgh_solver.h" +#include "region_fill.h" + + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn init_corner_node_masses_zero +/// +/// \brief a function to initialize corner and node masses to zero +/// +/// \param mesh is the simulation mesh +/// \param node_mass is the node mass +/// \param corner_mass is the corner mass +/// +///////////////////////////////////////////////////////////////////////////// +void SGH::init_corner_node_masses_zero(const mesh_t& mesh, + const DCArrayKokkos& node_mass, + const DCArrayKokkos& corner_mass) const +{ + + // calculate the nodal mass + FOR_ALL(node_gid, 0, mesh.num_nodes, { + node_mass(node_gid) = 0.0; + }); // end parallel over nodes + + FOR_ALL(corner_gid, 0, mesh.num_corners, { + corner_mass(corner_gid) = 0.0; + }); // end parallel over corners + +} // end setting masses equal to zero + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn fill_regions_sgh +/// +/// \brief a function to paint den, sie, vel, and mat_ids on the mesh +/// The arrays populated (on host and device) are: +/// elem_mat_id +/// GaussPoint_den +/// GaussPoint_sie +/// node_vel +/// +/// \param Materials holds the material models and global parameters +/// \param mesh is the simulation mesh +/// \param node_coords are the coordinates of the nodes +/// \param node_vel is the nodal velocity array +/// \param region_fills are the instructures to paint state on the mesh +/// \param voxel_elem_mat_id are the voxel values on a structured i,j,k mesh +/// \param GaussPoint_den is density at the GaussPoints on the mesh +/// \param GaussPoint_sie is specific internal energy at the GaussPoints on the mesh +/// \param elem_mat_id is the material id in an element +/// \param num_fills is number of fill instruction +/// \param num_elems is number of elements on the mesh +/// \param num_nodes is number of nodes on the mesh +/// \param rk_num_bins is number of time integration storage bins +/// +///////////////////////////////////////////////////////////////////////////// +void SGH::fill_regions_sgh(const Material_t& Materials, + const mesh_t& mesh, + const DCArrayKokkos & node_coords, + DCArrayKokkos & node_vel, + DCArrayKokkos & GaussPoint_den, + DCArrayKokkos & GaussPoint_sie, + DCArrayKokkos & elem_mat_id, + DCArrayKokkos & voxel_elem_mat_id, + const CArrayKokkos & region_fills, + const CArray & region_fills_host, + const size_t num_fills, + const size_t num_elems, + const size_t num_nodes, + const size_t rk_num_bins) const +{ + + + double voxel_dx, voxel_dy, voxel_dz; // voxel mesh resolution, set by input file + double orig_x, orig_y, orig_z; // origin of voxel elem center mesh, set by input file + size_t voxel_num_i, voxel_num_j, voxel_num_k; // num voxel elements in each direction, set by input file + + + // --------------------------------------------- + // copy to host, enum to read a voxel file + // --------------------------------------------- + + DCArrayKokkos read_voxel_file(num_fills); // check to see if readVoxelFile + + FOR_ALL(f_id, 0, num_fills, { + if (region_fills(f_id).volume == region::readVoxelFile) + { + read_voxel_file(f_id) = region::readVoxelFile; // read the voxel file + } + // add other mesh voxel files + else + { + read_voxel_file(f_id) = 0; + } + }); // end parallel for + read_voxel_file.update_host(); // copy to CPU if code is to read a file + Kokkos::fence(); + // --------------------------------------------- + + + // loop over the fill instructions + for (size_t f_id = 0; f_id < num_fills; f_id++) { + + // ---- + // voxel mesh setup + if (read_voxel_file.host(f_id) == region::readVoxelFile) + { + // read voxel mesh to get the values in the fcn interface + user_voxel_init(voxel_elem_mat_id, + voxel_dx, + voxel_dy, + voxel_dz, + orig_x, + orig_y, + orig_z, + voxel_num_i, + voxel_num_j, + voxel_num_k, + region_fills_host(f_id).scale_x, + region_fills_host(f_id).scale_y, + region_fills_host(f_id).scale_z, + region_fills_host(f_id).file_path); + + // copy values read from file to device + voxel_elem_mat_id.update_device(); + } // endif + // add else if for other mesh reads including STL-2-voxel + + + // parallel loop over elements in mesh + FOR_ALL(elem_gid, 0, num_elems, { + + // calculate the coordinates and radius of the element + double elem_coords_1D[3]; // note:initialization with a list won't work + ViewCArrayKokkos elem_coords(&elem_coords_1D[0], 3); + elem_coords(0) = 0.0; + elem_coords(1) = 0.0; + elem_coords(2) = 0.0; + + // get the coordinates of the element center (using rk_level=1 or node coords) + for (int node_lid = 0; node_lid < mesh.num_nodes_in_elem; node_lid++) { + elem_coords(0) += node_coords(1, mesh.nodes_in_elem(elem_gid, node_lid), 0); + elem_coords(1) += node_coords(1, mesh.nodes_in_elem(elem_gid, node_lid), 1); + if (mesh.num_dims == 3) { + elem_coords(2) += node_coords(1, mesh.nodes_in_elem(elem_gid, node_lid), 2); + } + else{ + elem_coords(2) = 0.0; + } + } // end loop over nodes in element + elem_coords(0) = (elem_coords(0) / mesh.num_nodes_in_elem); + elem_coords(1) = (elem_coords(1) / mesh.num_nodes_in_elem); + elem_coords(2) = (elem_coords(2) / mesh.num_nodes_in_elem); + + + // calc if we are to fill this element + size_t fill_this = fill_geometric_region(mesh, + voxel_elem_mat_id, + region_fills, + elem_coords, + voxel_dx, + voxel_dy, + voxel_dz, + orig_x, + orig_y, + orig_z, + voxel_num_i, + voxel_num_j, + voxel_num_k, + f_id); + + + // paint the material state on the element if fill_this=1 + if (fill_this == 1) { + + // default sgh paint + paint_gauss_den_sie(Materials, + mesh, + node_coords, + GaussPoint_den, + GaussPoint_sie, + elem_mat_id, + region_fills, + elem_coords, + elem_gid, + f_id); + + // add user defined paint here + // user_defined_sgh_state(); + + + // technically, not thread safe, but making it a separate loop created bad fill behavior + // loop over the nodes of this element and apply velocity + for (size_t node_lid = 0; node_lid < mesh.num_nodes_in_elem; node_lid++) { + + // get the mesh node index + size_t node_gid = mesh.nodes_in_elem(elem_gid, node_lid); + + // default sgh paint + paint_node_vel(region_fills, + node_vel, + node_coords, + node_gid, + mesh.num_dims, + f_id, + rk_num_bins); + + // add user defined paint here + // user_defined_vel_state(); + + } // end loop over the nodes in elem + + } // end if fill this + + }); // end FOR_ALL node loop + Kokkos::fence(); + + } // end for loop over fills + + + elem_mat_id.update_host(); + GaussPoint_den.update_host(); + GaussPoint_sie.update_host(); + node_vel.update_host(); + + Kokkos::fence(); + +} // end SGH fill regions + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn setup the SGH method +/// +/// \brief Allocate state, setup models, and fill mesh regions per the YAML input +/// +///////////////////////////////////////////////////////////////////////////// +void SGH::setup(SimulationParameters_t& SimulationParamaters, + Material_t& Materials, + mesh_t& mesh, + BoundaryCondition_t& Boundary, + State_t& State) +{ + + size_t num_fills = SimulationParamaters.region_fills.size(); + printf("Num Fills's = %zu\n", num_fills); + + // the number of elems and nodes in the mesh + const size_t num_elems = mesh.num_elems; + const size_t num_nodes = mesh.num_nodes; + + const size_t rk_num_bins = SimulationParamaters.dynamic_options.rk_num_bins; + + + // create temporary state fields + // Painting routine requires only 1 material per GaussPoint + DCArrayKokkos GaussPoint_den(num_elems); + DCArrayKokkos GaussPoint_sie(num_elems); + DCArrayKokkos elem_mat_id(num_elems); // the mat_id in the elem + + DCArrayKokkos voxel_elem_mat_id; // 1 or 0 if material exist, or it is the material_id + + + + // --------------------------------------------- + // fill den, sie, and velocity on the mesh + // --------------------------------------------- + fill_regions_sgh(Materials, + mesh, + State.node.coords, + State.node.vel, + GaussPoint_den, + GaussPoint_sie, + elem_mat_id, + voxel_elem_mat_id, + SimulationParamaters.region_fills, + SimulationParamaters.region_fills_host, + num_fills, + num_elems, + num_nodes, + rk_num_bins); + + + // note: the device and host side are updated in the above function + // --------------------------------------------- + + + // ---------------------------------------------------------------- + // Walk over the mesh and find dimensions of material storage arrays + // ---------------------------------------------------------------- + const size_t num_mats = Materials.num_mats; // the number of materials on the mesh + + // a counter for the Material index spaces + DCArrayKokkos num_elems_saved_for_mat(num_mats); + + for(int mat_id=0; mat_id (num_mats); + + State.MaterialPoints = CArray (num_mats); + State.MaterialCorners = CArray (num_mats); + // zones not needed with SGH + + + // for ALE SGH, add a buffer to num_elems_for_mat, like 10% of num_elems up to num_elems. + for(int mat_id=0; mat_id corner_areas(&corner_areas_array[0], 4); + // ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 4); + + // geometry::get_area_weights2D(corner_areas, elem_gid, node_coords, elem_node_gids); + + // // loop over the corners of the element and calculate the mass + // for (size_t corner_lid = 0; corner_lid < 4; corner_lid++) { + // size_t corner_gid = mesh.corners_in_elem(elem_gid, corner_lid); + // corner_mass(corner_gid) = corner_areas(corner_lid) * MaterialPoints.den(elem_gid); // node radius is added later + // } // end for over corners + // }); + // + // + // FOR_ALL(nodes_gid=0; nodes_gid ///////////////////////////////////////////////////////////////////////////// /// @@ -50,10 +45,7 @@ void SGH::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) + State_t& State) { std::cout << "In execute function in sgh solver" << std::endl; @@ -96,62 +88,82 @@ void SGH::execute(SimulationParameters_t& SimulationParamaters, CArrayKokkos node_extensive_mass(mesh.num_nodes); + + + std::cout << "Applying initial boundary conditions" << std::endl; + boundary_velocity(mesh, BoundaryConditions, State.node.vel, time_value); // Time value = 0.0; + + + // extensive energy tallies over the entire mesh double IE_t0 = 0.0; double KE_t0 = 0.0; double TE_t0 = 0.0; - double IE_sum = 0.0; - double KE_sum = 0.0; + double cached_pregraphics_dt = fuzz; - double IE_loc_sum = 0.0; - double KE_loc_sum = 0.0; + // calculate the extensive node mass, its key to 2D + calc_extensive_node_mass(node_extensive_mass, + State.node.coords, + State.node.mass, + mesh.num_dims, + mesh.num_nodes); - double cached_pregraphics_dt = fuzz; - // save the nodal mass - FOR_ALL(node_gid, 0, mesh.num_nodes, { - double radius = 1.0; - if (mesh.num_dims == 2) { - radius = node.coords(1, node_gid, 1); - } - node_extensive_mass(node_gid) = node.mass(node_gid) * radius; - }); // end parallel for + // the number of materials specified by the user input + const size_t num_mats = Materials.num_mats; // extensive IE - REDUCE_SUM(elem_gid, 0, mesh.num_elems, IE_loc_sum, { - IE_loc_sum += MaterialPoints.mass(elem_gid) * MaterialPoints.sie(1, elem_gid); - }, IE_sum); - IE_t0 = IE_sum; + for(size_t mat_id=0; mat_id tiny) { - node.mass(node_gid) = node_extensive_mass(node_gid) / node.coords(1, node_gid, 1); - } - }); // end parallel for over node_gid - Kokkos::fence(); - - FOR_ALL(node_bdy_gid, 0, mesh.num_bdy_nodes, { - size_t node_gid = mesh.bdy_nodes(node_bdy_gid); - - if (node.coords(1, node_gid, 1) < tiny) { - // node is on the axis - - for (size_t node_lid = 0; node_lid < mesh.num_nodes_in_node(node_gid); node_lid++) { - size_t node_neighbor_gid = mesh.nodes_in_node(node_gid, node_lid); - - // if the node is off the axis, use it's areal mass on the boundary - if (node.coords(1, node_neighbor_gid, 1) > tiny) { - node.mass(node_gid) = fmax(node.mass(node_gid), node.mass(node_neighbor_gid) / 2.0); - } - } // end for over neighboring nodes - } // end if - }); // end parallel for over elem_gid + calc_node_areal_mass(mesh, + State.node.coords, + State.node.mass, + node_extensive_mass, + tiny); + } // end of if 2D-RZ + + } // end of RK loop // increment the time @@ -408,7 +469,11 @@ void SGH::execute(SimulationParameters_t& SimulationParamaters, // write outputs if (write == 1) { printf("Writing outputs to file at %f \n", graphics_time); - mesh_writer.write_mesh(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters, time_value, graphics_times); + mesh_writer.write_mesh(mesh, + State, + SimulationParamaters, + time_value, + graphics_times); graphics_time = time_value + graphics_dt_ival; @@ -419,6 +484,7 @@ void SGH::execute(SimulationParameters_t& SimulationParamaters, if (time_value >= time_final) { break; } + } // end for cycle loop auto time_2 = std::chrono::high_resolution_clock::now(); @@ -431,41 +497,52 @@ void SGH::execute(SimulationParameters_t& SimulationParamaters, double KE_tend = 0.0; double TE_tend = 0.0; - IE_loc_sum = 0.0; - KE_loc_sum = 0.0; - IE_sum = 0.0; - KE_sum = 0.0; - // extensive IE - REDUCE_SUM(elem_gid, 0, mesh.num_elems, IE_loc_sum, { - IE_loc_sum += MaterialPoints.mass(elem_gid) * MaterialPoints.sie(1, elem_gid); - }, IE_sum); - IE_tend = IE_sum; + for(size_t mat_id=0; mat_id tensor) double abs_max_val = fmax(fabs(eig1), fabs(eig2)); return abs_max_val; } // end 2D max eignen value + + + +void calc_extensive_node_mass(const CArrayKokkos& node_extensive_mass, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass, + double num_dims, + double num_nodes) +{ + // save the nodal mass + FOR_ALL(node_gid, 0, num_nodes, { + + double radius = 1.0; + + if (num_dims == 2) { + radius = node_coords(1, node_gid, 1); + } + + node_extensive_mass(node_gid) = node_mass(node_gid) * radius; + }); // end parallel for +} // end function + + +// a function to tally the internal energy +double sum_domain_internal_energy(const DCArrayKokkos& MaterialPoints_mass, + const DCArrayKokkos& MaterialPoints_sie, + size_t num_mat_points) +{ + + double IE_sum = 0.0; + double IE_loc_sum; + + // loop over the material points and tally IE + REDUCE_SUM(matpt_lid, 0, num_mat_points, IE_loc_sum, { + IE_loc_sum += MaterialPoints_mass(matpt_lid) * MaterialPoints_sie(1,matpt_lid); + }, IE_sum); + Kokkos::fence(); + + + return IE_sum; +} // end function + +double sum_domain_kinetic_energy(const mesh_t& mesh, + const DCArrayKokkos& node_vel, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass) +{ + // extensive KE + double KE_sum = 0.0; + double KE_loc_sum; + + REDUCE_SUM(node_gid, 0, mesh.num_nodes, KE_loc_sum, { + double ke = 0; + + for (size_t dim = 0; dim < mesh.num_dims; dim++) { + ke += node_vel(1, node_gid, dim) * node_vel(1, node_gid, dim); // 1/2 at end + } // end for + + if (mesh.num_dims == 2) { + KE_loc_sum += node_mass(node_gid) * node_coords(1, node_gid, 1) * ke; + } + else{ + KE_loc_sum += node_mass(node_gid) * ke; + } + + }, KE_sum); + Kokkos::fence(); + + + return 0.5*KE_sum; +} // end function + + +// a function to tally the material point masses +double sum_domain_material_mass(const DCArrayKokkos& MaterialPoints_mass, + const size_t num_mat_points) +{ + + double mass_domain = 0.0; + double mass_loc_domain; + + REDUCE_SUM(matpt_lid, 0, num_mat_points, mass_loc_domain, { + + mass_loc_domain += MaterialPoints_mass(matpt_lid); + + }, mass_domain); + Kokkos::fence(); + + return mass_domain; +} // end function + + +double sum_domain_node_mass(const mesh_t& mesh, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass) +{ + + double mass_domain = 0.0; + double mass_loc_domain; + + REDUCE_SUM(node_gid, 0, mesh.num_nodes, mass_loc_domain, { + + if (mesh.num_dims == 2) { + mass_loc_domain += node_mass(node_gid) * node_coords(1, node_gid, 1); + } + else{ + mass_loc_domain += node_mass(node_gid); + } + + }, mass_domain); + Kokkos::fence(); + + + return mass_domain; +} // end function + + +// a function to calculate the 2D-RZ areal mass (rho A = m/R) +// for R=0, it is interpolated from off-axis +void calc_node_areal_mass(const mesh_t& mesh, + const DCArrayKokkos& node_coords, + const DCArrayKokkos& node_mass, + CArrayKokkos node_extensive_mass, + double tiny) +{ + + // calculate the nodal areal mass + FOR_ALL(node_gid, 0, mesh.num_nodes, { + node_mass(node_gid) = 0.0; + + if (node_coords(1, node_gid, 1) > tiny) { + node_mass(node_gid) = node_extensive_mass(node_gid) / node_coords(1, node_gid, 1); + } + }); // end parallel for over node_gid + Kokkos::fence(); + + // calculate the boundary areal mass + FOR_ALL(node_bdy_gid, 0, mesh.num_bdy_nodes, { + size_t node_gid = mesh.bdy_nodes(node_bdy_gid); + + if (node_coords(1, node_gid, 1) < tiny) { + // node is on the axis + + for (size_t node_lid = 0; node_lid < mesh.num_nodes_in_node(node_gid); node_lid++) { + size_t node_neighbor_gid = mesh.nodes_in_node(node_gid, node_lid); + + // if the node is off the axis, use it's areal mass on the boundary + if (node_coords(1, node_neighbor_gid, 1) > tiny) { + node_mass(node_gid) = fmax(node_mass(node_gid), node_mass(node_neighbor_gid) / 2.0); + } + } // end for over neighboring nodes + } // end if + }); // end parallel for over elem_gid + + return; +}// end function + + +// set the corner forces to zero +void set_corner_force_zero(const mesh_t& mesh, + const DCArrayKokkos& corner_force) +{ + + // set corner force to zero + FOR_ALL(corner_gid, 0, mesh.num_corners, { + for (size_t dim = 0; dim < mesh.num_dims; dim++) { + corner_force(corner_gid, dim) = 0.0; + } + }); // end parallel for corners +} // end function diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp index 1cfe57695..4714a9eb7 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/time_integration.cpp @@ -55,18 +55,21 @@ void SGH::rk_init(DCArrayKokkos& node_coords, DCArrayKokkos& MaterialPoints_stress, const size_t num_dims, const size_t num_elems, - const size_t num_nodes) const + const size_t num_nodes, + const size_t num_mat_points) const { + // save elem quantities - FOR_ALL(elem_gid, 0, num_elems, { + FOR_ALL(matpt_lid, 0, num_mat_points, { + // stress is always 3D even with 2D-RZ for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { - MaterialPoints_stress(0, elem_gid, i, j) = MaterialPoints_stress(1, elem_gid, i, j); + MaterialPoints_stress(0, matpt_lid, i, j) = MaterialPoints_stress(1, matpt_lid, i, j); } } // end for - MaterialPoints_sie(0, elem_gid) = MaterialPoints_sie(1, elem_gid); + MaterialPoints_sie(0, matpt_lid) = MaterialPoints_sie(1, matpt_lid); }); // end parallel for // save nodal quantities @@ -100,25 +103,31 @@ void SGH::rk_init(DCArrayKokkos& node_coords, /// ///////////////////////////////////////////////////////////////////////////// void SGH::get_timestep(mesh_t& mesh, - DCArrayKokkos& node_coords, - DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_sspd, - DCArrayKokkos& GaussPoints_vol, - double time_value, - const double graphics_time, - const double time_final, - const double dt_max, - const double dt_min, - const double dt_cfl, - double& dt, - const double fuzz) const + DCArrayKokkos& node_coords, + DCArrayKokkos& node_vel, + DCArrayKokkos& GaussPoints_vol, + DCArrayKokkos& MaterialPoints_sspd, + DCArrayKokkos& MaterialPoints_eroded, + DCArrayKokkos& MaterialToMeshMaps_elem, + size_t num_mat_elems, + double time_value, + const double graphics_time, + const double time_final, + const double dt_max, + const double dt_min, + const double dt_cfl, + double& dt, + const double fuzz) const { // increase dt by 10%, that is the largest dt value dt = dt * 1.1; double dt_lcl; double min_dt_calc; - REDUCE_MIN(elem_gid, 0, mesh.num_elems, dt_lcl, { + REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, { + + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + double coords0[24]; // element coords ViewCArrayKokkos coords(coords0, 8, 3); @@ -168,11 +177,15 @@ void SGH::get_timestep(mesh_t& mesh, } // local dt calc based on CFL - double dt_lcl_ = dt_cfl * dist_min / (MaterialPoints_sspd(elem_gid) + fuzz); + double dt_lcl_ = dt_cfl * dist_min / (MaterialPoints_sspd(mat_elem_lid) + fuzz); + + if (MaterialToMeshMaps_elem(mat_elem_lid) == true){ + dt_lcl_ = 1.0e32; // a huge time step as this element doesn't exist + } // make dt be in bounds - dt_lcl_ = fmin(dt_lcl_, dt_max); // make dt small than dt_max - dt_lcl_ = fmax(dt_lcl_, dt_min); // make dt larger than dt_min + dt_lcl_ = fmin(dt_lcl_, dt_max); // make dt small than dt_max + dt_lcl_ = fmax(dt_lcl_, dt_min); // make dt larger than dt_min if (dt_lcl_ < dt_lcl) { dt_lcl = dt_lcl_; @@ -210,25 +223,31 @@ void SGH::get_timestep(mesh_t& mesh, /// ///////////////////////////////////////////////////////////////////////////// void SGH::get_timestep2D(mesh_t& mesh, - DCArrayKokkos& node_coords, - DCArrayKokkos& node_vel, - DCArrayKokkos& MaterialPoints_sspd, - DCArrayKokkos& GaussPoints_vol, - double time_value, - const double graphics_time, - const double time_final, - const double dt_max, - const double dt_min, - const double dt_cfl, - double& dt, - const double fuzz) const + DCArrayKokkos& node_coords, + DCArrayKokkos& node_vel, + DCArrayKokkos& GaussPoints_vol, + DCArrayKokkos& MaterialPoints_sspd, + DCArrayKokkos& MaterialPoints_eroded, + DCArrayKokkos& MaterialToMeshMaps_elem, + size_t num_mat_elems, + double time_value, + const double graphics_time, + const double time_final, + const double dt_max, + const double dt_min, + const double dt_cfl, + double& dt, + const double fuzz) const { // increase dt by 10%, that is the largest dt value dt = dt * 1.1; double dt_lcl; double min_dt_calc; - REDUCE_MIN(elem_gid, 0, mesh.num_elems, dt_lcl, { + REDUCE_MIN(mat_elem_lid, 0, num_mat_elems, dt_lcl, { + + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + double coords0[8]; // element coords ViewCArrayKokkos coords(coords0, 4, 2); @@ -265,6 +284,11 @@ void SGH::get_timestep2D(mesh_t& mesh, // local dt calc based on CFL double dt_lcl_ = dt_cfl * dist_min / (MaterialPoints_sspd(elem_gid) + fuzz); + + if (MaterialToMeshMaps_elem(mat_elem_lid) == true){ + dt_lcl_ = 1.0e32; // a huge time step as this element doesn't exist + } + // make dt be in bounds dt_lcl_ = fmin(dt_lcl_, dt_max); // make dt small than dt_max dt_lcl_ = fmax(dt_lcl_, dt_min); // make dt larger than dt_min diff --git a/single-node-refactor/src/common/dynamic_options.h b/single-node-refactor/src/common/dynamic_options.h index 28be0f198..b9cbb69a7 100644 --- a/single-node-refactor/src/common/dynamic_options.h +++ b/single-node-refactor/src/common/dynamic_options.h @@ -31,9 +31,9 @@ WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. **********************************************************************************************/ - #ifndef FIERRO_DYNAMIC_OPTIONS_H #define FIERRO_DYNAMIC_OPTIONS_H + #include #include "matar.h" diff --git a/single-node-refactor/src/common/io_utils.h b/single-node-refactor/src/common/io_utils.h index c39d4769a..0f94daada 100644 --- a/single-node-refactor/src/common/io_utils.h +++ b/single-node-refactor/src/common/io_utils.h @@ -38,12 +38,73 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "mesh.h" #include "state.h" #include "simulation_parameters.h" +#include "region.h" #include #include #include #include + + + + +// ============================================================================== +// Functions to get 1D for an i,j,k layout mesh +// ============================================================================== +// inline int get_id(int i, int j, int k, int num_i, int num_j); + +// KOKKOS_INLINE_FUNCTION +// int get_id_device(int i, int j, int k, int num_i, int num_j); + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn get_id +/// +/// \brief This gives the index value of the point or the elem +/// +/// Assumes that the grid has an i,j,k structure +/// the elem = i + (j)*(num_points_i-1) + (k)*(num_points_i-1)*(num_points_j-1) +/// the point = i + (j)*num_points_i + (k)*num_points_i*num_points_j +/// +/// \param i index +/// \param j index +/// \param k index +/// \param Number of i indices +/// \param Number of j indices +/// +///////////////////////////////////////////////////////////////////////////// +inline int get_id(int i, int j, int k, int num_i, int num_j) +{ + return i + j * num_i + k * num_i * num_j; +} + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn get_id_device +/// +/// \brief This gives the index value of the point or the elem +/// +/// Assumes that the grid has an i,j,k structure +/// the elem = i + (j)*(num_points_i-1) + (k)*(num_points_i-1)*(num_points_j-1) +/// the point = i + (j)*num_points_i + (k)*num_points_i*num_points_j +/// +/// \param i index +/// \param j index +/// \param k index +/// \param Number of i indices +/// \param Number of j indices +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_INLINE_FUNCTION +int get_id_device(int i, int j, int k, int num_i, int num_j) +{ + return i + j * num_i + k * num_i * num_j; +} + + + ///////////////////////////////////////////////////////////////////////////// /// /// \class MeshReader @@ -97,7 +158,6 @@ class MeshReader /// ///////////////////////////////////////////////////////////////////////////// void read_mesh(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -111,7 +171,7 @@ class MeshReader // Check mesh file extension // and read based on extension - read_ensight_mesh(mesh, MaterialPoints, GaussPoints, node, corner, num_dims, rk_num_bins); + read_ensight_mesh(mesh, GaussPoints, node, corner, num_dims, rk_num_bins); } ///////////////////////////////////////////////////////////////////////////// @@ -129,7 +189,6 @@ class MeshReader /// ///////////////////////////////////////////////////////////////////////////// void read_ensight_mesh(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -208,7 +267,6 @@ class MeshReader // initialize elem variables mesh.initialize_elems(num_elem, num_dims); - MaterialPoints.initialize(rk_num_bins, num_elem, 3); // always 3D here, even for 2D GaussPoints.initialize(rk_num_bins, num_elem, 3); // always 3D here, even for 2D // for each cell read the list of associated nodes @@ -282,7 +340,6 @@ class MeshBuilder /// ///////////////////////////////////////////////////////////////////////////// void build_mesh(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -290,10 +347,10 @@ class MeshBuilder { if (SimulationParamaters.mesh_input.num_dims == 2) { if (SimulationParamaters.mesh_input.type == mesh_input::Cylinder) { - build_2d_polar(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters); + build_2d_polar(mesh, GaussPoints, node, corner, SimulationParamaters); } else if (SimulationParamaters.mesh_input.type == mesh_input::Box) { - build_2d_box(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters); + build_2d_box(mesh, GaussPoints, node, corner, SimulationParamaters); } else{ std::cout << "**** 2D MESH TYPE NOT SUPPORTED **** " << std::endl; @@ -306,7 +363,7 @@ class MeshBuilder } } else if (SimulationParamaters.mesh_input.num_dims == 3) { - build_3d_box(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters); + build_3d_box(mesh, GaussPoints, node, corner, SimulationParamaters); } else{ throw std::runtime_error("**** ONLY 2D RZ OR 3D MESHES ARE SUPPORTED ****"); @@ -327,7 +384,6 @@ class MeshBuilder /// ///////////////////////////////////////////////////////////////////////////// void build_2d_box(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -402,7 +458,6 @@ class MeshBuilder // intialize elem variables mesh.initialize_elems(num_elems, num_dim); - MaterialPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D GaussPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D // populate the elem center data structures @@ -460,7 +515,6 @@ class MeshBuilder /// ///////////////////////////////////////////////////////////////////////////// void build_2d_polar(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -538,7 +592,6 @@ class MeshBuilder // intialize elem variables mesh.initialize_elems(num_elems, num_dim); - MaterialPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D GaussPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D // populate the elem center data structures @@ -596,7 +649,6 @@ class MeshBuilder /// ///////////////////////////////////////////////////////////////////////////// void build_3d_box(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -685,7 +737,6 @@ class MeshBuilder // intialize elem variables mesh.initialize_elems(num_elems, num_dim); - MaterialPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D GaussPoints.initialize(rk_num_bins, num_elems, 3); // always 3D here, even for 2D // --- Build elems --- @@ -749,7 +800,6 @@ class MeshBuilder /// ///////////////////////////////////////////////////////////////////////////// void build_3d_HexN_box(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, GaussPoint_t& GaussPoints, node_t& node, corner_t& corner, @@ -758,27 +808,7 @@ class MeshBuilder printf(" ***** WARNING:: build_3d_HexN_box not yet implemented\n"); } - ///////////////////////////////////////////////////////////////////////////// - /// - /// \fn get_id - /// - /// \brief This gives the index value of the point or the elem - /// - /// Assumes that the grid has an i,j,k structure - /// the elem = i + (j)*(num_points_i-1) + (k)*(num_points_i-1)*(num_points_j-1) - /// the point = i + (j)*num_points_i + (k)*num_points_i*num_points_j - /// - /// \param i index - /// \param j index - /// \param k index - /// \param Number of i indices - /// \param Number of j indices - /// - ///////////////////////////////////////////////////////////////////////////// - int get_id(int i, int j, int k, int num_i, int num_j) const - { - return i + j * num_i + k * num_i * num_j; - } + ///////////////////////////////////////////////////////////////////////////// /// @@ -886,19 +916,24 @@ class MeshWriter /// ///////////////////////////////////////////////////////////////////////////// void write_mesh(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - node_t& node, - corner_t& corner, + State_t& State, SimulationParameters_t& SimulationParamaters, double time_value, CArray graphics_times) { if (SimulationParamaters.output_options.format == output_options::vtk) { - write_vtk(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters, time_value, graphics_times); + write_vtk(mesh, + State, + SimulationParamaters, + time_value, + graphics_times); } else if (SimulationParamaters.output_options.format == output_options::ensight) { - write_ensight(mesh, MaterialPoints, GaussPoints, node, corner, SimulationParamaters, time_value, graphics_times); + write_ensight(mesh, + State, + SimulationParamaters, + time_value, + graphics_times); } else{ std::cout << "**** MESH OUTPUT TYPE NOT SUPPORTED **** " << std::endl; @@ -928,30 +963,42 @@ class MeshWriter /// ///////////////////////////////////////////////////////////////////////////// void write_ensight(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - node_t& node, - corner_t& corner, + State_t& State, SimulationParameters_t& SimulationParamaters, double time_value, CArray graphics_times) { - // Update host data - MaterialPoints.den.update_host(); - MaterialPoints.pres.update_host(); - MaterialPoints.stress.update_host(); - MaterialPoints.sspd.update_host(); - MaterialPoints.sie.update_host(); - GaussPoints.vol.update_host(); - MaterialPoints.mass.update_host(); - GaussPoints.mat_id.update_host(); - - node.coords.update_host(); - node.vel.update_host(); - node.mass.update_host(); + + size_t num_mats = State.MaterialPoints.size(); + + // ---- Update host data ---- + + // material point values + for(int mat_id=0; mat_id vec_fields(num_nodes, num_vec_vars, 3); for (size_t node_gid = 0; node_gid < num_nodes; node_gid++) { // position, var 0 - vec_fields(node_gid, 0, 0) = node.coords.host(1, node_gid, 0); - vec_fields(node_gid, 0, 1) = node.coords.host(1, node_gid, 1); + vec_fields(node_gid, 0, 0) = State.node.coords.host(1, node_gid, 0); + vec_fields(node_gid, 0, 1) = State.node.coords.host(1, node_gid, 1); if (num_dims == 2) { vec_fields(node_gid, 0, 2) = 0.0; } else{ - vec_fields(node_gid, 0, 2) = node.coords.host(1, node_gid, 2); + vec_fields(node_gid, 0, 2) = State.node.coords.host(1, node_gid, 2); } // position, var 1 - vec_fields(node_gid, 1, 0) = node.vel.host(1, node_gid, 0); - vec_fields(node_gid, 1, 1) = node.vel.host(1, node_gid, 1); + vec_fields(node_gid, 1, 0) = State.node.vel.host(1, node_gid, 0); + vec_fields(node_gid, 1, 1) = State.node.vel.host(1, node_gid, 1); if (num_dims == 2) { vec_fields(node_gid, 1, 2) = 0.0; } else{ - vec_fields(node_gid, 1, 2) = node.vel.host(1, node_gid, 2); + vec_fields(node_gid, 1, 2) = State.node.vel.host(1, node_gid, 2); } } // end for loop over vertices @@ -1094,16 +1165,16 @@ class MeshWriter // write all components of the point coordinates for (int node_gid = 0; node_gid < num_nodes; node_gid++) { - fprintf(out[0], "%12.5e\n", node.coords.host(1, node_gid, 0)); + fprintf(out[0], "%12.5e\n", State.node.coords.host(1, node_gid, 0)); } for (int node_gid = 0; node_gid < num_nodes; node_gid++) { - fprintf(out[0], "%12.5e\n", node.coords.host(1, node_gid, 1)); + fprintf(out[0], "%12.5e\n", State.node.coords.host(1, node_gid, 1)); } for (int node_gid = 0; node_gid < num_nodes; node_gid++) { if (num_dims == 3) { - fprintf(out[0], "%12.5e\n", node.coords.host(1, node_gid, 2)); + fprintf(out[0], "%12.5e\n", State.node.coords.host(1, node_gid, 2)); } else{ fprintf(out[0], "%12.5e\n", 0.0); @@ -1261,7 +1332,7 @@ class MeshWriter ///////////////////////////////////////////////////////////////////////////// /// - /// \fn write_ensight + /// \fn write_vtk /// /// \brief Writes a vtk output file /// @@ -1275,10 +1346,7 @@ class MeshWriter /// ///////////////////////////////////////////////////////////////////////////// void write_vtk(mesh_t& mesh, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - node_t& node, - corner_t& corner, + State_t& State, SimulationParameters_t& SimulationParamaters, double time_value, CArray graphics_times) @@ -1286,6 +1354,8 @@ class MeshWriter // Not yet supported throw std::runtime_error("**** VTK OUTPUT TYPE NOT YET SUPPORTED ****"); } -}; +}; // end class + + #endif // end Header Guard \ No newline at end of file diff --git a/single-node-refactor/src/common/material.h b/single-node-refactor/src/common/material.h index 2b82de1c9..85ec0f9db 100644 --- a/single-node-refactor/src/common/material.h +++ b/single-node-refactor/src/common/material.h @@ -170,6 +170,9 @@ struct MaterialEnums_t // Strength model type: none, or increment- or state-based model::StrengthType StrengthType = model::noStrengthType; + // Erosion model type: none or basis + model::ErosionModels ErosionModels = model::noErosion; + }; // end boundary condition enums @@ -190,7 +193,7 @@ struct MaterialFunctions_t // Equation of state (EOS) function pointers void (*calc_pressure)(const DCArrayKokkos& MaterialPoints_pres, const DCArrayKokkos& MaterialPoints_stress, - const size_t MaterialPoints_gid, + const size_t MaterialPoints_lid, const size_t mat_id, const DCArrayKokkos& MaterialPoints_state_vars, const DCArrayKokkos& MaterialPoints_sspd, @@ -200,7 +203,7 @@ struct MaterialFunctions_t void (*calc_sound_speed)(const DCArrayKokkos& MaterialPoints_pres, const DCArrayKokkos& MaterialPoints_stress, - const size_t MaterialPoints_gid, + const size_t MaterialPoints_lid, const size_t mat_id, const DCArrayKokkos& MaterialPoints_state_vars, const DCArrayKokkos& MaterialPoints_sspd, @@ -208,19 +211,20 @@ struct MaterialFunctions_t const double sie, const RaggedRightArrayKokkos &eos_global_vars) = NULL; + // -- Strength -- // Material strength model function pointers void (*calc_stress)(const DCArrayKokkos& MaterialPoints_pres, const DCArrayKokkos& MaterialPoints_stress, - const size_t MaterialPoints_gid, + const size_t MaterialPoints_lid, const size_t mat_id, const DCArrayKokkos& MaterialPoints_state_vars, const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const ViewCArrayKokkos& vel_grad, - const ViewCArrayKokkos& MaterialPoints_node_gids, + const ViewCArrayKokkos& elem_node_gids, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, const double vol, @@ -230,21 +234,18 @@ struct MaterialFunctions_t // -- Erosion -- - size_t void_mat_id; ///< eroded elements get this mat_id double erode_tension_val; ///< tension threshold to initiate erosion double erode_density_val; ///< density threshold to initiate erosion // above should be removed, they go in CArrayKokkos erosion_global_vars; - void (*erode)(const DCArrayKokkos& MaterialPoints_pres, + void (*erode)(const DCArrayKokkos& MaterialPoints_eroded, const DCArrayKokkos& MaterialPoints_stress, - const DCArrayKokkos& MaterialPoints_eroded, - const DCArrayKokkos& MaterialPoints_mat_id, - const size_t MaterialPoints_gid, - const size_t void_mat_id, + const double MaterialPoint_pres, + const double MaterialPoint_den, + const double MaterialPoint_sie, + const double MaterialPoint_sspd, const double erode_tension_val, const double erode_density_val, - const DCArrayKokkos& MaterialPoints_sspd, - const DCArrayKokkos& MaterialPoints_den, - const double sie) = NULL; + const size_t mat_point_lid) = NULL; double q1 = 1.0; ///< acoustic coefficient in Riemann solver for compression double q1ex = 1.3333; ///< acoustic coefficient in Riemann solver for expansion @@ -282,14 +283,8 @@ struct Material_t{ RaggedRightArrayKokkos eos_global_vars; ///< Array of global variables for the EOS CArrayKokkos num_eos_global_vars; - RaggedRightArrayKokkos eos_state_vars; ///< Array of state (in each element) variables for the EOS - CArrayKokkos num_eos_state_vars; - RaggedRightArrayKokkos strength_global_vars; ///< Array of global variables for the strength model CArrayKokkos num_strength_global_vars; - - RaggedRightArrayKokkos strength_state_vars; ///< Array of state (in each element) variables for the strength - CArrayKokkos num_strength_state_vars; RaggedRightArrayKokkos failure_global_vars; ///< Array of global variables for the failure model diff --git a/single-node-refactor/src/common/mesh.h b/single-node-refactor/src/common/mesh.h index a491619ed..6ef4448a4 100644 --- a/single-node-refactor/src/common/mesh.h +++ b/single-node-refactor/src/common/mesh.h @@ -240,7 +240,7 @@ struct mesh_t size_t num_lob_gauss_in_elem; ///< Number of Gauss Lobatto points in an element DCArrayKokkos nodes_in_elem; ///< Nodes in an element - CArrayKokkos corners_in_elem; ///< Corners in an element + CArrayKokkos corners_in_elem; ///< Corners in an element -- this can just be a functor RaggedRightArrayKokkos elems_in_elem; ///< Elements connected to an element CArrayKokkos num_elems_in_elem; ///< Number of elements connected to an element @@ -407,7 +407,7 @@ struct mesh_t size_t j = count_saved_corners_in_node(node_gid); // Save corner index to this node_gid - size_t corner_gid = node_lid + elem_gid * num_nodes_in_elem; + size_t corner_gid = node_lid + elem_gid * num_nodes_in_elem; // this can be a functor corners_in_node(node_gid, j) = corner_gid; elems_in_node(node_gid, j) = elem_gid; // save the elem_gid diff --git a/single-node-refactor/src/common/region.h b/single-node-refactor/src/common/region.h index 7c8a83da7..d589a440d 100644 --- a/single-node-refactor/src/common/region.h +++ b/single-node-refactor/src/common/region.h @@ -96,7 +96,7 @@ struct reg_fill_t double radius2 = 0.0; ///< Outer radius to fill for sphere // initial conditions - init_conds::init_velocity_conds velocity; ///< Initial conditions for this region WARNING: Currently unimplemented + init_conds::init_velocity_conds velocity; ///< Initial conditions for this region // velocity coefficients by component double u = 0.0; ///< U component of velocity diff --git a/single-node-refactor/src/common/region_fill.h b/single-node-refactor/src/common/region_fill.h new file mode 100644 index 000000000..3dd96ed8c --- /dev/null +++ b/single-node-refactor/src/common/region_fill.h @@ -0,0 +1,922 @@ +/********************************************************************************************** +� 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ +#ifndef REGION_FILL_H +#define REGION_FILL_H + +#include "matar.h" +#include "mesh.h" +#include "material.h" +#include "state.h" +#include "simulation_parameters.h" +#include "region.h" +#include "io_utils.h" + +#include + +// ============================================================================== +// Functions to read voxel mesh +// ============================================================================== +void user_voxel_init(DCArrayKokkos& elem_values, + double& dx, + double& dy, + double& dz, + double& orig_x, + double& orig_y, + double& orig_z, + size_t& voxel_num_i, + size_t& voxel_num_j, + size_t& voxel_num_k, + double scale_x, + double scale_y, + double scale_z, + std::string mesh_file); + +// ----------------------------------------------------------------------------- +// The function to read a voxel vtk file from Dream3d and intialize the mesh +// ------------------------------------------------------------------------------ +void user_voxel_init(DCArrayKokkos& elem_values, + double& dx, + double& dy, + double& dz, + double& orig_x, + double& orig_y, + double& orig_z, + size_t& num_elems_i, + size_t& num_elems_j, + size_t& num_elems_k, + double scale_x, + double scale_y, + double scale_z, + std::string mesh_file) +{ + std::string MESH = mesh_file; // user specified + + std::ifstream in; // FILE *in; + in.open(MESH); + + // check to see of a mesh was supplied when running the code + if (in) + { + printf("\nReading the 3D voxel mesh: "); + std::cout << MESH << std::endl; + } + else + { + std::cout << "\n\n**********************************\n\n"; + std::cout << " ERROR:\n"; + std::cout << " Voxel vtk input does not exist \n"; + std::cout << "**********************************\n\n" << std::endl; + std::exit(EXIT_FAILURE); + } // end if + + size_t i; // used for writing information to file + size_t point_id; // the global id for the point + size_t elem_id; // the global id for the elem + size_t this_point; // a local id for a point in a elem (0:7 for a Hexahedral elem) + + size_t num_points_i; + size_t num_points_j; + size_t num_points_k; + + size_t num_dims = 3; + + std::string token; + + bool found = false; + + // look for POINTS + i = 0; + while (found == false) { + std::string str; + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split(str, delimiter); + + // looking for the following text: + // POINTS %d float + if (v[0] == "DIMENSIONS") + { + num_points_i = std::stoi(v[1]); + num_points_j = std::stoi(v[2]); + num_points_k = std::stoi(v[3]); + printf("Num voxel nodes read in = %zu, %zu, %zu\n", num_points_i, num_points_j, num_points_k); + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find POINTS \n"); + break; + } // end if + + i++; + } // end while + + found = false; + + int num_points = num_points_i * num_points_j * num_points_k; + CArray pt_coords_x(num_points_i); + CArray pt_coords_y(num_points_j); + CArray pt_coords_z(num_points_k); + + while (found == false) { + std::string str; + std::string str0; + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split(str, delimiter); + + // looking for the following text: + if (v[0] == "X_COORDINATES") + { + size_t num_saved = 0; + + while (num_saved < num_points_i - 1) { + // get next line + std::getline(in, str0); + + // remove starting and trailing spaces + str = trim(str0); + std::vector v_coords = split(str, delimiter); + + // loop over the contents of the vector v_coords + for (size_t this_point = 0; this_point < v_coords.size(); this_point++) + { + pt_coords_x(num_saved) = scale_x*std::stod(v_coords[this_point]); + num_saved++; + } // end for + } // end while + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find X_COORDINATES \n"); + break; + } // end if + + i++; + } // end while + found = false; + + while (found == false) { + std::string str; + std::string str0; + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split(str, delimiter); + + // looking for the following text: + if (v[0] == "Y_COORDINATES") + { + size_t num_saved = 0; + + while (num_saved < num_points_j - 1) { + // get next line + std::getline(in, str0); + + // remove starting and trailing spaces + str = trim(str0); + std::vector v_coords = split(str, delimiter); + + // loop over the contents of the vector v_coords + for (size_t this_point = 0; this_point < v_coords.size(); this_point++) + { + pt_coords_y(num_saved) = scale_y*std::stod(v_coords[this_point]); + num_saved++; + } // end for + } // end while + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find Y_COORDINATES \n"); + break; + } // end if + + i++; + } // end while + found = false; + + while (found == false) { + std::string str; + std::string str0; + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split(str, delimiter); + + // looking for the following text: + if (v[0] == "Z_COORDINATES") + { + size_t num_saved = 0; + + while (num_saved < num_points_k - 1) { + // get next line + std::getline(in, str0); + + // remove starting and trailing spaces + str = trim(str0); + std::vector v_coords = split(str, delimiter); + + // loop over the contents of the vector v_coords + for (size_t this_point = 0; this_point < v_coords.size(); this_point++) + { + pt_coords_z(num_saved) = scale_z*std::stod(v_coords[this_point]); + num_saved++; + } // end for + } // end while + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find Z_COORDINATES \n"); + break; + } // end if + + i++; + } // end while + found = false; + + size_t num_elems; + num_elems_i = num_points_i - 1; + num_elems_j = num_points_j - 1; + num_elems_k = num_points_k - 1; + + // center to center distance between first and last elem along each edge + double Lx = (pt_coords_x(num_points_i - 2) - pt_coords_x(0)); + double Ly = (pt_coords_y(num_points_j - 2) - pt_coords_y(0)); + double Lz = (pt_coords_z(num_points_k - 2) - pt_coords_z(0)); + + // spacing between elems + dx = Lx / ((double) num_elems_i); + dy = Ly / ((double) num_elems_j); + dz = Lz / ((double) num_elems_k); + + // element mesh origin + orig_x = 0.5 * (pt_coords_x(0) + pt_coords_x(1)), + orig_y = 0.5 * (pt_coords_y(0) + pt_coords_y(1)), + orig_z = 0.5 * (pt_coords_z(0) + pt_coords_z(1)), + + // look for CELLS + i = 0; + while (found == false) { + std::string str; + std::getline(in, str); + + std::string delimiter = " "; + std::vector v = split(str, delimiter); + + // looking for the following text: + // CELLS num_elems size + if (v[0] == "CELL_DATA") + { + num_elems = std::stoi(v[1]); + printf("Num voxel elements read in %zu\n", num_elems); + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find CELL_DATA \n"); + break; + } // end if + + i++; + } // end while + found = false; + + // allocate memory for element voxel values + elem_values = DCArrayKokkos(num_elems); + + // reading the cell data + while (found == false) { + std::string str; + std::string str0; + + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split(str, delimiter); + + // looking for the following text: + if (v[0] == "LOOKUP_TABLE") + { + size_t num_saved = 0; + + while (num_saved < num_elems - 1) { + // get next line + std::getline(in, str0); + + // remove starting and trailing spaces + str = trim(str0); + std::vector v_values = split(str, delimiter); + + // loop over the contents of the vector v_coords + for (size_t this_elem = 0; this_elem < v_values.size(); this_elem++) + { + // save integers (0 or 1) to host side + elem_values.host(num_saved) = std::stoi(v_values[this_elem]); + num_saved++; + } // end for + + // printf(" done with one row of data \n"); + } // end while + + found = true; + } // end if + + if (i > 1000) + { + printf("ERROR: Failed to find LOOKUP_TABLE data \n"); + break; + } // end if + + i++; + } // end while + found = false; + + printf("\n"); + + in.close(); +} // end routine +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn fill_geometric_region +/// +/// \brief a function to calculate whether to fill this element based on the +/// input instructions. The output is +/// = 0 then no, do not fill this element +/// = 1 then yes, fill this element +/// +/// \param mesh is the simulation mesh +/// \param node_coords is the nodal position array +/// \param voxel_elem_mat_id are the voxel values on a structured i,j,k mesh +/// \param region_fills are the instructures to paint state on the mesh +/// \param mesh_coords is the geometric center of the element or a node coordinates +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +size_t fill_geometric_region(const mesh_t& mesh, + const DCArrayKokkos& voxel_elem_mat_id, + const CArrayKokkos& region_fills, + const ViewCArrayKokkos & mesh_coords, + const double voxel_dx, + const double voxel_dy, + const double voxel_dz, + const double orig_x, + const double orig_y, + const double orig_z, + const size_t voxel_num_i, + const size_t voxel_num_j, + const size_t voxel_num_k, + const size_t f_id){ + + // default is not to fill the element + size_t fill_this = 0; + + + // for shapes with an origin (e.g., sphere and circle), accounting for the origin + double dist_x = mesh_coords(0) - region_fills(f_id).origin[0]; + double dist_y = mesh_coords(1) - region_fills(f_id).origin[1]; + double dist_z = mesh_coords(2) - region_fills(f_id).origin[2]; + + // spherical radius + double radius = sqrt(dist_x * dist_x + + dist_y * dist_y + + dist_z * dist_z); + + // cylindrical radius + double radius_cyl = sqrt(dist_x * dist_x + + dist_y * dist_y); + + + // check to see if this element should be filled + switch (region_fills(f_id).volume) { + case region::global: + { + fill_this = 1; + break; + } + case region::box: + { + + double x_lower_bound = region_fills(f_id).x1; + double x_upper_bound = region_fills(f_id).x2; + + double y_lower_bound = region_fills(f_id).y1; + double y_upper_bound = region_fills(f_id).y2; + + double z_lower_bound = region_fills(f_id).z1; + double z_upper_bound = region_fills(f_id).z2; + + + if (mesh_coords(0) >= x_lower_bound && mesh_coords(0) <= x_upper_bound && + mesh_coords(1) >= y_lower_bound && mesh_coords(1) <= y_upper_bound && + mesh_coords(2) >= z_lower_bound && mesh_coords(2) <= z_upper_bound) { + fill_this = 1; + } + break; + } + case region::cylinder: + { + if (radius_cyl >= region_fills(f_id).radius1 + && radius_cyl <= region_fills(f_id).radius2) { + fill_this = 1; + } + break; + } + case region::sphere: + { + if (radius >= region_fills(f_id).radius1 + && radius <= region_fills(f_id).radius2) { + fill_this = 1; + } + break; + } + + case region::readVoxelFile: + { + + fill_this = 0; // default is no, don't fill it + + // find the closest element in the voxel mesh to this element + double i0_real = (mesh_coords(0) - orig_x - region_fills(f_id).origin[0]) / (voxel_dx); + double j0_real = (mesh_coords(1) - orig_y - region_fills(f_id).origin[1]) / (voxel_dy); + double k0_real = (mesh_coords(2) - orig_z - region_fills(f_id).origin[2]) / (voxel_dz); + + int i0 = (int)i0_real; + int j0 = (int)j0_real; + int k0 = (int)k0_real; + + // look for the closest element in the voxel mesh + int elem_id0 = get_id_device(i0, j0, k0, voxel_num_i, voxel_num_j); + + // if voxel mesh overlaps this mesh, then fill it if =1 + if (elem_id0 < voxel_elem_mat_id.size() && elem_id0 >= 0 && + i0 >= 0 && j0 >= 0 && k0 >= 0 && + i0 < voxel_num_i && j0 < voxel_num_j && k0 < voxel_num_k) { + // voxel mesh elem values = 0 or 1 + fill_this = voxel_elem_mat_id(elem_id0); // values from file + + } // end if + + break; + + } // end case + case region::no_volume: + { + fill_this = 0; // default is no, don't fill it + + break; + } + default: + { + fill_this = 0; // default is no, don't fill it + + break; + } + + } // end of switch + + + return fill_this; + +} // end function + + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn paint_gauss_den_sie +/// +/// \brief a function to paint den and sie on the Gauss points of the mesh +/// +/// \param Materials holds the material models and global parameters +/// \param mesh is the simulation mesh +/// \param node_coords are the node coordinates of the element +/// \param GaussPoint_den is density at the GaussPoints on the mesh +/// \param GaussPoint_sie is specific internal energy at the GaussPoints on the mesh +/// \param elem_mat_id is the material id in an element +/// \param region_fills are the instructures to paint state on the mesh +/// \param elem_coords is the geometric center of the element +/// \param elem_gid is the element global mesh index +/// \param f_id is fill instruction +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +void paint_gauss_den_sie(const Material_t& Materials, + const mesh_t& mesh, + const DCArrayKokkos & node_coords, + const DCArrayKokkos & GaussPoint_den, + const DCArrayKokkos & GaussPoint_sie, + const DCArrayKokkos & elem_mat_id, + const CArrayKokkos& region_fills, + const ViewCArrayKokkos elem_coords, + const double elem_gid, + const size_t f_id){ + + // the material id + size_t mat_id = region_fills(f_id).material_id; + + // --- material_id in elem --- + elem_mat_id(elem_gid) = mat_id; + + // loop over the Gauss points in the element + { + + const size_t gauss_gid = elem_gid; // 1 gauss point per element + + // add test problem state setups here + if (region_fills(f_id).velocity == init_conds::tg_vortex) { + + GaussPoint_den(gauss_gid) = 1.0; + + // note: elem_coords are the gauss_coords, higher quadrature requires ref elem data + double pres = 0.25 * (cos(2.0 * PI * elem_coords(0)) + + cos(2.0 * PI * elem_coords(1)) ) + 1.0; + + // p = rho*ie*(gamma - 1) + // makes sure index 0 matches the gamma in the gamma law function + double gamma = Materials.eos_global_vars(mat_id,0); + GaussPoint_sie(gauss_gid) = + pres / (GaussPoint_den(gauss_gid) * (gamma - 1.0)); + } // end + // add user initialization here + else{ + + // --- density --- + GaussPoint_den(gauss_gid) = region_fills(f_id).den; + + // --- specific internal energy --- + GaussPoint_sie(gauss_gid) = region_fills(f_id).sie; + + } // end if + + } // end loop over gauss points in element' + + // done setting the element state + +} // end function + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn paint_node_vel +/// +/// \brief a function to paint a velocity field on the nodes of the mesh +/// +/// \param mesh is the simulation mesh +/// \param node_vel is the nodal velocity array +/// \param node_coords are the coordinates of the nodes +/// \param elem_gid is the element global mesh index +/// \param f_id is fill instruction +/// \param rk_num_bins is time integration storage level +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +void paint_node_vel(const CArrayKokkos& region_fills, + const DCArrayKokkos& node_vel, + const DCArrayKokkos& node_coords, + const double node_gid, + const double num_dims, + const size_t f_id, + const size_t rk_num_bins){ + + // save velocity at all rk_levels + for(size_t rk_level=0; rk_level 1.0e-14) { + dir[dim] /= (radius_val); + } + else{ + dir[dim] = 0.0; + } + } // end for + + node_vel(rk_level, node_gid, 0) = region_fills(f_id).speed * dir[0]; + node_vel(rk_level, node_gid, 1) = region_fills(f_id).speed * dir[1]; + if (num_dims == 3) { + node_vel(rk_level, node_gid, 2) = 0.0; + } + + break; + } + case init_conds::spherical: + { + // Setting up spherical + double dir[3]; + dir[0] = 0.0; + dir[1] = 0.0; + dir[2] = 0.0; + double radius_val = 0.0; + + for (int dim = 0; dim < 3; dim++) { + dir[dim] = node_coords(rk_level, node_gid, dim); + radius_val += node_coords(rk_level, node_gid, dim) * node_coords(rk_level, node_gid, dim); + } // end for + radius_val = sqrt(radius_val); + + for (int dim = 0; dim < 3; dim++) { + if (radius_val > 1.0e-14) { + dir[dim] /= (radius_val); + } + else{ + dir[dim] = 0.0; + } + } // end for + + node_vel(rk_level, node_gid, 0) = region_fills(f_id).speed * dir[0]; + node_vel(rk_level, node_gid, 1) = region_fills(f_id).speed * dir[1]; + if (num_dims == 3) { + node_vel(rk_level, node_gid, 2) = region_fills(f_id).speed * dir[2]; + } + + break; + } + case init_conds::radial_linear: + { + printf("**** Radial_linear initial conditions not yet supported ****\n"); + break; + } + case init_conds::spherical_linear: + { + printf("**** spherical_linear initial conditions not yet supported ****\n"); + break; + } + case init_conds::tg_vortex: + { + node_vel(rk_level, node_gid, 0) = sin(PI * node_coords(rk_level, node_gid, 0)) * + cos(PI * node_coords(rk_level, node_gid, 1)); + node_vel(rk_level, node_gid, 1) = -1.0 * cos(PI * node_coords(rk_level, node_gid, 0)) * + sin(PI * node_coords(rk_level, node_gid, 1)); + if (num_dims == 3) { + node_vel(rk_level, node_gid, 2) = 0.0; + } + + break; + } + + case init_conds::no_ic_vel: + { + // no velocity + node_vel(rk_level, node_gid, 0) = 0.0; + node_vel(rk_level, node_gid, 1) = 0.0; + if (num_dims == 3) { + node_vel(rk_level, node_gid, 2) = 0.0; + } + + break; + } + default: + { + // no velocity + node_vel(rk_level, node_gid, 0) = 0.0; + node_vel(rk_level, node_gid, 1) = 0.0; + if (num_dims == 3) { + node_vel(rk_level, node_gid, 2) = 0.0; + } + + break; + } + } // end of switch + + } // end loop over rk_num_bins + + + // done setting the velocity +} // end function paint_node_vel + + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn init_press_sspd_stress +/// +/// \brief a function to initialize pressure, sound speed and stress +/// +/// \param Materials holds the material models and global parameters +/// \param mesh is the simulation mesh +/// \param GaussPoint_den is density at the GaussPoints on the mesh +/// \param GaussPoint_pres is pressure at the GaussPoints on the mesh +/// \param GaussPoint_stress is stress at the GaussPoints on the mesh +/// \param GaussPoint_sspd is sound speed at the GaussPoints on the mesh +/// \param GaussPoint_sie is specific internal energy at the GaussPoints on the mesh +/// \param GaussPoint_statev are the state variables at the GaussPoints on the mesh +/// \param num_mat_pts is the number of material points for mat_id +/// \param mat_id is material id +/// \param rk_num_bins is number of time integration storage bins +/// +///////////////////////////////////////////////////////////////////////////// +void init_press_sspd_stress(const Material_t& Materials, + const mesh_t& mesh, + const DCArrayKokkos& MaterialPoints_den, + const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const DCArrayKokkos& MaterialPoints_sspd, + const DCArrayKokkos& MaterialPoints_sie, + const DCArrayKokkos& MaterialPoints_statev, + const size_t rk_num_bins, + const size_t num_mat_pts, + const size_t mat_id){ + + + // ------- + // the call to the model initialization goes here + // ------- + + // --- pressure and sound speed --- + // loop over the material points + FOR_ALL(mat_point_lid, 0, num_mat_pts, { + + // --- Pressure --- + Materials.MaterialFunctions(mat_id).calc_pressure( + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(0, mat_point_lid), + Materials.eos_global_vars); + + // --- Sound Speed --- + Materials.MaterialFunctions(mat_id).calc_sound_speed( + MaterialPoints_pres, + MaterialPoints_stress, + mat_point_lid, + mat_id, + MaterialPoints_statev, + MaterialPoints_sspd, + MaterialPoints_den(mat_point_lid), + MaterialPoints_sie(0, mat_point_lid), + Materials.eos_global_vars); + }); // end pressure and sound speed + + + // --- stress tensor --- + for(size_t rk_level=0; rk_level& node_coords, + const DCArrayKokkos& node_mass, + const DCArrayKokkos& corner_mass, + const DCArrayKokkos& MaterialPoints_mass, + const DCArrayKokkos& MaterialToMeshMaps_elem, + const size_t num_mat_elems){ + + + FOR_ALL(mat_elem_lid, 0, num_mat_elems, { + + // get elem gid + size_t elem_gid = MaterialToMeshMaps_elem(mat_elem_lid); + + // calculate the fraction of matpt mass to scatter to each corner + double corner_frac = 1.0/((double)mesh.num_nodes_in_elem); // =1/8 + + // partion the mass to the corners + for(size_t corner_lid=0; corner_lid& node_coords, + const DCArrayKokkos& node_mass, + const DCArrayKokkos& corner_mass){ + + + FOR_ALL(node_gid, 0, mesh.num_nodes, { + for (size_t corner_lid = 0; corner_lid < mesh.num_corners_in_node(node_gid); corner_lid++) { + + size_t corner_gid = mesh.corners_in_node(node_gid, corner_lid); + + node_mass(node_gid) += corner_mass(corner_gid); + } // end for elem_lid + }); // end parallel loop over nodes in the mesh + +} // end function calculate node mass + + + + +#endif \ No newline at end of file diff --git a/single-node-refactor/src/common/simulation_parameters.h b/single-node-refactor/src/common/simulation_parameters.h index 0287ba764..f5330a5fe 100644 --- a/single-node-refactor/src/common/simulation_parameters.h +++ b/single-node-refactor/src/common/simulation_parameters.h @@ -63,7 +63,7 @@ struct SimulationParameters_t std::vector solver_inputs; ///< Solvers to use during the simulation CArrayKokkos region_fills; ///< Region data for simulation mesh, set the initial conditions - + CArray region_fills_host; ///< Region data on CPU, set the initial conditions }; // simulation_parameters_t diff --git a/single-node-refactor/src/common/state.h b/single-node-refactor/src/common/state.h index 39057b8ff..bdce7f43a 100644 --- a/single-node-refactor/src/common/state.h +++ b/single-node-refactor/src/common/state.h @@ -38,6 +38,12 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. using namespace mtr; + + + + + + ///////////////////////////////////////////////////////////////////////////// /// /// \struct node_t @@ -60,53 +66,57 @@ struct node_t }; // end method }; // end node_t + + ///////////////////////////////////////////////////////////////////////////// /// -/// \struct Zone_t +/// \struct GaussPoint_t /// -/// \brief Stores state information associated with zone index space +/// \brief Stores state information associated with the Gauss point /// ///////////////////////////////////////////////////////////////////////////// -struct Zone_t +struct GaussPoint_t { - DCArrayKokkos sie; ///< coefficients for the sie polynomial field + //const size_t num_bins = 3; - // initialization method (num_rk_storage_bins, num_nodes, num_dims) - void initialize(size_t num_rk, size_t num_nodes, size_t num_dims) + DCArrayKokkos vol; ///< GaussPoint volume + DCArrayKokkos div; ///< GaussPoint divergence of velocity + + + // initialization method (num_rk_storage_bins, num_cells, num_dims) + void initialize(size_t num_rk, size_t num_gauss_pnts, size_t num_dims) { - this->sie = DCArrayKokkos(num_rk, num_nodes, num_dims, "zone_sie"); - }; // end method + this->vol = DCArrayKokkos(num_gauss_pnts, "gauss_point_volume"); + this->div = DCArrayKokkos(num_gauss_pnts, "gauss_point_div"); -}; // end zone_t + }; // end method +}; // end GuassPoint_t ///////////////////////////////////////////////////////////////////////////// /// -/// \struct GaussPoint_t +/// \struct MaterialtoMeshMap_t /// -/// \brief Stores state information associated with the Gauss point +/// \brief Stores state information associated with maps from material to mesh maps /// ///////////////////////////////////////////////////////////////////////////// -struct GaussPoint_t +struct MaterialToMeshMap_t { - DCArrayKokkos vol; ///< GAussPoint volume - DCArrayKokkos div; ///< GAussPoint divergence of velocity - - DCArrayKokkos mat_id; ///< MaterialPoint material index + size_t num_material_elems; ///< returns the exact number of matpts - DCArrayKokkos eroded; ///< MaterialPoint eroded or not + DCArrayKokkos elem; ///< returns the elem for this material - // initialization method (num_rk_storage_bins, num_cells, num_dims) - void initialize(size_t num_rk, size_t num_elems, size_t num_dims) + // initialization method for FE-SGH and MPM methods (max number of elems needed) + void initialize(size_t num_elem_max) { - this->vol = DCArrayKokkos(num_elems, "gauss_point_volume"); - this->div = DCArrayKokkos(num_elems, "gauss_point_div"); - this->mat_id = DCArrayKokkos(num_elems, "gauss_point_mat_id"); - this->eroded = DCArrayKokkos(num_elems, "gauss_point_eroded"); + this->elem = DCArrayKokkos(num_elem_max, "material_pt_to_elem"); }; // end method -}; // end GuassPoint_t + +}; // end MaterialtoMeshMaps_t + + ///////////////////////////////////////////////////////////////////////////// @@ -118,27 +128,93 @@ struct GaussPoint_t ///////////////////////////////////////////////////////////////////////////// struct MaterialPoint_t { + size_t num_material_points; ///< the actual number of material points, omitting the buffer + DCArrayKokkos den; ///< MaterialPoint density DCArrayKokkos pres; ///< MaterialPoint pressure DCArrayKokkos stress; ///< MaterialPoint stress DCArrayKokkos sspd; ///< MaterialPoint sound speed - DCArrayKokkos sie; ///< MaterialPoint specific internal energy DCArrayKokkos mass; ///< MaterialPoint mass - DCArrayKokkos statev; ///< MaterialPoint state variable + + DCArrayKokkos eroded; ///< MaterialPoint eroded or not flag - // initialization method (num_rk_storage_bins, num_cells, num_dims) - void initialize(size_t num_rk, size_t num_elems, size_t num_dims) + DCArrayKokkos sie; ///< coefficients for the sie in strong form, only used in some methods e.g., FE-SGH and MPM + + // Material Models are stored on Material points + DCArrayKokkos statev; // a place holder to get things to compile + DCArrayKokkos eos_state_vars; ///< Array of state variables for the EOS + DCArrayKokkos strength_state_vars; ///< Array of state variables for the strength + + + // initialization method (num_rk_storage_bins, num_pts_max, num_dims) + void initialize(size_t num_rk, size_t num_pts_max, size_t num_dims) + { + this->den = DCArrayKokkos(num_pts_max, "material_point_density"); + this->pres = DCArrayKokkos(num_pts_max, "material_point_pressure"); + this->stress = DCArrayKokkos(num_rk, num_pts_max, num_dims, num_dims, "material_point_stress"); + this->sspd = DCArrayKokkos(num_pts_max, "material_point_sspd"); + this->mass = DCArrayKokkos(num_pts_max, "material_point_mass"); + this->sie = DCArrayKokkos(num_rk, num_pts_max, "material_point_sie"); + this->eroded = DCArrayKokkos(num_pts_max, "material_point_eroded"); + }; // end method + + // initialization method for arbitrary-order FE (num_rk_storage_bins, num_pts_max, num_dims) + void initialize_Pn(size_t num_rk, size_t num_pts_max, size_t num_dims) { - this->den = DCArrayKokkos(num_elems, "material_point_density"); - this->pres = DCArrayKokkos(num_elems, "material_point_pressure"); - this->stress = DCArrayKokkos(num_rk, num_elems, num_dims, num_dims, "material_point_stress"); - this->sspd = DCArrayKokkos(num_elems, "material_point_sspd"); - this->sie = DCArrayKokkos(num_rk, num_elems, "material_point_sie"); // only used with DG - this->mass = DCArrayKokkos(num_elems, "material_point_mass"); + this->den = DCArrayKokkos(num_pts_max, "material_point_density"); + this->pres = DCArrayKokkos(num_pts_max, "material_point_pressure"); + this->stress = DCArrayKokkos(num_rk, num_pts_max, num_dims, num_dims, "material_point_stress"); + this->sspd = DCArrayKokkos(num_pts_max, "material_point_sspd"); + this->mass = DCArrayKokkos(num_pts_max, "material_point_mass"); }; // end method }; // end MaterialPoint + +///////////////////////////////////////////////////////////////////////////// +/// +/// \struct MaterialZone_t +/// +/// \brief Stores state information associated with zone index space +/// +///////////////////////////////////////////////////////////////////////////// +struct MaterialZone_t +{ + size_t num_material_zones; ///< the actual number of material zones, omitting the buffer + + DCArrayKokkos sie; ///< coefficients for the sie polynomial field + + // initialization method for arbitrary-order FE (num_rk_storage_bins, num_zones) + void initialize_Pn(size_t num_rk, size_t num_zones_max) + { + this->sie = DCArrayKokkos(num_rk, num_zones_max, "material_zone_sie"); + }; // end method + +}; // end MaterialZone_t + + +///////////////////////////////////////////////////////////////////////////// +/// +/// \struct MaterialCorner_t +/// +/// \brief Stores state information associated with a material in an element corner +/// +///////////////////////////////////////////////////////////////////////////// +struct MaterialCorner_t +{ + size_t num_material_corners; ///< the actual number of material corners, omitting the buffer + + DCArrayKokkos force; ///< Corner force for the material + + // initialization method (num_corners, num_dims) + void initialize(size_t num_corners_max, size_t num_dims) + { + this->force = DCArrayKokkos(num_corners_max, num_dims, "material_corner_force"); + }; // end method +}; // end material corner + + + ///////////////////////////////////////////////////////////////////////////// /// /// \struct corner_t @@ -159,4 +235,187 @@ struct corner_t }; // end method }; // end corner_t + +///////////////////////////////////////////////////////////////////////////// +/// +/// \struct map for getting the corners in the material index +/// +/// \brief Stores state information associated with material corner index space +/// +///////////////////////////////////////////////////////////////////////////// +struct corners_in_mat_t +{ + private: + size_t num_corners_in_elem_; + public: + corners_in_mat_t() { + }; + + corners_in_mat_t(const size_t num_corners_in_elem_inp) { + this->num_corners_in_elem_ = num_corners_in_elem_inp; + }; + + // return global corner index for given local corner index in a material storage + size_t host(const size_t mat_storage_lid, const size_t corner_lid) const + { + return mat_storage_lid * num_corners_in_elem_ + corner_lid; + }; + + // Return the global corner ID given a material storage gloabl ID and a local corner ID + KOKKOS_INLINE_FUNCTION + size_t operator()(const size_t mat_storage_lid, const size_t corner_lid) const + { + return mat_storage_lid * num_corners_in_elem_ + corner_lid; + }; +}; + +///////////////////////////////////////////////////////////////////////////// +/// +/// \struct maps for high-order FE methods +/// +/// \brief Stores state information associated with other mesh index spaces +/// +///////////////////////////////////////////////////////////////////////////// + +struct zones_in_mat_t +{ + private: + size_t num_zones_in_elem_; + public: + zones_in_mat_t() { + }; + + zones_in_mat_t(const size_t num_zones_in_elem_inp) { + this->num_zones_in_elem_ = num_zones_in_elem_inp; + }; + + // return global zone index for given local zone index in a material storage + size_t host(const size_t mat_storage_lid, const size_t zone_lid) const + { + return mat_storage_lid * num_zones_in_elem_ + zone_lid; + }; + + // Return the global zone ID given a material storage gloabl ID and a local zone ID + KOKKOS_INLINE_FUNCTION + size_t operator()(const size_t mat_storage_lid, const size_t zone_lid) const + { + return mat_storage_lid * num_zones_in_elem_ + zone_lid; + }; +}; + +// if material points are defined strictly internal to the element. +struct legendre_in_mat_t +{ + private: + size_t num_leg_gauss_in_elem_; + public: + legendre_in_mat_t() { + }; + + legendre_in_mat_t(const size_t num_leg_gauss_in_elem_inp) { + this->num_leg_gauss_in_elem_ = num_leg_gauss_in_elem_inp; + }; + + // return global gauss index for given local gauss index in a material storage + size_t host(const size_t mat_storage_lid, const size_t leg_gauss_lid) const + { + return mat_storage_lid * num_leg_gauss_in_elem_ + leg_gauss_lid; + }; + + // Return the global gauss ID given a material storage gloabl ID and a local gauss ID + KOKKOS_INLINE_FUNCTION + size_t operator()(const size_t mat_storage_lid, const size_t leg_gauss_lid) const + { + return mat_storage_lid * num_leg_gauss_in_elem_ + leg_gauss_lid; + }; +}; + +/// if material points are defined at element interfaces, e.g., for nodal DG +struct lobatto_in_mat_t +{ + private: + size_t num_lob_gauss_in_elem_; + public: + lobatto_in_mat_t() { + }; + + lobatto_in_mat_t(const size_t num_lob_gauss_in_elem_inp) { + this->num_lob_gauss_in_elem_ = num_lob_gauss_in_elem_inp; + }; + + // return global gauss index for given local gauss index in a material storage + size_t host(const size_t mat_storage_lid, const size_t lob_gauss_lid) const + { + return mat_storage_lid * num_lob_gauss_in_elem_ + lob_gauss_lid; + }; + + // Return the global gauss ID given a material storage ID and a local gauss ID + KOKKOS_INLINE_FUNCTION + size_t operator()(const size_t mat_storage_lid, const size_t lob_gauss_lid) const + { + return mat_storage_lid * num_lob_gauss_in_elem_ + lob_gauss_lid; + }; +}; + + +// the local id for material points in elem +struct points_in_mat_t +{ + private: + size_t num_points_in_elem_; + public: + points_in_mat_t() { + }; + + points_in_mat_t(const size_t num_points_in_elem_inp) { + this->num_points_in_elem_ = num_points_in_elem_inp; + }; + + // return global gauss index for given local gauss index in a material storage + size_t host(const size_t mat_storage_lid, const size_t points_lid) const + { + return mat_storage_lid * num_points_in_elem_ + points_lid; + }; + + // Return the global gauss ID given a material storage gloabl ID and a local gauss ID + KOKKOS_INLINE_FUNCTION + size_t operator()(const size_t mat_storage_lid, const size_t points_lid) const + { + return mat_storage_lid * num_points_in_elem_ + points_lid; + }; +}; + +///////////////////////////////////////////////////////////////////////////// +/// +/// \struct state_t +/// +/// \brief Stores all state +/// +///////////////////////////////////////////////////////////////////////////// +struct State_t +{ + // --------------------------------------------------------------------- + // state data on mesh declarations + // --------------------------------------------------------------------- + node_t node; + GaussPoint_t GaussPoints; + corner_t corner; + + // --------------------------------------------------------------------- + // material to mesh maps + // --------------------------------------------------------------------- + CArray MaterialToMeshMaps; ///< access as MaterialToMeshMaps(mat_id).elem(mat_storage_lid) + corners_in_mat_t corners_in_mat_elem; ///< access the corner mat lid using (mat_elem_lid, corn_lid) + points_in_mat_t points_in_mat_elem; ///< for accessing e.g., guass points mat lid with arbitrary-order FE + zones_in_mat_t zones_in_mat_elem; ///< for accessing sub-zones mat lid with arbitrary-order FE + + // --------------------------------------------------------------------- + // material state, compressed, and sequentially accessed + // --------------------------------------------------------------------- + CArray MaterialPoints; ///< access as MaterialPoints(mat_id).var(mat_pt) + CArray MaterialCorners; ///< access as MaterialCorners(mat_id).var(mat_corner), not used with MPM + CArray MaterialZones; ///< access as MaterialZones(mat_id).var(mat_zone), only used with arbitrary-order FE + +}; // end state_t + #endif diff --git a/single-node-refactor/src/driver.h b/single-node-refactor/src/driver.h index 6afe4115f..b3a38a219 100644 --- a/single-node-refactor/src/driver.h +++ b/single-node-refactor/src/driver.h @@ -44,70 +44,6 @@ #include "state.h" -void fill_regions(CArrayKokkos&, - CArray &, - Material_t&, - mesh_t&, - node_t&, - MaterialPoint_t&, - GaussPoint_t&, - corner_t&); - -// ============================================================================== -// Functions to read voxel mesh -// ============================================================================== -void user_voxel_init(DCArrayKokkos& elem_values, - double& dx, - double& dy, - double& dz, - double& orig_x, - double& orig_y, - double& orig_z, - size_t& voxel_num_i, - size_t& voxel_num_j, - size_t& voxel_num_k, - double scale_x, - double scale_y, - double scale_z, - std::string mesh_file); - -///////////////////////////////////////////////////////////////////////////// -/// -/// \fn get_id_device -/// -/// \brief This gives the index value of the point or the elem -/// -/// Assumes that the grid has an i,j,k structure -/// the elem = i + (j)*(num_points_i-1) + (k)*(num_points_i-1)*(num_points_j-1) -/// the point = i + (j)*num_points_i + (k)*num_points_i*num_points_j -/// -/// \param i index -/// \param j index -/// \param k index -/// \param Number of i indices -/// \param Number of j indices -/// -///////////////////////////////////////////////////////////////////////////// -KOKKOS_FUNCTION -int get_id_device(int i, int j, int k, int num_i, int num_j) -{ - return i + j * num_i + k * num_i * num_j; -} - -// for string delimiter parsing -std::vector split(std::string s, std::string delimiter); - -// retrieves multiple values between [ ] -std::vector extract_list(std::string str); - -const std::string WHITESPACE = " "; - -std::string ltrim(const std::string& s); - -std::string rtrim(const std::string& s); - -std::string trim(const std::string& s); - class Driver @@ -141,12 +77,9 @@ class Driver mesh_t mesh; // --------------------------------------------------------------------- - // state data type declarations + // state data type declaration // --------------------------------------------------------------------- - node_t node; - GaussPoint_t GaussPoints; - MaterialPoint_t MaterialPoints; - corner_t corner; + State_t State; int num_solvers = 0; @@ -183,19 +116,17 @@ class Driver std::cout << "Mesh file path: " << SimulationParamaters.mesh_input.file_path << std::endl; mesh_reader.set_mesh_file(SimulationParamaters.mesh_input.file_path.data()); mesh_reader.read_mesh(mesh, - MaterialPoints, - GaussPoints, - node, - corner, + State.GaussPoints, + State.node, + State.corner, num_dims, SimulationParamaters.dynamic_options.rk_num_bins); } else if (SimulationParamaters.mesh_input.source == mesh_input::generate) { mesh_builder.build_mesh(mesh, - MaterialPoints, - GaussPoints, - node, - corner, + State.GaussPoints, + State.node, + State.corner, SimulationParamaters); } else{ @@ -209,44 +140,29 @@ class Driver // --- calculate bdy sets ---// mesh.init_bdy_sets(num_bcs); - tag_bdys(BoundaryConditions, mesh, node.coords); + tag_bdys(BoundaryConditions, mesh, State.node.coords); mesh.build_boundry_node_sets(mesh); // Calculate element volume - geometry::get_vol(GaussPoints.vol, node.coords, mesh); + geometry::get_vol(State.GaussPoints.vol, State.node.coords, mesh); - // Create memory for state variables - //size_t max_num_vars = 0; - //size_t num_materials = Materials.num_eos_global_vars.size(); - //for (size_t mat_id=0; mat_id(mesh.num_elems, max_num_vars); // WARNING: HACK - - // --- apply the fill instructions over the Elements---// - Kokkos::fence(); - - //fill_regions(); - fill_regions(SimulationParamaters.region_fills, - SimulationParamaters.region_fills_host, - Materials, - mesh, - node, - MaterialPoints, - GaussPoints, - corner); - - - // --- Move the following sovler setup to yaml parsing routine - // Create solvers for (int solver_id = 0; solver_id < SimulationParamaters.solver_inputs.size(); solver_id++) { + if (SimulationParamaters.solver_inputs[solver_id].method == solver_input::SGH) { - SGH* sgh_solver = new SGH(); // , mesh, node, MaterialPoints, corner - sgh_solver->initialize(SimulationParamaters, Materials, BoundaryConditions); + + SGH* sgh_solver = new SGH(); + + sgh_solver->initialize(SimulationParamaters, + Materials, + mesh, + BoundaryConditions, + State); + solvers.push_back(sgh_solver); - } - } + } // end if SGH solver + + } // end for loop over solvers } // end initialize @@ -260,17 +176,17 @@ class Driver void setup() { std::cout << "Inside driver setup" << std::endl; + + // allocate state, setup models, and apply fill instructions for (auto& solver : solvers) { solver->setup(SimulationParamaters, Materials, - BoundaryConditions, mesh, - node, - MaterialPoints, - GaussPoints, - corner); - } - } + BoundaryConditions, + State); + } // end for over solvers + + } // end setup function of driver ///////////////////////////////////////////////////////////////////////////// /// @@ -287,11 +203,8 @@ class Driver Materials, BoundaryConditions, mesh, - node, - MaterialPoints, - GaussPoints, - corner); - } + State); + } // loop over solvers } ///////////////////////////////////////////////////////////////////////////// @@ -323,835 +236,3 @@ class Driver }; // end driver class - ///////////////////////////////////////////////////////////////////////////// - /// - /// \fn fill_regions - /// - /// \brief Fills mesh regions based on YAML input - /// - ///////////////////////////////////////////////////////////////////////////// - void fill_regions(CArrayKokkos& region_fills, - CArray ®ion_fills_host, - Material_t& Materials, - mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) - { - int num_fills = region_fills.size(); - printf("Num Fills's = %d\n", num_fills); - - // --------------------------------------------- - // variables from a voxel file - // --------------------------------------------- - DCArrayKokkos voxel_elem_mat_id; // 1 or 0 if material exist, or it is the material_id - double voxel_dx, voxel_dy, voxel_dz; // voxel mesh resolution, set by input file - double orig_x, orig_y, orig_z; // origin of voxel elem center mesh, set by input file - size_t voxel_num_i, voxel_num_j, voxel_num_k; // num voxel elements in each direction, set by input file - - DCArrayKokkos read_voxel_file(num_fills); // check to see if readVoxelFile - FOR_ALL(f_id, 0, num_fills, { - if (region_fills(f_id).volume == region::readVoxelFile) - { - read_voxel_file(f_id) = region::readVoxelFile; // read the voxel file - } - // add other mesh voxel files - else - { - read_voxel_file(f_id) = 0; - } - }); // end parallel for - read_voxel_file.update_host(); // copy to CPU if code is to read a file - Kokkos::fence(); - // --------------------------------------------- - - - - for (int f_id = 0; f_id < num_fills; f_id++) { - - // ---- - // voxel mesh setup - if (read_voxel_file.host(f_id) == region::readVoxelFile) - { - // read voxel mesh to get the values in the fcn interface - user_voxel_init(voxel_elem_mat_id, - voxel_dx, voxel_dy, voxel_dz, - orig_x, orig_y, orig_z, - voxel_num_i, voxel_num_j, voxel_num_k, - region_fills_host(f_id).scale_x, - region_fills_host(f_id).scale_y, - region_fills_host(f_id).scale_z, - region_fills_host(f_id).file_path); - - // copy values read from file to device - voxel_elem_mat_id.update_device(); - } // endif - // add else if for other mesh reads including STL-2-voxel - - - - int num_elems = mesh.num_elems; - // parallel loop over elements in mesh - FOR_ALL(elem_gid, 0, num_elems, { - for (int rk_level = 0; rk_level < 2; rk_level++) { - - // Set erosion flag to false - GaussPoints.eroded(elem_gid) = false; - - // calculate the coordinates and radius of the element - double elem_coords[3]; // note:initialization with a list won't work - elem_coords[0] = 0.0; - elem_coords[1] = 0.0; - elem_coords[2] = 0.0; - - // get the coordinates of the element center - for (int node_lid = 0; node_lid < mesh.num_nodes_in_elem; node_lid++) { - elem_coords[0] += node.coords(rk_level, mesh.nodes_in_elem(elem_gid, node_lid), 0); - elem_coords[1] += node.coords(rk_level, mesh.nodes_in_elem(elem_gid, node_lid), 1); - if (mesh.num_dims == 3) { - elem_coords[2] += node.coords(rk_level, mesh.nodes_in_elem(elem_gid, node_lid), 2); - } - else{ - elem_coords[2] = 0.0; - } - } // end loop over nodes in element - elem_coords[0] = (elem_coords[0] / mesh.num_nodes_in_elem); - elem_coords[1] = (elem_coords[1] / mesh.num_nodes_in_elem); - elem_coords[2] = (elem_coords[2] / mesh.num_nodes_in_elem); - - // for shapes with an origin (e.g., sphere and circle), accounting for the origin - double dist_x = elem_coords[0] - region_fills(f_id).origin[0]; - double dist_y = elem_coords[1] - region_fills(f_id).origin[1]; - double dist_z = elem_coords[2] - region_fills(f_id).origin[2]; - - // spherical radius - double radius = sqrt(dist_x * dist_x + - dist_y * dist_y + - dist_z * dist_z); - - // cylindrical radius - double radius_cyl = sqrt(dist_x * dist_x + - dist_y * dist_y); - - // default is not to fill the element - size_t fill_this = 0; - - // check to see if this element should be filled - switch (region_fills(f_id).volume) { - case region::global: - { - fill_this = 1; - break; - } - case region::box: - { - - double x_lower_bound = region_fills(f_id).x1; - double x_upper_bound = region_fills(f_id).x2; - - double y_lower_bound = region_fills(f_id).y1; - double y_upper_bound = region_fills(f_id).y2; - - double z_lower_bound = region_fills(f_id).z1; - double z_upper_bound = region_fills(f_id).z2; - - - if (elem_coords[0] >= x_lower_bound && elem_coords[0] <= x_upper_bound && - elem_coords[1] >= y_lower_bound && elem_coords[1] <= y_upper_bound && - elem_coords[2] >= z_lower_bound && elem_coords[2] <= z_upper_bound) { - fill_this = 1; - } - break; - } - case region::cylinder: - { - if (radius_cyl >= region_fills(f_id).radius1 - && radius_cyl <= region_fills(f_id).radius2) { - fill_this = 1; - } - break; - } - case region::sphere: - { - if (radius >= region_fills(f_id).radius1 - && radius <= region_fills(f_id).radius2) { - fill_this = 1; - } - break; - } - - case region::readVoxelFile: - { - - fill_this = 0; // default is no, don't fill it - - // find the closest element in the voxel mesh to this element - double i0_real = (elem_coords[0] - orig_x - region_fills(f_id).origin[0]) / (voxel_dx); - double j0_real = (elem_coords[1] - orig_y - region_fills(f_id).origin[1]) / (voxel_dy); - double k0_real = (elem_coords[2] - orig_z - region_fills(f_id).origin[2]) / (voxel_dz); - - int i0 = (int)i0_real; - int j0 = (int)j0_real; - int k0 = (int)k0_real; - - // look for the closest element in the voxel mesh - int elem_id0 = get_id_device(i0, j0, k0, voxel_num_i, voxel_num_j); - - // if voxel mesh overlaps this mesh, then fill it if =1 - if (elem_id0 < voxel_elem_mat_id.size() && elem_id0 >= 0 && - i0 >= 0 && j0 >= 0 && k0 >= 0 && - i0 < voxel_num_i && j0 < voxel_num_j && k0 < voxel_num_k) { - // voxel mesh elem values = 0 or 1 - fill_this = voxel_elem_mat_id(elem_id0); // values from file - - } // end if - - break; - - } // end case - case region::no_volume: - { - fill_this = 0; // default is no, don't fill it - - break; - } - default: - { - fill_this = 0; // default is no, don't fill it - - break; - } - - } // end of switch - - // paint the material state on the element - if (fill_this == 1) { - // density - MaterialPoints.den(elem_gid) = region_fills(f_id).den; - - // mass - MaterialPoints.mass(elem_gid) = MaterialPoints.den(elem_gid) * GaussPoints.vol(elem_gid); - - // specific internal energy - MaterialPoints.sie(rk_level, elem_gid) = region_fills(f_id).sie; - - GaussPoints.mat_id(elem_gid) = region_fills(f_id).material_id; - - size_t mat_id = GaussPoints.mat_id(elem_gid); // short name - - // get state_vars from the input file or read them in - if (false) { // Materials.MaterialEnums(mat_id).strength_setup == model_init::user_init) { - // use the values read from a file to get elem state vars - // for (size_t var = 0; var < Materials.num_eos_state_vars(mat_id); var++) { - // MaterialPoints.statev(elem_gid, var) = file_state_vars(mat_id, elem_gid, var); - // } // end for - } - else{ - // use the values in the input file - // set state vars for the region where mat_id resides - - //int num_eos_global_vars = Materials.eos_global_vars.stride(mat_id); // ragged-right storage - - - //for (size_t var = 0; var < num_eos_global_vars; var++) { - // MaterialPoints.statev(elem_gid, var) = Materials.eos_global_vars(mat_id,var); // state_vars(mat_id, var); // WARNING: HACK - //} // end for - - - } // end logical on type - - // --- stress tensor --- - // always 3D even for 2D-RZ - for (size_t i = 0; i < 3; i++) { - for (size_t j = 0; j < 3; j++) { - MaterialPoints.stress(rk_level, elem_gid, i, j) = 0.0; - } - } // end for - - // --- Calculate Pressure and Stress --- - - // --- Pressure --- - Materials.MaterialFunctions(mat_id).calc_pressure( - MaterialPoints.pres, - MaterialPoints.stress, - elem_gid, - GaussPoints.mat_id(elem_gid), - MaterialPoints.statev, - MaterialPoints.sspd, - MaterialPoints.den(elem_gid), - MaterialPoints.sie(rk_level, elem_gid), - Materials.eos_global_vars); - - // --- Sound Speed --- - Materials.MaterialFunctions(mat_id).calc_sound_speed( - MaterialPoints.pres, - MaterialPoints.stress, - elem_gid, - GaussPoints.mat_id(elem_gid), - MaterialPoints.statev, - MaterialPoints.sspd, - MaterialPoints.den(elem_gid), - MaterialPoints.sie(rk_level, elem_gid), - Materials.eos_global_vars); - - // loop over the nodes of this element and apply velocity - for (size_t node_lid = 0; node_lid < mesh.num_nodes_in_elem; node_lid++) { - // get the mesh node index - size_t node_gid = mesh.nodes_in_elem(elem_gid, node_lid); - - // --- Velocity --- - switch (region_fills(f_id).velocity) { - case init_conds::cartesian: - { - node.vel(rk_level, node_gid, 0) = region_fills(f_id).u; - node.vel(rk_level, node_gid, 1) = region_fills(f_id).v; - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = region_fills(f_id).w; - } - - break; - } - // radial in the (x,y) plane where x=r*cos(theta) and y=r*sin(theta) - case init_conds::radial: - { - // Setting up radial - double dir[2]; - dir[0] = 0.0; - dir[1] = 0.0; - double radius_val = 0.0; - - for (int dim = 0; dim < 2; dim++) { - dir[dim] = node.coords(rk_level, node_gid, dim); - radius_val += node.coords(rk_level, node_gid, dim) * node.coords(rk_level, node_gid, dim); - } // end for - radius_val = sqrt(radius_val); - - for (int dim = 0; dim < 2; dim++) { - if (radius_val > 1.0e-14) { - dir[dim] /= (radius_val); - } - else{ - dir[dim] = 0.0; - } - } // end for - - node.vel(rk_level, node_gid, 0) = region_fills(f_id).speed * dir[0]; - node.vel(rk_level, node_gid, 1) = region_fills(f_id).speed * dir[1]; - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = 0.0; - } - - break; - } - case init_conds::spherical: - { - // Setting up spherical - double dir[3]; - dir[0] = 0.0; - dir[1] = 0.0; - dir[2] = 0.0; - double radius_val = 0.0; - - for (int dim = 0; dim < 3; dim++) { - dir[dim] = node.coords(rk_level, node_gid, dim); - radius_val += node.coords(rk_level, node_gid, dim) * node.coords(rk_level, node_gid, dim); - } // end for - radius_val = sqrt(radius_val); - - for (int dim = 0; dim < 3; dim++) { - if (radius_val > 1.0e-14) { - dir[dim] /= (radius_val); - } - else{ - dir[dim] = 0.0; - } - } // end for - - node.vel(rk_level, node_gid, 0) = region_fills(f_id).speed * dir[0]; - node.vel(rk_level, node_gid, 1) = region_fills(f_id).speed * dir[1]; - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = region_fills(f_id).speed * dir[2]; - } - - break; - } - case init_conds::radial_linear: - { - printf("**** Radial_linear initial conditions not yet supported ****\n"); - break; - } - case init_conds::spherical_linear: - { - printf("**** spherical_linear initial conditions not yet supported ****\n"); - break; - } - case init_conds::tg_vortex: - { - node.vel(rk_level, node_gid, 0) = sin(PI * node.coords(rk_level, node_gid, 0)) * cos(PI * node.coords(rk_level, node_gid, 1)); - node.vel(rk_level, node_gid, 1) = -1.0 * cos(PI * node.coords(rk_level, node_gid, 0)) * sin(PI * node.coords(rk_level, node_gid, 1)); - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = 0.0; - } - - break; - } - - case init_conds::no_ic_vel: - { - // no velocity - node.vel(rk_level, node_gid, 0) = 0.0; - node.vel(rk_level, node_gid, 1) = 0.0; - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = 0.0; - } - - break; - } - default: - { - // no velocity - node.vel(rk_level, node_gid, 0) = 0.0; - node.vel(rk_level, node_gid, 1) = 0.0; - if (mesh.num_dims == 3) { - node.vel(rk_level, node_gid, 2) = 0.0; - } - - break; - } - } // end of switch - } // end loop over nodes of element - - if (region_fills(f_id).velocity == init_conds::tg_vortex) { - MaterialPoints.pres(elem_gid) = 0.25 * (cos(2.0 * PI * elem_coords[0]) + cos(2.0 * PI * elem_coords[1]) ) + 1.0; - - // p = rho*ie*(gamma - 1) - double gamma = Materials.eos_global_vars(mat_id,0); // makes sure it matches the gamma in the gamma law function - MaterialPoints.sie(rk_level, elem_gid) = - MaterialPoints.pres(elem_gid) / (region_fills(f_id).den * (gamma - 1.0)); - } // end if - - } // end if fill - } // end RK loop - }); // end FOR_ALL element loop - Kokkos::fence(); - } // end for loop over fills - - // // calculate the corner massess if 2D - // if (mesh.num_dims == 2) { - // FOR_ALL(elem_gid, 0, mesh.num_elems, { - // // facial area of the corners - // double corner_areas_array[4]; - - // ViewCArrayKokkos corner_areas(&corner_areas_array[0], 4); - // ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 4); - - // geometry::get_area_weights2D(corner_areas, elem_gid, node_coords, elem_node_gids); - - // // loop over the corners of the element and calculate the mass - // for (size_t corner_lid = 0; corner_lid < 4; corner_lid++) { - // size_t corner_gid = mesh.corners_in_elem(elem_gid, corner_lid); - // corner_mass(corner_gid) = corner_areas(corner_lid) * MaterialPoints.den(elem_gid); // node radius is added later - // } // end for over corners - // }); - // } // end of - - // calculate the nodal mass - FOR_ALL(node_gid, 0, mesh.num_nodes, { - node.mass(node_gid) = 0.0; - - if (mesh.num_dims == 3) { - for (size_t elem_lid = 0; elem_lid < mesh.num_corners_in_node(node_gid); elem_lid++) { - size_t elem_gid = mesh.elems_in_node(node_gid, elem_lid); - node.mass(node_gid) += 1.0 / 8.0 * MaterialPoints.mass(elem_gid); - } // end for elem_lid - } // end if dims=3 - else{ - // 2D-RZ - for (size_t corner_lid = 0; corner_lid < mesh.num_corners_in_node(node_gid); corner_lid++) { - size_t corner_gid = mesh.corners_in_node(node_gid, corner_lid); - node.mass(node_gid) += corner.mass(corner_gid); // sans the radius so it is areal node mass - - corner.mass(corner_gid) *= node.coords(1, node_gid, 1); // true corner mass now - } // end for elem_lid - } // end else - }); // end FOR_ALL - Kokkos::fence(); - - } // end fill regions - - -// ----------------------------------------------------------------------------- -// The function to read a voxel vtk file from Dream3d and intialize the mesh -// ------------------------------------------------------------------------------ -void user_voxel_init(DCArrayKokkos& elem_values, - double& dx, - double& dy, - double& dz, - double& orig_x, - double& orig_y, - double& orig_z, - size_t& num_elems_i, - size_t& num_elems_j, - size_t& num_elems_k, - double scale_x, - double scale_y, - double scale_z, - std::string mesh_file) -{ - std::string MESH = mesh_file; // user specified - - std::ifstream in; // FILE *in; - in.open(MESH); - - // check to see of a mesh was supplied when running the code - if (in) - { - printf("\nReading the 3D voxel mesh: "); - std::cout << MESH << std::endl; - } - else - { - std::cout << "\n\n**********************************\n\n"; - std::cout << " ERROR:\n"; - std::cout << " Voxel vtk input does not exist \n"; - std::cout << "**********************************\n\n" << std::endl; - std::exit(EXIT_FAILURE); - } // end if - - size_t i; // used for writing information to file - size_t point_id; // the global id for the point - size_t elem_id; // the global id for the elem - size_t this_point; // a local id for a point in a elem (0:7 for a Hexahedral elem) - - size_t num_points_i; - size_t num_points_j; - size_t num_points_k; - - size_t num_dims = 3; - - std::string token; - - bool found = false; - - // look for POINTS - i = 0; - while (found == false) { - std::string str; - std::string delimiter = " "; - std::getline(in, str); - std::vector v = split(str, delimiter); - - // looking for the following text: - // POINTS %d float - if (v[0] == "DIMENSIONS") - { - num_points_i = std::stoi(v[1]); - num_points_j = std::stoi(v[2]); - num_points_k = std::stoi(v[3]); - printf("Num voxel nodes read in = %zu, %zu, %zu\n", num_points_i, num_points_j, num_points_k); - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find POINTS \n"); - break; - } // end if - - i++; - } // end while - - found = false; - - int num_points = num_points_i * num_points_j * num_points_k; - CArray pt_coords_x(num_points_i); - CArray pt_coords_y(num_points_j); - CArray pt_coords_z(num_points_k); - - while (found == false) { - std::string str; - std::string str0; - std::string delimiter = " "; - std::getline(in, str); - std::vector v = split(str, delimiter); - - // looking for the following text: - if (v[0] == "X_COORDINATES") - { - size_t num_saved = 0; - - while (num_saved < num_points_i - 1) { - // get next line - std::getline(in, str0); - - // remove starting and trailing spaces - str = trim(str0); - std::vector v_coords = split(str, delimiter); - - // loop over the contents of the vector v_coords - for (size_t this_point = 0; this_point < v_coords.size(); this_point++) - { - pt_coords_x(num_saved) = scale_x*std::stod(v_coords[this_point]); - num_saved++; - } // end for - } // end while - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find X_COORDINATES \n"); - break; - } // end if - - i++; - } // end while - found = false; - - while (found == false) { - std::string str; - std::string str0; - std::string delimiter = " "; - std::getline(in, str); - std::vector v = split(str, delimiter); - - // looking for the following text: - if (v[0] == "Y_COORDINATES") - { - size_t num_saved = 0; - - while (num_saved < num_points_j - 1) { - // get next line - std::getline(in, str0); - - // remove starting and trailing spaces - str = trim(str0); - std::vector v_coords = split(str, delimiter); - - // loop over the contents of the vector v_coords - for (size_t this_point = 0; this_point < v_coords.size(); this_point++) - { - pt_coords_y(num_saved) = scale_y*std::stod(v_coords[this_point]); - num_saved++; - } // end for - } // end while - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find Y_COORDINATES \n"); - break; - } // end if - - i++; - } // end while - found = false; - - while (found == false) { - std::string str; - std::string str0; - std::string delimiter = " "; - std::getline(in, str); - std::vector v = split(str, delimiter); - - // looking for the following text: - if (v[0] == "Z_COORDINATES") - { - size_t num_saved = 0; - - while (num_saved < num_points_k - 1) { - // get next line - std::getline(in, str0); - - // remove starting and trailing spaces - str = trim(str0); - std::vector v_coords = split(str, delimiter); - - // loop over the contents of the vector v_coords - for (size_t this_point = 0; this_point < v_coords.size(); this_point++) - { - pt_coords_z(num_saved) = scale_z*std::stod(v_coords[this_point]); - num_saved++; - } // end for - } // end while - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find Z_COORDINATES \n"); - break; - } // end if - - i++; - } // end while - found = false; - - size_t num_elems; - num_elems_i = num_points_i - 1; - num_elems_j = num_points_j - 1; - num_elems_k = num_points_k - 1; - - // center to center distance between first and last elem along each edge - double Lx = (pt_coords_x(num_points_i - 2) - pt_coords_x(0)); - double Ly = (pt_coords_y(num_points_j - 2) - pt_coords_y(0)); - double Lz = (pt_coords_z(num_points_k - 2) - pt_coords_z(0)); - - // spacing between elems - dx = Lx / ((double) num_elems_i); - dy = Ly / ((double) num_elems_j); - dz = Lz / ((double) num_elems_k); - - // element mesh origin - orig_x = 0.5 * (pt_coords_x(0) + pt_coords_x(1)), - orig_y = 0.5 * (pt_coords_y(0) + pt_coords_y(1)), - orig_z = 0.5 * (pt_coords_z(0) + pt_coords_z(1)), - - // look for CELLS - i = 0; - while (found == false) { - std::string str; - std::getline(in, str); - - std::string delimiter = " "; - std::vector v = split(str, delimiter); - - // looking for the following text: - // CELLS num_elems size - if (v[0] == "CELL_DATA") - { - num_elems = std::stoi(v[1]); - printf("Num voxel elements read in %zu\n", num_elems); - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find CELL_DATA \n"); - break; - } // end if - - i++; - } // end while - found = false; - - // allocate memory for element voxel values - elem_values = DCArrayKokkos(num_elems); - - // reading the cell data - while (found == false) { - std::string str; - std::string str0; - - std::string delimiter = " "; - std::getline(in, str); - std::vector v = split(str, delimiter); - - // looking for the following text: - if (v[0] == "LOOKUP_TABLE") - { - size_t num_saved = 0; - - while (num_saved < num_elems - 1) { - // get next line - std::getline(in, str0); - - // remove starting and trailing spaces - str = trim(str0); - std::vector v_values = split(str, delimiter); - - // loop over the contents of the vector v_coords - for (size_t this_elem = 0; this_elem < v_values.size(); this_elem++) - { - // save integers (0 or 1) to host side - elem_values.host(num_saved) = std::stoi(v_values[this_elem]); - num_saved++; - } // end for - - // printf(" done with one row of data \n"); - } // end while - - found = true; - } // end if - - if (i > 1000) - { - printf("ERROR: Failed to find LOOKUP_TABLE data \n"); - break; - } // end if - - i++; - } // end while - found = false; - - printf("\n"); - - in.close(); -} // end routine - - -// Code from stackover flow for string delimiter parsing -std::vector split(std::string s, std::string delimiter) -{ - size_t pos_start = 0, pos_end, delim_len = delimiter.length(); - std::string token; - std::vector res; - - while ((pos_end = s.find(delimiter, pos_start)) != std::string::npos) { - token = s.substr(pos_start, pos_end - pos_start); - pos_start = pos_end + delim_len; - res.push_back(token); - } - - res.push_back(s.substr(pos_start)); - return res; -} // end of split - - -// retrieves multiple values between [ ] -std::vector extract_list(std::string str) -{ - // replace '[' with a space and ']' with a space - std::replace(str.begin(), str.end(), '[', ' '); - std::replace(str.begin(), str.end(), ']', ' '); - - std::vector str_values; - std::vector values; - - // exact the str values into a vector - str_values = split(str, ","); - - // convert the text values into double values - for (auto& word : str_values) - { - values.push_back(atof(word.c_str()) ); - } // end for - - return values; -} // end of extract_list - - -std::string ltrim(const std::string& s) -{ - size_t start = s.find_first_not_of(WHITESPACE); - return (start == std::string::npos) ? "" : s.substr(start); -} - - -std::string rtrim(const std::string& s) -{ - size_t end = s.find_last_not_of(WHITESPACE); - return (end == std::string::npos) ? "" : s.substr(0, end + 1); -} - -std::string trim(const std::string& s) -{ - return rtrim(ltrim(s)); -} \ No newline at end of file diff --git a/single-node-refactor/src/input/CMakeLists.txt b/single-node-refactor/src/input/CMakeLists.txt index 560e4d741..a3015a0e3 100644 --- a/single-node-refactor/src/input/CMakeLists.txt +++ b/single-node-refactor/src/input/CMakeLists.txt @@ -18,16 +18,21 @@ endif() include_directories(../common) -set(SRC_Files -parse_yaml.cpp -Yaml.cpp) +set(YAML_SRC_Files +${CMAKE_CURRENT_SOURCE_DIR}/parse_yaml.cpp +${CMAKE_CURRENT_SOURCE_DIR}/parse_yaml.h +${CMAKE_CURRENT_SOURCE_DIR}/Yaml.cpp +${CMAKE_CURRENT_SOURCE_DIR}/Yaml.hpp +PARENT_SCOPE +) # set(SGH_Solver_SRC src/sgh_solver.cpp ) message("\n ****** ADDING PARSE_YAML LIBRARY ******** \n ") -add_library(parse_yaml ${SRC_Files} ) -target_link_libraries(parse_yaml matar Kokkos::kokkos) + +# add_library(parse_yaml ${SRC_Files} ) +# target_link_libraries(parse_yaml matar Kokkos::kokkos) # target_include_directories(sgh_solver PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include) diff --git a/single-node-refactor/src/input/parse_yaml.cpp b/single-node-refactor/src/input/parse_yaml.cpp index 447886306..fb57cb857 100644 --- a/single-node-refactor/src/input/parse_yaml.cpp +++ b/single-node-refactor/src/input/parse_yaml.cpp @@ -105,6 +105,65 @@ std::vector exact_array_values(std::string s, std::string delimiter return res; } // end of extract_array_values +// Code from stackover flow for string delimiter parsing +std::vector split(std::string s, std::string delimiter) +{ + size_t pos_start = 0, pos_end, delim_len = delimiter.length(); + std::string token; + std::vector res; + + while ((pos_end = s.find(delimiter, pos_start)) != std::string::npos) { + token = s.substr(pos_start, pos_end - pos_start); + pos_start = pos_end + delim_len; + res.push_back(token); + } + + res.push_back(s.substr(pos_start)); + return res; +} // end of split + + +// retrieves multiple values between [ ] +std::vector extract_list(std::string str) +{ + // replace '[' with a space and ']' with a space + std::replace(str.begin(), str.end(), '[', ' '); + std::replace(str.begin(), str.end(), ']', ' '); + + std::vector str_values; + std::vector values; + + // exact the str values into a vector + str_values = split(str, ","); + + // convert the text values into double values + for (auto& word : str_values) + { + values.push_back(atof(word.c_str()) ); + } // end for + + return values; +} // end of extract_list + + +std::string ltrim(const std::string& s) +{ + size_t start = s.find_first_not_of(WHITESPACE); + return (start == std::string::npos) ? "" : s.substr(start); +} + + +std::string rtrim(const std::string& s) +{ + size_t end = s.find_last_not_of(WHITESPACE); + return (end == std::string::npos) ? "" : s.substr(0, end + 1); +} + +std::string trim(const std::string& s) +{ + return rtrim(ltrim(s)); +} + // ================================================================================= // Print a yaml file to 6 levels // ================================================================================= @@ -313,7 +372,10 @@ void parse_yaml(Yaml::Node& root, SimulationParameters_t& SimulationParamaters, std::cout << "Parsing YAML regions:" << std::endl; } // parse the region yaml text into a vector of region_fills - parse_regions(root, SimulationParamaters.region_fills, SimulationParamaters.region_fills_host); + parse_regions(root, + SimulationParamaters.region_fills, + SimulationParamaters.region_fills_host); + if (VERBOSE) { printf("\n"); @@ -559,11 +621,6 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) if (VERBOSE) { std::cout << "\tsource = " << source << std::endl; } - - if (mesh_input.source == mesh_input::generate && !mesh_input.file_path.empty()) { - std::cout << "ERROR: When the mesh source is set to generate, a mesh file cannot be passed in" << std::endl; - exit(0); - } } else{ std::cout << "ERROR: invalid mesh option input in YAML file: " << source << std::endl; @@ -612,18 +669,10 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) std::cout << "\tfile_path = " << path << std::endl; } - mesh_input.file_path = path; - - if (mesh_input.source == mesh_input::file && mesh_input.file_path.empty()) { - std::cout << "ERROR: When the mesh source is a file, a file_path must be set to point to the mesh file" << std::endl; - std::cout << "A mesh can either be generated or read in from a file, but not both" << std::endl; - } + // absolute path to file or local to the director where exe is run + mesh_input.file_path = path; - if (mesh_input.source == mesh_input::generate) { - std::cout << "ERROR: When the mesh source is set to generate, a mesh file cannot be passed in" << std::endl; - exit(0); - } - } // file path + } // end file path // Origin for the mesh else if (a_word.compare("origin") == 0) { std::string origin = root["mesh_options"][a_word].As(); @@ -799,13 +848,38 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) } throw std::runtime_error("**** Mesh Not Understood ****"); } + + + + // ----------------------------------------------- + // check for consistency in input settings + + if (mesh_input.source == mesh_input::file && mesh_input.file_path.empty()) { + std::cout << "ERROR: When the mesh source is a file, a file_path must be set to point to the mesh file" << std::endl; + std::cout << "A mesh can either be generated or read in from a file, but not both" << std::endl; + } + + if (mesh_input.source == mesh_input::generate) { + + // check to see if a file path was set + if(mesh_input.file_path.size()>0){ + std::cout << "ERROR: When the mesh source is set to generate, a mesh file cannot be passed in" << std::endl; + exit(0); + } + + } + + // ----------------------------------------------- + + } // end user_mesh_inputs } // end of parse mesh options // ================================================================================= // Parse Output options // ================================================================================= -void parse_output_options(Yaml::Node& root, output_options_t& output_options) +void parse_output_options(Yaml::Node& root, + output_options_t& output_options) { Yaml::Node& out_opts = root["output_options"]; @@ -900,7 +974,8 @@ void parse_output_options(Yaml::Node& root, output_options_t& output_options) // ================================================================================= void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills, - CArray ®ion_fills_host) + CArray& region_fills_host) + { Yaml::Node& region_yaml = root["regions"]; @@ -909,7 +984,6 @@ void parse_regions(Yaml::Node& root, region_fills = CArrayKokkos(num_regions , "sim_param.region_fills"); region_fills_host = CArray(num_regions); - // loop over the fill regions specified for (int reg_id = 0; reg_id < num_regions; reg_id++) { // read the variables names @@ -1325,7 +1399,9 @@ void parse_regions(Yaml::Node& root, std::cout << "\tfile_path = " << path << std::endl; } + // absolute path to file or local to the director where exe is run region_fills_host(reg_id).file_path = path; // saving the absolute file path + } // end file path // @@ -1402,7 +1478,12 @@ void parse_regions(Yaml::Node& root, "********************************************************************************************\n"); } }); - } + + } // end if + + // ----------------------------------------------- + + } // end loop over regions } // end of function to parse region @@ -1533,6 +1614,8 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) RUN({ Materials.MaterialEnums(mat_id).EOSType = model::decoupledEOSType; }); + + Materials.MaterialEnums.host(mat_id).EOSType = model::decoupledEOSType; break; case model::coupledEOSType: @@ -1540,12 +1623,14 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) RUN({ Materials.MaterialEnums(mat_id).EOSType = model::coupledEOSType; }); + Materials.MaterialEnums.host(mat_id).EOSType = model::coupledEOSType; break; default: RUN({ Materials.MaterialEnums(mat_id).EOSType = model::noEOSType; }); + Materials.MaterialEnums.host(mat_id).EOSType = model::noEOSType; std::cout << "ERROR: No valid EOS type input " << std::endl; std::cout << "Valid EOS types are: " << std::endl; @@ -1639,6 +1724,7 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) RUN({ Materials.MaterialEnums(mat_id).StrengthType = model::noStrengthType; }); + Materials.MaterialEnums.host(mat_id).StrengthType = model::noStrengthType; if (VERBOSE) { std::cout << "\tstrength_model_type_type = " << strength_model_type << std::endl; } @@ -1648,6 +1734,7 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) RUN({ Materials.MaterialEnums(mat_id).StrengthType = model::incrementBased; }); + Materials.MaterialEnums.host(mat_id).StrengthType = model::incrementBased; if (VERBOSE) { std::cout << "\tstrength_model_type = " << strength_model_type << std::endl; @@ -1657,6 +1744,8 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) RUN({ Materials.MaterialEnums(mat_id).StrengthType = model::stateBased; }); + Materials.MaterialEnums.host(mat_id).StrengthType = model::stateBased; + std::cout << "ERROR: state_based models not yet defined: " << std::endl; throw std::runtime_error("**** ERROR: state_based models not yet defined ****"); if (VERBOSE) { @@ -1723,12 +1812,16 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) // erosion_model_map[erosion_model] returns enum value, e.g., model::erosion switch(erosion_model_map[erosion_model]){ case model::basicErosion: + Materials.MaterialEnums.host(mat_id).ErosionModels = model::basicErosion; RUN({ + Materials.MaterialEnums(mat_id).ErosionModels = model::basicErosion; Materials.MaterialFunctions(mat_id).erode = &BasicErosionModel::erode; }); break; case model::noErosion: + Materials.MaterialEnums.host(mat_id).ErosionModels = model::noErosion; RUN({ + Materials.MaterialEnums(mat_id).ErosionModels = model::noErosion; Materials.MaterialFunctions(mat_id).erode = &NoErosionModel::erode; }); break; @@ -1750,15 +1843,6 @@ void parse_materials(Yaml::Node& root, Material_t& Materials) } // end if } // erosion model variables - else if (a_word.compare("void_mat_id") == 0) { - double void_mat_id = root["materials"][mat_id]["material"]["void_mat_id"].As(); - if (VERBOSE) { - std::cout << "\tvoid_mat_id = " << void_mat_id << std::endl; - } - RUN({ - Materials.MaterialFunctions(mat_id).void_mat_id = void_mat_id; - }); - } // blank_mat_id else if (a_word.compare("erode_tension_val") == 0) { double erode_tension_val = root["materials"][mat_id]["material"]["erode_tension_val"].As(); if (VERBOSE) { diff --git a/single-node-refactor/src/input/parse_yaml.h b/single-node-refactor/src/input/parse_yaml.h index af11ef464..2121e7f72 100644 --- a/single-node-refactor/src/input/parse_yaml.h +++ b/single-node-refactor/src/input/parse_yaml.h @@ -74,6 +74,25 @@ static bool DoesPathExist(const std::string& s) return (stat(s.c_str(), &buffer) == 0); } +// ============================================================================== +// Functions to parse strings +// ============================================================================== +// for string delimiter parsing +std::vector split(std::string s, std::string delimiter); + +// retrieves multiple values between [ ] +std::vector extract_list(std::string str); + +const std::string WHITESPACE = " "; + +std::string ltrim(const std::string& s); + +std::string rtrim(const std::string& s); + +std::string trim(const std::string& s); + + + // for string delimiter parsing std::vector exact_array_values(std::string s, std::string delimiter); @@ -106,7 +125,9 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input); void parse_output_options(Yaml::Node& root, output_options_t& output_options); // parse the region text -void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills, CArray ®ion_fills_host); +void parse_regions(Yaml::Node& root, + CArrayKokkos& region_fills, + CArray& region_fills_host); // parse the region text void parse_materials(Yaml::Node& root, Material_t& Materials); diff --git a/single-node-refactor/src/material_models/eos/gamma_law_eos.h b/single-node-refactor/src/material_models/eos/gamma_law_eos.h index fb6e37080..f7f5ae415 100644 --- a/single-node-refactor/src/material_models/eos/gamma_law_eos.h +++ b/single-node-refactor/src/material_models/eos/gamma_law_eos.h @@ -71,11 +71,11 @@ namespace GammaLawGasEOSModel { }; // host side function - static void initialize(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, + static void initialize(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, const size_t mat_id, - const DCArrayKokkos& elem_sspd, + const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const RaggedRightArrayKokkos &eos_global_vars, @@ -88,12 +88,12 @@ namespace GammaLawGasEOSModel { } // end func KOKKOS_FUNCTION - static void calc_pressure(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, + static void calc_pressure(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const RaggedRightArrayKokkos &eos_global_vars) @@ -102,19 +102,19 @@ namespace GammaLawGasEOSModel { double gamma = eos_global_vars(mat_id, VarNames::gamma); // pressure - elem_pres(elem_gid) = (gamma - 1.0) * sie * den; + MaterialPoints_pres(mat_pt_lid) = (gamma - 1.0) * sie * den; return; } // end func KOKKOS_FUNCTION - static void calc_sound_speed(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, + static void calc_sound_speed(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const RaggedRightArrayKokkos &eos_global_vars) @@ -125,7 +125,7 @@ namespace GammaLawGasEOSModel { // sound speed - elem_sspd(elem_gid) = fmax(sqrt(gamma * (gamma - 1.0) * sie), csmin); + MaterialPoints_sspd(mat_pt_lid) = fmax(sqrt(gamma * (gamma - 1.0) * sie), csmin); return; diff --git a/single-node-refactor/src/material_models/eos/no_eos.h b/single-node-refactor/src/material_models/eos/no_eos.h index 9f71f3b19..6dd9f889f 100644 --- a/single-node-refactor/src/material_models/eos/no_eos.h +++ b/single-node-refactor/src/material_models/eos/no_eos.h @@ -56,29 +56,29 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace NoEOSModel { KOKKOS_FUNCTION - static void calc_pressure(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie, - const RaggedRightArrayKokkos &eos_global_vars) + static void calc_pressure(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, + const size_t mat_id, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, + const double den, + const double sie, + const RaggedRightArrayKokkos &eos_global_vars) { return; } // end func KOKKOS_FUNCTION - static void calc_sound_speed(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie, - const RaggedRightArrayKokkos &eos_global_vars) + static void calc_sound_speed(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, + const size_t mat_id, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, + const double den, + const double sie, + const RaggedRightArrayKokkos &eos_global_vars) { return; diff --git a/single-node-refactor/src/material_models/eos/user_defined_eos.h b/single-node-refactor/src/material_models/eos/user_defined_eos.h index 10a7f102b..d7132022f 100644 --- a/single-node-refactor/src/material_models/eos/user_defined_eos.h +++ b/single-node-refactor/src/material_models/eos/user_defined_eos.h @@ -60,12 +60,12 @@ namespace UserDefinedEOSModel { KOKKOS_FUNCTION - static void calc_pressure(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, + static void calc_pressure(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const RaggedRightArrayKokkos &eos_global_vars) @@ -82,12 +82,12 @@ namespace UserDefinedEOSModel } // end for user_eos_model KOKKOS_FUNCTION - static void calc_sound_speed(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, + static void calc_sound_speed(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, const double den, const double sie, const RaggedRightArrayKokkos &eos_global_vars) @@ -130,36 +130,36 @@ namespace UserDefinedEOSModel namespace NotionalEOSModel { KOKKOS_FUNCTION - static void calc_pressure(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie, - const RaggedRightArrayKokkos &eos_global_vars) + static void calc_pressure(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, + const size_t mat_id, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, + const double den, + const double sie, + const RaggedRightArrayKokkos &eos_global_vars) { // pressure of a void is 0 - elem_pres(elem_gid) = 0.0; + MaterialPoints_pres(mat_pt_lid) = 0.0; return; } // end func KOKKOS_FUNCTION - static void calc_sound_speed(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie, - const RaggedRightArrayKokkos &eos_global_vars) + static void calc_sound_speed(const DCArrayKokkos& MaterialPoints_pres, + const DCArrayKokkos& MaterialPoints_stress, + const size_t mat_pt_lid, + const size_t mat_id, + const DCArrayKokkos& MaterialPoints_state_vars, + const DCArrayKokkos& MaterialPoints_sspd, + const double den, + const double sie, + const RaggedRightArrayKokkos &eos_global_vars) { // sound speed of a void is 0, machine small must be used for CFL calculation - elem_sspd(elem_gid) = 1.0e-32; + MaterialPoints_sspd(mat_pt_lid) = 1.0e-32; return; } // end func diff --git a/single-node-refactor/src/material_models/eos/void_eos.h b/single-node-refactor/src/material_models/eos/void_eos.h index 745f59e61..d736d8aba 100644 --- a/single-node-refactor/src/material_models/eos/void_eos.h +++ b/single-node-refactor/src/material_models/eos/void_eos.h @@ -58,7 +58,7 @@ namespace VoidEOSModel { KOKKOS_FUNCTION static void calc_pressure(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, @@ -67,7 +67,7 @@ namespace VoidEOSModel { const RaggedRightArrayKokkos &eos_global_vars) { // pressure of a void is 0 - elem_pres(elem_gid) = 0.0; + elem_pres(mat_pt_lid) = 0.0; return; } // end func @@ -75,7 +75,7 @@ namespace VoidEOSModel { KOKKOS_FUNCTION static void calc_sound_speed(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, @@ -85,7 +85,7 @@ namespace VoidEOSModel { { // sound speed of a void is 0, machine small must be used for CFL calculation - elem_sspd(elem_gid) = 1.0e-32; + elem_sspd(mat_pt_lid) = 1.0e-32; return; } // end func diff --git a/single-node-refactor/src/material_models/erosion/basic_erosion.h b/single-node-refactor/src/material_models/erosion/basic_erosion.h index e609f16a6..6f1f0dd0b 100644 --- a/single-node-refactor/src/material_models/erosion/basic_erosion.h +++ b/single-node-refactor/src/material_models/erosion/basic_erosion.h @@ -37,31 +37,26 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // ----------------------------------------------------------------------------- -// This is the basic element Erosion model +// This is the basic material pt Erosion model // ------------------------------------------------------------------------------ namespace BasicErosionModel { KOKKOS_FUNCTION - static void erode (const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const DCArrayKokkos& elem_eroded, - const DCArrayKokkos& elem_mat_id, - const size_t elem_gid, - const size_t void_mat_id, + static void erode (const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& MaterialPoints_stress, + const double MaterialPoint_pres, + const double MaterialPoint_den, + const double MaterialPoint_sie, + const double MaterialPoint_sspd, const double erode_tension_val, const double erode_density_val, - const DCArrayKokkos& elem_sspd, - const DCArrayKokkos& elem_den, - const double sie) + const size_t mat_point_lid) { - // simple model based on tension and density - if (elem_pres(elem_gid) <= erode_tension_val - || elem_den(elem_gid) <= erode_density_val) { + if (MaterialPoint_pres <= erode_tension_val || MaterialPoint_den <= erode_density_val) { - elem_mat_id(elem_gid) = void_mat_id; - elem_eroded(elem_gid) = true; + MaterialPoints_eroded(mat_point_lid) = true; } // end if diff --git a/single-node-refactor/src/material_models/erosion/no_erosion.h b/single-node-refactor/src/material_models/erosion/no_erosion.h index 45c6d976a..4190c84c5 100644 --- a/single-node-refactor/src/material_models/erosion/no_erosion.h +++ b/single-node-refactor/src/material_models/erosion/no_erosion.h @@ -42,17 +42,15 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace NoErosionModel { KOKKOS_FUNCTION - static void erode (const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const DCArrayKokkos& elem_eroded, - const DCArrayKokkos& elem_mat_id, - const size_t elem_gid, - const size_t void_mat_id, + static void erode (const DCArrayKokkos& MaterialPoints_eroded, + const DCArrayKokkos& MaterialPoints_stress, + const double MaterialPoint_pres, + const double MaterialPoint_den, + const double MaterialPoint_sie, + const double MaterialPoint_sspd, const double erode_tension_val, const double erode_density_val, - const DCArrayKokkos& elem_sspd, - const DCArrayKokkos& elem_den, - const double sie) + const size_t mat_point_lid) { return; diff --git a/single-node-refactor/src/material_models/strength/no_strength.h b/single-node-refactor/src/material_models/strength/no_strength.h index 957c9c37f..35ca884e4 100644 --- a/single-node-refactor/src/material_models/strength/no_strength.h +++ b/single-node-refactor/src/material_models/strength/no_strength.h @@ -67,7 +67,7 @@ namespace UserDefinedStrengthModel { KOKKOS_FUNCTION static void calc_stress(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, @@ -102,7 +102,7 @@ namespace NoStrengthModel { KOKKOS_FUNCTION static void calc_stress(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, @@ -146,7 +146,7 @@ namespace NotionalStrengthModel { KOKKOS_FUNCTION static void calc_stress(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, diff --git a/single-node-refactor/src/material_models/strength/user_defined_strength.h b/single-node-refactor/src/material_models/strength/user_defined_strength.h index 8922cb007..c194c8a1a 100644 --- a/single-node-refactor/src/material_models/strength/user_defined_strength.h +++ b/single-node-refactor/src/material_models/strength/user_defined_strength.h @@ -67,7 +67,7 @@ namespace UserDefinedStrengthModel { KOKKOS_FUNCTION static void calc_stress(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, @@ -120,7 +120,7 @@ namespace NotionalStrengthModel { KOKKOS_FUNCTION static void calc_stress(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, - const size_t elem_gid, + const size_t mat_pt_lid, const size_t mat_id, const DCArrayKokkos& elem_state_vars, const DCArrayKokkos& elem_sspd, diff --git a/single-node-refactor/src/solver.h b/single-node-refactor/src/solver.h index 1c107e64e..6ebbcf0a0 100644 --- a/single-node-refactor/src/solver.h +++ b/single-node-refactor/src/solver.h @@ -43,6 +43,7 @@ #include "material.h" #include "region.h" #include "boundary_conditions.h" +#include "dynamic_options.h" struct SimulationParameters_t; @@ -56,25 +57,21 @@ class Solver virtual void initialize(SimulationParameters_t& SimulationParamaters, Material_t& Materials, - BoundaryCondition_t& Boundary) const = 0; + mesh_t& mesh, + BoundaryCondition_t& Boundary, + State_t& State) const = 0; virtual void setup(SimulationParameters_t& SimulationParamaters, Material_t& Materials, - BoundaryCondition_t& Boundary, mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) const = 0; + BoundaryCondition_t& Boundary, + State_t& State) = 0; virtual void execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, - BoundaryCondition_t& Boundary, + BoundaryCondition_t& BoundaryConditions, mesh_t& mesh, - node_t& node, - MaterialPoint_t& MaterialPoints, - GaussPoint_t& GaussPoints, - corner_t& corner) = 0; + State_t& State) = 0; virtual void finalize(SimulationParameters_t& SimulationParamaters, Material_t& Materials,