diff --git a/apps/multiphysics/src/CMakeLists.txt b/apps/multiphysics/src/CMakeLists.txt index 21c47a527..95cbc76d9 100755 --- a/apps/multiphysics/src/CMakeLists.txt +++ b/apps/multiphysics/src/CMakeLists.txt @@ -61,7 +61,8 @@ include_directories(material_models/artificial_viscosity) include_directories(material_models/eos) include_directories(material_models/erosion) include_directories(material_models/fracture) -include_directories(material_models/strength) +include_directories(material_models/transient_strength) +include_directories(material_models/quasi_static_strength) include_directories(material_models/equilibration) include_directories(material_models/equilibration/include) include_directories(material_models/tabular/include) @@ -75,8 +76,10 @@ add_subdirectory(material_models/tabular) # Boundary conditions message("\n ****** ADDING BOUNDARY CONDITIONS ******** \n ") include_directories(boundary_conditions/velocity) +include_directories(boundary_conditions/displacement) include_directories(boundary_conditions/temperature) include_directories(boundary_conditions/stress) +include_directories(boundary_conditions/quasi_static_stress) diff --git a/apps/multiphysics/src/Solvers/SGH_solver_3D/include/sgh_solver_3D.hpp b/apps/multiphysics/src/Solvers/SGH_solver_3D/include/sgh_solver_3D.hpp index 9e68bdb08..121739847 100644 --- a/apps/multiphysics/src/Solvers/SGH_solver_3D/include/sgh_solver_3D.hpp +++ b/apps/multiphysics/src/Solvers/SGH_solver_3D/include/sgh_solver_3D.hpp @@ -199,7 +199,8 @@ class SGH3D : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, swage::Mesh& mesh, - State_t& State) override; + State_t& State, + elements::fe_ref_elem_t& ref_elem) override; ///////////////////////////////////////////////////////////////////////////// /// diff --git a/apps/multiphysics/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp b/apps/multiphysics/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp index b1ef9e4df..7c583010b 100644 --- a/apps/multiphysics/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp +++ b/apps/multiphysics/src/Solvers/SGH_solver_3D/src/sgh_execute.cpp @@ -59,7 +59,8 @@ void SGH3D::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, swage::Mesh& mesh, - State_t& State) + State_t& State, + elements::fe_ref_elem_t& ref_elem) { // Get MPI ranks and num ranks diff --git a/apps/multiphysics/src/Solvers/SGH_solver_rz/include/sgh_solver_rz.hpp b/apps/multiphysics/src/Solvers/SGH_solver_rz/include/sgh_solver_rz.hpp index 3271bd074..43ebc6d98 100644 --- a/apps/multiphysics/src/Solvers/SGH_solver_rz/include/sgh_solver_rz.hpp +++ b/apps/multiphysics/src/Solvers/SGH_solver_rz/include/sgh_solver_rz.hpp @@ -164,7 +164,8 @@ class SGHRZ : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, swage::Mesh& mesh, - State_t& State) override; + State_t& State, + elements::fe_ref_elem_t& ref_elem) override; ///////////////////////////////////////////////////////////////////////////// /// diff --git a/apps/multiphysics/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp b/apps/multiphysics/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp index 58f5250e6..23d1c7fa8 100644 --- a/apps/multiphysics/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp +++ b/apps/multiphysics/src/Solvers/SGH_solver_rz/src/sgh_execute_rz.cpp @@ -54,7 +54,8 @@ void SGHRZ::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, swage::Mesh& mesh, - State_t& State) + State_t& State, + elements::fe_ref_elem_t& ref_elem) { double fuzz = SimulationParamaters.DynamicOptions.fuzz; diff --git a/apps/multiphysics/src/Solvers/SGTM_solver_3D/include/sgtm_solver_3D.hpp b/apps/multiphysics/src/Solvers/SGTM_solver_3D/include/sgtm_solver_3D.hpp index 2941b56a9..cdf1e69eb 100644 --- a/apps/multiphysics/src/Solvers/SGTM_solver_3D/include/sgtm_solver_3D.hpp +++ b/apps/multiphysics/src/Solvers/SGTM_solver_3D/include/sgtm_solver_3D.hpp @@ -202,7 +202,8 @@ class SGTM3D : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, swage::Mesh& mesh, - State_t& State) override; + State_t& State, + elements::fe_ref_elem_t& ref_elem) override; ///////////////////////////////////////////////////////////////////////////// /// diff --git a/apps/multiphysics/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp b/apps/multiphysics/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp index 3eed376d2..e421032db 100644 --- a/apps/multiphysics/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp +++ b/apps/multiphysics/src/Solvers/SGTM_solver_3D/src/sgtm_execute.cpp @@ -56,7 +56,8 @@ void SGTM3D::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, swage::Mesh& mesh, - State_t& State) + State_t& State, + elements::fe_ref_elem_t& ref_elem) { double fuzz = SimulationParamaters.DynamicOptions.fuzz; diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/CMakeLists.txt b/apps/multiphysics/src/Solvers/TLQS_solver_3D/CMakeLists.txt index 9719b742e..c4ce80983 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/CMakeLists.txt +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/CMakeLists.txt @@ -20,9 +20,14 @@ message("\n ****** ADDING TLQS SOURCE FILES ******** \n ") set(TLQS_3D_SRC_Files ${CMAKE_CURRENT_SOURCE_DIR}/src/boundary.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/elem_arrays.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/cgm_functions.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/tlqs_execute.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/time_integration.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/tlqs_setup.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/tlqs_initialize.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/post_process.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/qr_solver.cpp +${CMAKE_CURRENT_SOURCE_DIR}/src/chebyshev_smoothing.cpp PARENT_SCOPE ) \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp index 1711368a9..7b4cf68fd 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/include/tlqs_solver_3D.hpp @@ -58,6 +58,9 @@ namespace TLQS3D_State static const std::vector required_node_state = { node_state::coords, + node_state::coords_t0, + node_state::displacement, + node_state::mass }; // Gauss point state to be initialized for the TLQS solver @@ -70,14 +73,11 @@ namespace TLQS3D_State static const std::vector required_material_pt_state = { material_pt_state::density, - material_pt_state::pressure, material_pt_state::stress, - material_pt_state::specific_internal_energy, - material_pt_state::sound_speed, + material_pt_state::strain, material_pt_state::mass, material_pt_state::volume_fraction, material_pt_state::eroded_flag, - material_pt_state::shear_modulii }; // Material corner state to be initialized for the TLQS solver @@ -98,27 +98,14 @@ namespace TLQS3D_State // Node state that must be filled (setup) for the SGH solver static const std::vector required_fill_node_state = { - fill_node_state::velocity + //fill_node_state::displacement }; // Material point state that must be filled (setup) for the SGH solver // option A - static const std::vector required_optA_fill_material_pt_state = + static const std::vector required_fill_material_pt_state = { - fill_gauss_state::density, - fill_gauss_state::specific_internal_energy - }; - // option B - static const std::vector required_optB_fill_material_pt_state = - { - fill_gauss_state::density, - fill_gauss_state::internal_energy - }; - // option C - static const std::vector required_optC_fill_material_pt_state = - { - fill_gauss_state::density, - fill_gauss_state::stress + fill_gauss_state::density }; // ------------------------------------- @@ -127,13 +114,13 @@ namespace TLQS3D_State ///////////////////////////////////////////////////////////////////////////// /// -/// \class SGTM +/// \class TLQS /// -/// \brief Class for containing functions required to perform staggered grid -/// thermomechanical 3D Cartesian meshes. +/// \brief Class for containing functions required to perform Total Lagrangian +/// quasi-static 3D Cartesian meshes. /// /// This class contains the requisite functions requited to perform -/// staggered grid thermomechanical heat transfer. +/// Total Lagrangian quasi-static mechanics. /// ///////////////////////////////////////////////////////////////////////////// class TLQS3D : public Solver @@ -202,7 +189,8 @@ class TLQS3D : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, swage::Mesh& mesh, - State_t& State) override; + State_t& State, + elements::fe_ref_elem_t& ref_elem) override; ///////////////////////////////////////////////////////////////////////////// /// @@ -224,10 +212,15 @@ class TLQS3D : public Solver // **** Functions defined in boundary.cpp **** // - void boundary_position(const swage::Mesh& mesh, - const BoundaryCondition_t& BoundaryConditions, - MPICArrayKokkos& node_vel, - const double time_value) const; + void boundary_displacement(const swage::Mesh& mesh, + const BoundaryCondition_t& BoundaryConditions, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end) const; // **** Functions defined in time_integration.cpp **** // void timestep_init( @@ -244,7 +237,170 @@ class TLQS3D : public Solver const size_t num_nodes, const size_t mat_id) const; - + // **** Functions defined in elem_arrays.cpp **** // + + // inputs: material_matrix, mesh map for nodes in the element, element id, node reference coordinates, node displacements, basis gradients wrt master element + // outputs: displacement gradient, inverse of jacobian, jacobian determinant, 2nd PK stress in current configuration + KOKKOS_FUNCTION + static void get_gradients( + const double material_matrix[6][6], + ViewCArrayKokkos & nodes_in_elem, + const MPICArrayKokkos & coords_t0, + const MPICArrayKokkos & displacement, + const CArrayKokkos & displacement_step, + ViewCArrayKokkos & gauss_point_grad_basis, + double grad_u[3][3], + double inv_J[3][3], + double& det_J, + double PK2_curr_config[6] + ); + + // inputs: material_matrix, displacement gradient, inverse Jacobian, basis gradients wrt master element, current PK2 stress + // outputs: updated element stiffness matrix, updated element force vector + KOKKOS_FUNCTION + static void tally_elem_arrays( + const double material_matrix[6][6], + const double grad_u[3][3], + const double inv_J[3][3], + ViewCArrayKokkos & gauss_point_grad_basis, + double gauss_point_weight, + const double PK2_curr_config[6], + ViewCArrayKokkos & Kel, + ViewCArrayKokkos & Fel, + const double vol_frac + ); + + // **** Functions defined in cgm_functions.cpp **** // + + // inputs: mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, F_elem, K_elem, displacement_iter + // outputs: initial cgm residual: r0 + void get_r0( + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& K_elem, + const CArrayKokkos& displacement_iter, + const CArrayKokkos& r0 + ); + + // inputs: mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, K_elem, rk, p + // outputs: alpha for cgm + double get_alpha( + const size_t num_nodes, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const CArrayKokkos& K_elem, + const double rktrk, + const CArrayKokkos& p + ); + + void get_rkp1( + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const CArrayKokkos& K_elem, + const CArrayKokkos& rk, + const CArrayKokkos& p, + const double alpha, + const CArrayKokkos& rkp1 + ); + + // **** Functions defined in post_process.cpp **** // + + // postprocessing function: updates stress and strain material point fields + KOKKOS_FUNCTION + static void post_process( + const double material_matrix[6][6], + ViewCArrayKokkos & nodes_in_elem, + const MPICArrayKokkos & coords_t0, + const MPICArrayKokkos & displacement, + ViewCArrayKokkos & gauss_point_grad_basis, + ViewCArrayKokkos & stress, + ViewCArrayKokkos & strain + ); + + // **** Functions defined in qr_solver.cpp **** // + void QR_backsub(const CArrayKokkos &R, + const CArrayKokkos &y, + DCArrayKokkos &x + ); + + void QR_decompose(const CArrayKokkos &A, + FArrayKokkos &Q, + CArrayKokkos &R + ); + + double QR_determinant(const FArrayKokkos &Q, + const CArrayKokkos &R + ); + + void QR_solver(const CArrayKokkos &A, + const CArrayKokkos &b, + DCArrayKokkos &x + ); + + // **** Functions defined in additive_schwarz_preconditioning.cpp **** // + /* void get_z0( + const size_t num_nodes, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const RaggedRightArrayKokkos& elems_in_node, + const CArrayKokkos& K_elem, + const CArrayKokkos& rk, + CArrayKokkos& zk + ); + + void get_zkp1( + const size_t num_nodes, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const RaggedRightArrayKokkos& elems_in_node, + const CArrayKokkos& K_elem, + const CArrayKokkos& rk, + CArrayKokkos& zk + ); */ + + // **** Functions defined in chebyshev_smoothing.cpp **** // + void apply_chebyshev_preconditioner(const CArrayKokkos& rk, + const CArrayKokkos& zkp1, + const CArrayKokkos& D_inv, + const CArrayKokkos& zk, + const CArrayKokkos& delta_z, + const CArrayKokkos& temporary, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const double alpha, + const double beta, + const int degree + ); + + void get_diagonal_inverse(CArrayKokkos& D_inv, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem + ); + + void get_chebyshev_bounds(double& alpha, + double& beta, + const CArrayKokkos& D_inv, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + CArrayKokkos& v_scratch, + CArrayKokkos& w_scratch, + const int max_iters + ); + }; diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/boundary.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/boundary.cpp index c8f85140e..78da6c4bd 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/boundary.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/boundary.cpp @@ -41,44 +41,53 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. /// /// \fn boundary_position /// -/// \brief Evolves the boundary according to a given velocity +/// \brief Evolves the boundary according to a given displacement /// /// \param mesh The simulation mesh /// \param BoundaryConditions Boundary contains arrays of information about BCs -/// \param node_vel The nodal velocity array +/// \param node_disp The nodal displacement array /// \param time_value The current simulation time /// ///////////////////////////////////////////////////////////////////////////// -void TLQS3D::boundary_position(const swage::Mesh& mesh, +void TLQS3D::boundary_displacement(const swage::Mesh& mesh, const BoundaryCondition_t& BoundaryConditions, - MPICArrayKokkos& node_vel, - const double time_value) const + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end) const { - // size_t num_pos_bdy_sets = BoundaryConditions.num_pos_bdy_sets_in_solver.host(this->solver_id); + size_t num_disp_bdy_sets = BoundaryConditions.num_disp_bdy_sets_in_solver.host(this->solver_id); + //std::cout << "NUM DISP BDY SETS" << num_disp_bdy_sets << std::endl; + // Loop over the displacement boundary sets + for (size_t bc_lid = 0; bc_lid < num_disp_bdy_sets; bc_lid++) { + + size_t bdy_set = BoundaryConditions.disp_bdy_sets_in_solver.host(this->solver_id, bc_lid); + //std::cout << "NUM BDY NODES IN SET" << mesh.num_bdy_nodes_in_set.host(bdy_set) << std::endl << std::endl; + // Loop over boundary nodes in a boundary set + FOR_ALL(bdy_node_lid, 0, mesh.num_bdy_nodes_in_set.host(bdy_set), { + // get the global index for this node on the boundary + size_t bdy_node_gid = mesh.bdy_nodes_in_set(bdy_set, bdy_node_lid); - // // Loop over the velocity boundary sets - // for (size_t bc_lid = 0; bc_lid < num_pos_bdy_sets; bc_lid++) { - - // size_t bdy_set = BoundaryConditions.vel_bdy_sets_in_solver.host(this->solver_id, bc_lid); - - // // Loop over boundary nodes in a boundary set - // FOR_ALL(bdy_node_lid, 0, mesh.num_bdy_nodes_in_set.host(bdy_set), { - // // get the global index for this node on the boundary - // size_t bdy_node_gid = mesh.bdy_nodes_in_set(bdy_set, bdy_node_lid); - - // // evaluate velocity on this boundary node - // BoundaryConditions.BoundaryConditionFunctions(bdy_set).position( - // mesh, - // BoundaryConditions.BoundaryConditionEnums, - // BoundaryConditions.velocity_bc_global_vars, - // BoundaryConditions.bc_state_vars, - // node_vel, - // time_value, - // 1, // rk_stage isn't used - // bdy_node_gid, - // bdy_set); - // }); // end for bdy_node_lid - // } // end for bdy_set + // evaluate displacement on this boundary node + BoundaryConditions.BoundaryConditionFunctions(bdy_set).displacement( + mesh, + BoundaryConditions.BoundaryConditionEnums, + BoundaryConditions.displacement_bc_global_vars, + BoundaryConditions.bc_state_vars, + K_elem, + F_elem, + displacement_step, + dt, + time_value, + time_start, + time_end, + bdy_node_gid, + bdy_set); + }); // end for bdy_node_lid + } // end for bdy_set return; -} // end boundary_position function \ No newline at end of file +} // end boundary_displacement function \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp new file mode 100644 index 000000000..9e6f19880 --- /dev/null +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/cgm_functions.cpp @@ -0,0 +1,182 @@ +/********************************************************************************************** +� 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 "tlqs_solver_3D.hpp" + +void TLQS3D::get_r0( + const size_t num_nodes, + const RaggedRightArrayKokkos & elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos & nodes_in_elem, + const CArrayKokkos & F_elem, + const CArrayKokkos & K_elem, + const CArrayKokkos & displacement_iter, + const CArrayKokkos & r0 +) +{ + // getting r0 = (02F - 01F) - K * displacement_iter + FOR_ALL(node_gid, 0, num_nodes, { + const size_t num_elems_in_node = elems_in_node.stride(node_gid); + + for (size_t p = 0; p < 3; p++) { + const size_t global_dof = 3 * node_gid + p; + double val = 0.0; + + // Sum contributions from all elements containing this node + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); + + // Find local index of this node within the element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p; + + // F_elem contribution + val += F_elem(elem_gid, local_dof); + //std::cout << "F_ELEM: " << F_elem(elem_gid, local_dof) << std::endl; + + // Subtract K_elem * displacement_iter + for (size_t b = 0; b < num_nodes_in_elem; b++) { + const size_t node_gid_b = nodes_in_elem(elem_gid, b); + for (size_t q = 0; q < 3; q++) { + const size_t local_dof_b = 3 * b + q; + const size_t global_dof_b = 3 * node_gid_b + q; + val -= K_elem(elem_gid, local_dof, local_dof_b) * displacement_iter(global_dof_b); + //std::cout << "K_ELEM: " << K_elem(elem_gid, local_dof, local_dof_b) << std::endl; + } + } + } + + r0(global_dof) = val; + } + }); + Kokkos::fence(); +} // end get_r0 + +double TLQS3D::get_alpha( + const size_t num_nodes, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const CArrayKokkos& K_elem, + const double rktrk, + const CArrayKokkos& p) +{ + + // denominator: p^T * K * p + // first compute Kp = K * p via assembly-free matvec + // then dot with p + double ptkp = 0.0; + double loc_ptkp = 0.0; + FOR_REDUCE_SUM(elem_gid, 0, K_elem.dims(0), loc_ptkp, { + + for (size_t a = 0; a < num_nodes_in_elem; a++) { + const size_t node_gid_a = nodes_in_elem(elem_gid, a); + for (size_t p_dir = 0; p_dir < 3; p_dir++) { + const size_t local_dof_a = 3 * a + p_dir; + const size_t global_dof_a = 3 * node_gid_a + p_dir; + + double Kp_val = 0.0; + for (size_t b = 0; b < num_nodes_in_elem; b++) { + const size_t node_gid_b = nodes_in_elem(elem_gid, b); + for (size_t q = 0; q < 3; q++) { + const size_t local_dof_b = 3 * b + q; + const size_t global_dof_b = 3 * node_gid_b + q; + Kp_val += K_elem(elem_gid, local_dof_a, local_dof_b) * p(global_dof_b); + } + } + loc_ptkp += p(global_dof_a) * Kp_val; + } + } + }, ptkp); + Kokkos::fence(); + //std::cout << "PTKP: " << ptkp << std::endl; + + return rktrk / (ptkp+1E-16); +} // end get_alpha + +void TLQS3D::get_rkp1( + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const CArrayKokkos& K_elem, + const CArrayKokkos& rk, + const CArrayKokkos& p, + const double alpha, + const CArrayKokkos& rkp1) +{ + // r_{k+1} = r_k - alpha * K * p + FOR_ALL(node_gid, 0, num_nodes, { + const size_t num_elems_in_node = elems_in_node.stride(node_gid); + + for (size_t p_dir = 0; p_dir < 3; p_dir++) { + const size_t global_dof = 3 * node_gid + p_dir; + double Kp_val = 0.0; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); + + // find local index of this node within the element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p_dir; + + for (size_t b = 0; b < num_nodes_in_elem; b++) { + const size_t node_gid_b = nodes_in_elem(elem_gid, b); + for (size_t q = 0; q < 3; q++) { + const size_t local_dof_b = 3 * b + q; + const size_t global_dof_b = 3 * node_gid_b + q; + Kp_val += K_elem(elem_gid, local_dof, local_dof_b) * p(global_dof_b); + } + } + } + + rkp1(global_dof) = rk(global_dof) - alpha * Kp_val; + } + }); + Kokkos::fence(); +} // end get_rkp1 \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp new file mode 100644 index 000000000..cb76438ea --- /dev/null +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/chebyshev_smoothing.cpp @@ -0,0 +1,292 @@ +/********************************************************************************************** +© 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 "tlqs_solver_3D.hpp" + +/** + * @brief Applies a matrix-free Chebyshev polynomial preconditioner using a thread-safe, + * node-based gathering approach (no atomic operations required). + */ +void TLQS3D::apply_chebyshev_preconditioner(const CArrayKokkos& rk, + const CArrayKokkos& zkp1, + const CArrayKokkos& D_inv, + const CArrayKokkos& zk, + const CArrayKokkos& delta_z, + const CArrayKokkos& temporary, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + const double alpha, + const double beta, + const int degree) +{ + const size_t total_dofs = 3 * num_nodes; + + // Compute Chebyshev parameters based on spectral bounds + const double d = (beta + alpha) / 2.0; + const double c = (beta - alpha) / 2.0; + + // --- Step 1: Initialize the 3-term recurrence (Iteration k = 0) --- + FOR_ALL(i, 0, total_dofs, { + zk(i) = 0.0; + delta_z(i) = (1.0 / d) * D_inv(i) * rk(i); + zk(i) += delta_z(i); + }); + MATAR_FENCE(); + + double rho_prev = c / (2.0 * d); + + // --- Step 2: Recurrence Loop (Iteration k = 1 to degree-1) --- + for (int k = 1; k < degree; ++k) { + double rho_k = 1.0 / (2.0 * d / c - rho_prev); + double gamma_k = rho_k * 2.0 / c; + + // --- THREAD-SAFE MATRIX-FREE MULTIPLICATION (temporary = K * zk) --- + // Parallelized by node; each thread computes its exact global DOF value. + FOR_ALL(node_gid, 0, num_nodes, { + const size_t num_elems_in_node = elems_in_node.stride(node_gid); + + for (size_t p = 0; p < 3; p++) { + const size_t global_dof = 3 * node_gid + p; + double val = 0.0; + + // Sum contributions from all elements containing this global node + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); + + // Find local index of this node within the current element + size_t local_node_lid = num_nodes_in_elem; + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p; + + // Contract K_elem row with the incoming Chebyshev vector (zk) + for (size_t b = 0; b < num_nodes_in_elem; b++) { + const size_t node_gid_b = nodes_in_elem(elem_gid, b); + + for (size_t q = 0; q < 3; q++) { + const size_t local_dof_b = 3 * b + q; + const size_t global_dof_b = 3 * node_gid_b + q; + + val += K_elem(elem_gid, local_dof, local_dof_b) * zk(global_dof_b); + } + } + } + + // Directly assign to scratch array without atomics or clearing passes + temporary(global_dof) = val; + } + }); + MATAR_FENCE(); + // ----------------------------------------------------------------- + + // Perform the vector updates using MATAR 1D indexing + FOR_ALL(i, 0, total_dofs, { + delta_z(i) = rho_k * delta_z(i) + gamma_k * D_inv(i) * (rk(i) - temporary(i)); + zk(i) += delta_z(i); + }); + MATAR_FENCE(); + + rho_prev = rho_k; + } + + // --- Step 3: Finalize Output --- + FOR_ALL(i, 0, total_dofs, { + zkp1(i) = zk(i); + }); + MATAR_FENCE(); +} + +/** + * @brief Computes the inverse of the global stiffness matrix diagonal (D_inv) + * using a thread-safe, node-based assembly approach. + * * @param D_inv Output vector tracking 1.0 / (K_ii + epsilon) + * @param K_elem Element stiffness arrays (num_elems, 3*npe, 3*npe) + * @param num_nodes Total number of global nodes in the mesh + * @param elems_in_node Ragged array mapping global node ID to local elements + * @param num_nodes_in_elem Number of nodes per element (e.g., 64 for cubic hex) + * @param nodes_in_elem Array mapping element ID to its global node IDs + */ +void TLQS3D::get_diagonal_inverse(CArrayKokkos& D_inv, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem) +{ + FOR_ALL(node_gid, 0, num_nodes, { + const size_t num_elems_in_node = elems_in_node.stride(node_gid); + + for (size_t p_dir = 0; p_dir < 3; p_dir++) { + const size_t global_dof = 3 * node_gid + p_dir; + double diag = 0.0; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); + + // Find local index of this node within the element + size_t local_node_lid = num_nodes_in_elem; + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p_dir; + + // Accumulate diagonal entries (local_dof, local_dof) + diag += K_elem(elem_gid, local_dof, local_dof); + } + + // Invert with safety epsilon protection + D_inv(global_dof) = 1.0 / (diag + 1e-16); + } + }); + MATAR_FENCE(); +} + +/** + * @brief Computes the spectral bounds (alpha and beta) for the Chebyshev preconditioner + * using a performance-optimized Power Iteration loop with fused kernels and proper MATAR reductions. + */ +void TLQS3D::get_chebyshev_bounds(double& alpha, + double& beta, + const CArrayKokkos& D_inv, + const CArrayKokkos& K_elem, + const size_t num_nodes, + const RaggedRightArrayKokkos& elems_in_node, + const size_t num_nodes_in_elem, + const DCArrayKokkos& nodes_in_elem, + CArrayKokkos& v_scratch, + CArrayKokkos& w_scratch, + const int max_iters) +{ + const size_t total_dofs = 3 * num_nodes; + double lambda_max = 0.0; + + // Initialize initial guess vector v to 1.0 + FOR_ALL(i, 0, total_dofs, { + v_scratch(i) = 1.0; + }); + MATAR_FENCE(); + + // --- Power Iteration Loop --- + for (int iter = 0; iter < max_iters; ++iter) { + + // 1. FUSED KERNEL: Compute Matrix-Free Product AND Scale by D_inv directly + FOR_ALL(node_gid, 0, num_nodes, { + const size_t num_elems_in_node = elems_in_node.stride(node_gid); + + for (size_t p = 0; p < 3; p++) { + const size_t global_dof = 3 * node_gid + p; + double val = 0.0; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = elems_in_node(node_gid, elem_lid); + + size_t local_node_lid = num_nodes_in_elem; + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (nodes_in_elem(elem_gid, a) == node_gid) { + local_node_lid = a; + break; + } + } + + const size_t local_dof = 3 * local_node_lid + p; + + for (size_t b = 0; b < num_nodes_in_elem; b++) { + const size_t node_gid_b = nodes_in_elem(elem_gid, b); + + for (size_t q = 0; q < 3; q++) { + const size_t local_dof_b = 3 * b + q; + const size_t global_dof_b = 3 * node_gid_b + q; + + val += K_elem(elem_gid, local_dof, local_dof_b) * v_scratch(global_dof_b); + } + } + } + + // Fused operation: Scale the accumulated stiffness action by D_inv + w_scratch(global_dof) = val * D_inv(global_dof); + } + }); + MATAR_FENCE(); + + // 2. Compute Dot Products for Rayleigh Quotient using 6-arg FOR_REDUCE_SUM + double v_dot_w = 0.0; + double v_dot_v = 0.0; + + double local_v_dot_w = 0.0; + FOR_REDUCE_SUM(i, 0, total_dofs, local_v_dot_w, { + local_v_dot_w += v_scratch(i) * w_scratch(i); + }, v_dot_w); + + double local_v_dot_v = 0.0; + FOR_REDUCE_SUM(i, 0, total_dofs, local_v_dot_v, { + local_v_dot_v += v_scratch(i) * v_scratch(i); + }, v_dot_v); + MATAR_FENCE(); + + // Calculate current estimate of the maximum eigenvalue + lambda_max = v_dot_w / (v_dot_v + 1e-16); + + // 3. Normalize w vector to update our guess v: v = w / ||w||_2 + double w_norm2 = 0.0; + double local_w_norm2 = 0.0; + FOR_REDUCE_SUM(i, 0, total_dofs, local_w_norm2, { + local_w_norm2 += w_scratch(i) * w_scratch(i); + }, w_norm2); + MATAR_FENCE(); + + double inv_norm = 1.0 / (sqrt(w_norm2) + 1e-16); + + FOR_ALL(i, 0, total_dofs, { + v_scratch(i) = w_scratch(i) * inv_norm; + }); + MATAR_FENCE(); + } + + // --- 4. Apply Heuristic Bounding Box --- + beta = 1.1 * lambda_max; + alpha = beta / 10.0; +} \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp new file mode 100644 index 000000000..8fe059143 --- /dev/null +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/elem_arrays.cpp @@ -0,0 +1,303 @@ +/********************************************************************************************** +© 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 "tlqs_solver_3D.hpp" + + +KOKKOS_FUNCTION +void TLQS3D::get_gradients( + const double material_matrix[6][6], + ViewCArrayKokkos & nodes_in_elem, + const MPICArrayKokkos & coords_t0, + const MPICArrayKokkos & displacement, + const CArrayKokkos & displacement_step, + ViewCArrayKokkos & gauss_point_grad_basis, + double grad_u[3][3], + double inv_J[3][3], + double& det_J, + double PK2_curr_config[6]) +{ + // allocate and initialize Jacobian + // dX_i / dxi_j + double J[3][3]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + J[i][j] = 0.0; + } + } + + // *************************************************** + // FIX THIS TO LOOP CORRECTLY FOR C DATA LAYOUT + // *************************************************** + // get Jacobian at the mat point + for (int k = 0; k < gauss_point_grad_basis.dims(0); k++) { + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + J[j][i] += coords_t0(nodes_in_elem(k), j) * gauss_point_grad_basis(k, i); + } + } + } + + // get det(J) + det_J = J[0][0]*(J[1][1]*J[2][2] - J[2][1]*J[1][2]) - J[0][1]*(J[1][0]*J[2][2] - J[2][0]*J[1][2]) + J[0][2]*(J[1][0]*J[2][1] - J[2][0]*J[1][1]); + //std::cout << "DET_J: " << det_J << std::endl; + // get inv(J) + double adjoint[3][3]; + adjoint[0][0] = J[1][1]*J[2][2] - J[2][1]*J[1][2]; + adjoint[0][1] = -(J[0][1]*J[2][2] - J[2][1]*J[0][2]); + adjoint[0][2] = J[0][1]*J[1][2] - J[1][1]*J[0][2]; + adjoint[1][0] = -(J[1][0]*J[2][2] - J[2][0]*J[1][2]); + adjoint[1][1] = J[0][0]*J[2][2] - J[2][0]*J[0][2]; + adjoint[1][2] = -(J[0][0]*J[1][2] - J[1][0]*J[0][2]); + adjoint[2][0] = J[1][0]*J[2][1] - J[2][0]*J[1][1]; + adjoint[2][1] = -(J[0][0]*J[2][1] - J[2][0]*J[0][1]); + adjoint[2][2] = J[0][0]*J[1][1] - J[1][0]*J[0][1]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + inv_J[i][j] = adjoint[i][j] / det_J; + } + } + + // get grad(displacement) + // du_i / dX_j + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + grad_u[i][j] = 0.0; + } + } + + for (int k = 0; k < gauss_point_grad_basis.dims(0); k++) { + const size_t node_gid = nodes_in_elem(k); + + const double dpsig_k0 = inv_J[0][0]*gauss_point_grad_basis(k,0) + inv_J[0][1]*gauss_point_grad_basis(k,1) + inv_J[0][2]*gauss_point_grad_basis(k,2); + const double dpsig_k1 = inv_J[1][0]*gauss_point_grad_basis(k,0) + inv_J[1][1]*gauss_point_grad_basis(k,1) + inv_J[1][2]*gauss_point_grad_basis(k,2); + const double dpsig_k2 = inv_J[2][0]*gauss_point_grad_basis(k,0) + inv_J[2][1]*gauss_point_grad_basis(k,1) + inv_J[2][2]*gauss_point_grad_basis(k,2); + + for (int j = 0; j < 3; j++) { + const double u_total = displacement(node_gid, j) + displacement_step(3*node_gid + j); + grad_u[j][0] += u_total * dpsig_k0; + grad_u[j][1] += u_total * dpsig_k1; + grad_u[j][2] += u_total * dpsig_k2; + } + } + + // get second PK stress of current configuration + // *************************************************** + // WARNING: CURRENTLY ASSUMES ISOTROPIC LINEAR ELASTIC + // WARNING: NEED TO PUT THIS INTO A SEPARATE FUNCTION + // WARNING: POINTER CALLED CALC QUASI_STATIC STRESS + // *************************************************** + double current_strain[6]; // [Exx Eyy Ezz Eyz Exz Exy] + current_strain[0] = grad_u[0][0] + 0.5 * (grad_u[0][0]*grad_u[0][0] + grad_u[1][0]*grad_u[1][0] + grad_u[2][0]*grad_u[2][0]); + current_strain[1] = grad_u[1][1] + 0.5 * (grad_u[0][1]*grad_u[0][1] + grad_u[1][1]*grad_u[1][1] + grad_u[2][1]*grad_u[2][1]); + current_strain[2] = grad_u[2][2] + 0.5 * (grad_u[0][2]*grad_u[0][2] + grad_u[1][2]*grad_u[1][2] + grad_u[2][2]*grad_u[2][2]); + current_strain[3] = 0.5 * (grad_u[1][2] + grad_u[2][1] + grad_u[0][1]*grad_u[0][2] + grad_u[1][1]*grad_u[1][2] + grad_u[2][1]*grad_u[2][2]); + current_strain[4] = 0.5 * (grad_u[0][2] + grad_u[2][0] + grad_u[0][0]*grad_u[0][2] + grad_u[1][0]*grad_u[1][2] + grad_u[2][0]*grad_u[2][2]); + current_strain[5] = 0.5 * (grad_u[0][1] + grad_u[1][0] + grad_u[0][0]*grad_u[0][1] + grad_u[1][0]*grad_u[1][1] + grad_u[2][0]*grad_u[2][1]); + + // Normal stresses + PK2_curr_config[0] = material_matrix[0][0] * current_strain[0] + material_matrix[0][1] * current_strain[1] + material_matrix[0][2] * current_strain[2]; // Sxx + PK2_curr_config[1] = material_matrix[1][0] * current_strain[0] + material_matrix[1][1] * current_strain[1] + material_matrix[1][2] * current_strain[2]; // Syy + PK2_curr_config[2] = material_matrix[2][0] * current_strain[0] + material_matrix[2][1] * current_strain[1] + material_matrix[2][2] * current_strain[2]; // Szz + + // Shear stresses + PK2_curr_config[3] = material_matrix[3][3] * current_strain[3]; // Syz + PK2_curr_config[4] = material_matrix[4][4] * current_strain[4]; // Sxz + PK2_curr_config[5] = material_matrix[5][5] * current_strain[5]; // Sxy + +} // end get_gradients + +KOKKOS_FUNCTION +void TLQS3D::tally_elem_arrays( + const double material_matrix[6][6], + const double grad_u[3][3], + const double inv_J[3][3], + ViewCArrayKokkos & gauss_point_grad_basis, + double gauss_point_weight, + const double PK2_curr_config[6], + ViewCArrayKokkos & Kel, + ViewCArrayKokkos & Fel, + const double vol_frac) +{ + const int num_nodes = gauss_point_grad_basis.dims(0); + + // Unpack grad_u for readability + const double ux = grad_u[0][0]; + const double uy = grad_u[0][1]; + const double uz = grad_u[0][2]; + const double vx = grad_u[1][0]; + const double vy = grad_u[1][1]; + const double vz = grad_u[1][2]; + const double wx = grad_u[2][0]; + const double wy = grad_u[2][1]; + const double wz = grad_u[2][2]; + + // Build symmetric 3x3 stress tensor for K2 (extracted from Voigt PK2) + double S[3][3]; + S[0][0] = PK2_curr_config[0]; + S[0][1] = PK2_curr_config[5]; + S[0][2] = PK2_curr_config[4]; + S[1][0] = PK2_curr_config[5]; + S[1][1] = PK2_curr_config[1]; + S[1][2] = PK2_curr_config[3]; + S[2][0] = PK2_curr_config[4]; + S[2][1] = PK2_curr_config[3]; + S[2][2] = PK2_curr_config[2]; + + // temp arrays for forming element matrix + double B1_a[6][3]; // B1 values for node a + double CT_matmul_B1_a[6][3]; // C^T * B1_a = C * B1_a (C symmetric), using transpose for better memory access pattern + double S_mul_glob_grad_a[3]; // only needs to be size 3 because of the sparsity and repeating nature of [S]9x9 and B2 + double B1_b[6][3]; // B1 values for node b + + // looping through each node to avoid dynamic allocations + // outer loop of node a + for (int a = 0; a < num_nodes; a++) { + + // dpsig for node a: inv_J * grad_basis(a) + const double glob_grad_a_0 = inv_J[0][0]*gauss_point_grad_basis(a,0) + inv_J[0][1]*gauss_point_grad_basis(a,1) + inv_J[0][2]*gauss_point_grad_basis(a,2); + const double glob_grad_a_1 = inv_J[1][0]*gauss_point_grad_basis(a,0) + inv_J[1][1]*gauss_point_grad_basis(a,1) + inv_J[1][2]*gauss_point_grad_basis(a,2); + const double glob_grad_a_2 = inv_J[2][0]*gauss_point_grad_basis(a,0) + inv_J[2][1]*gauss_point_grad_basis(a,1) + inv_J[2][2]*gauss_point_grad_basis(a,2); + + // columns of B1 for node a + B1_a[0][0] = glob_grad_a_0*(1+ux); + B1_a[0][1] = glob_grad_a_0*vx; + B1_a[0][2] = glob_grad_a_0*wx; + + B1_a[1][0] = glob_grad_a_1*uy; + B1_a[1][1] = glob_grad_a_1*(1+vy); + B1_a[1][2] = glob_grad_a_1*wy; + + B1_a[2][0] = glob_grad_a_2*uz; + B1_a[2][1] = glob_grad_a_2*vz; + B1_a[2][2] = glob_grad_a_2*(1+wz); + + B1_a[3][0] = glob_grad_a_1*uz + glob_grad_a_2*uy; + B1_a[3][1] = glob_grad_a_1*vz + glob_grad_a_2*(1+vy); + B1_a[3][2] = glob_grad_a_1*(1+wz) + glob_grad_a_2*wy; + + B1_a[4][0] = glob_grad_a_2*(1+ux) + glob_grad_a_0*uz; + B1_a[4][1] = glob_grad_a_0*vz + glob_grad_a_2*vx; + B1_a[4][2] = glob_grad_a_0*(1+wz) + glob_grad_a_2*wx; + + B1_a[5][0] = glob_grad_a_1*(1+ux) + glob_grad_a_0*uy; + B1_a[5][1] = glob_grad_a_0*(1+vy) + glob_grad_a_1*vx; + B1_a[5][2] = glob_grad_a_0*wy + glob_grad_a_1*wx; + + // Precompute C^T * B1_a + // Used in inner loop as: K1(3a+p, 3b+q) = sum_m CtB1_a[m][p] * B1_b[m][q] + for (int i = 0; i < 6; i++) { + for (int j = 0; j < 3; j++) { + CT_matmul_B1_a[i][j] = 0.0; + } + } + for (int m = 0; m < 6; m++) { + for (int k = 0; k < 6; k++) { + for (int p = 0; p < 3; p++) { + CT_matmul_B1_a[m][p] += material_matrix[m][k] * B1_a[k][p]; + } + } + } + + // K2 simplification: The geometric stiffness block between nodes a and b is a + // diagonal 3x3 matrix where the non-zero terms equal grad(a) · S · grad(b). + // Precomputing S · grad(a) to save operations in the inner loop. + // Since S is symmetric, (Sg)⋅h=g⋅S⋅h + S_mul_glob_grad_a[0] = S[0][0]*glob_grad_a_0 + S[0][1]*glob_grad_a_1 + S[0][2]*glob_grad_a_2; + S_mul_glob_grad_a[1] = S[1][0]*glob_grad_a_0 + S[1][1]*glob_grad_a_1 + S[1][2]*glob_grad_a_2; + S_mul_glob_grad_a[2] = S[2][0]*glob_grad_a_0 + S[2][1]*glob_grad_a_1 + S[2][2]*glob_grad_a_2; + + // Fel: Fel(3a+p) -= weight * sum_k B1_a[k][p] * PK2[k] + // minus because Fel = boundarcy forces - current force in the element aka Fel = F02 - F01 + for (int p = 0; p < 3; p++) { + double fel_val = 0.0; + for (int k = 0; k < 6; k++) { + fel_val += B1_a[k][p] * PK2_curr_config[k]; + } + Fel(3*a + p) -= gauss_point_weight * fel_val * vol_frac; + } + + // Inner loop over b for Kel + for (int b = 0; b < num_nodes; b++) { + + // dpsig for node b + const double glob_grad_b_0 = inv_J[0][0]*gauss_point_grad_basis(b,0) + inv_J[0][1]*gauss_point_grad_basis(b,1) + inv_J[0][2]*gauss_point_grad_basis(b,2); + const double glob_grad_b_1 = inv_J[1][0]*gauss_point_grad_basis(b,0) + inv_J[1][1]*gauss_point_grad_basis(b,1) + inv_J[1][2]*gauss_point_grad_basis(b,2); + const double glob_grad_b_2 = inv_J[2][0]*gauss_point_grad_basis(b,0) + inv_J[2][1]*gauss_point_grad_basis(b,1) + inv_J[2][2]*gauss_point_grad_basis(b,2); + + // B1_b[6][3]: columns of B1 for node b + double B1_b[6][3]; + B1_b[0][0] = glob_grad_b_0*(1+ux); + B1_b[0][1] = glob_grad_b_0*vx; + B1_b[0][2] = glob_grad_b_0*wx; + + B1_b[1][0] = glob_grad_b_1*uy; + B1_b[1][1] = glob_grad_b_1*(1+vy); + B1_b[1][2] = glob_grad_b_1*wy; + + B1_b[2][0] = glob_grad_b_2*uz; + B1_b[2][1] = glob_grad_b_2*vz; + B1_b[2][2] = glob_grad_b_2*(1+wz); + + B1_b[3][0] = glob_grad_b_1*uz + glob_grad_b_2*uy; + B1_b[3][1] = glob_grad_b_1*vz + glob_grad_b_2*(1+vy); + B1_b[3][2] = glob_grad_b_1*(1+wz) + glob_grad_b_2*wy; + + B1_b[4][0] = glob_grad_b_2*(1+ux) + glob_grad_b_0*uz; + B1_b[4][1] = glob_grad_b_0*vz + glob_grad_b_2*vx; + B1_b[4][2] = glob_grad_b_0*(1+wz) + glob_grad_b_2*wx; + + B1_b[5][0] = glob_grad_b_1*(1+ux) + glob_grad_b_0*uy; + B1_b[5][1] = glob_grad_b_0*(1+vy) + glob_grad_b_1*vx; + B1_b[5][2] = glob_grad_b_0*wy + glob_grad_b_1*wx; + + // K2 scalar: dpsig_a · S3 · dpsig_b (same value along p==q diagonal) + const double k2_scalar = S_mul_glob_grad_a[0]*glob_grad_b_0 + S_mul_glob_grad_a[1]*glob_grad_b_1 + S_mul_glob_grad_a[2]*glob_grad_b_2; + + // Accumulate 3x3 block into Kel + for (int p = 0; p < 3; p++) { + for (int q = 0; q < 3; q++) { + double k1_val = 0.0; + for (int m = 0; m < 6; m++) { + k1_val += CT_matmul_B1_a[m][p] * B1_b[m][q]; + } + const double k2_val = (p == q) ? k2_scalar : 0.0; + Kel(3*a+p, 3*b+q) += gauss_point_weight * (k1_val + k2_val) * vol_frac; + } + } + } + + } + +} // end tally_elem_arrays \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/post_process.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/post_process.cpp new file mode 100644 index 000000000..55bf435af --- /dev/null +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/post_process.cpp @@ -0,0 +1,152 @@ +/********************************************************************************************** +© 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 "tlqs_solver_3D.hpp" + + +KOKKOS_FUNCTION +void TLQS3D::post_process( + const double material_matrix[6][6], + ViewCArrayKokkos & nodes_in_elem, + const MPICArrayKokkos & coords_t0, + const MPICArrayKokkos & displacement, + ViewCArrayKokkos & gauss_point_grad_basis, + ViewCArrayKokkos & stress, + ViewCArrayKokkos & strain) +{ + // allocate and initialize Jacobian + // dX_i / dxi_j + double J[3][3]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + J[i][j] = 0.0; + } + } + + // *************************************************** + // FIX THIS TO LOOP CORRECTLY FOR C DATA LAYOUT + // *************************************************** + // get Jacobian at the mat point + for (int k = 0; k < gauss_point_grad_basis.dims(0); k++) { + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + J[j][i] += coords_t0(nodes_in_elem(k), j) * gauss_point_grad_basis(k, i); + } + } + } + + // get det(J) + double det_J = J[0][0]*(J[1][1]*J[2][2] - J[2][1]*J[1][2]) - J[0][1]*(J[1][0]*J[2][2] - J[2][0]*J[1][2]) + J[0][2]*(J[1][0]*J[2][1] - J[2][0]*J[1][1]); + //std::cout << "DET_J: " << det_J << std::endl; + // get inv(J) + double adjoint[3][3]; + adjoint[0][0] = J[1][1]*J[2][2] - J[2][1]*J[1][2]; + adjoint[0][1] = -(J[0][1]*J[2][2] - J[2][1]*J[0][2]); + adjoint[0][2] = J[0][1]*J[1][2] - J[1][1]*J[0][2]; + adjoint[1][0] = -(J[1][0]*J[2][2] - J[2][0]*J[1][2]); + adjoint[1][1] = J[0][0]*J[2][2] - J[2][0]*J[0][2]; + adjoint[1][2] = -(J[0][0]*J[1][2] - J[1][0]*J[0][2]); + adjoint[2][0] = J[1][0]*J[2][1] - J[2][0]*J[1][1]; + adjoint[2][1] = -(J[0][0]*J[2][1] - J[2][0]*J[0][1]); + adjoint[2][2] = J[0][0]*J[1][1] - J[1][0]*J[0][1]; + + double inv_J[3][3]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + inv_J[i][j] = adjoint[i][j] / det_J; + } + } + + double grad_u[3][3]; + // get grad(displacement) + // du_i / dX_j + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + grad_u[i][j] = 0.0; + } + } + + for (int k = 0; k < gauss_point_grad_basis.dims(0); k++) { + const size_t node_gid = nodes_in_elem(k); + + const double dpsig_k0 = inv_J[0][0]*gauss_point_grad_basis(k,0) + inv_J[0][1]*gauss_point_grad_basis(k,1) + inv_J[0][2]*gauss_point_grad_basis(k,2); + const double dpsig_k1 = inv_J[1][0]*gauss_point_grad_basis(k,0) + inv_J[1][1]*gauss_point_grad_basis(k,1) + inv_J[1][2]*gauss_point_grad_basis(k,2); + const double dpsig_k2 = inv_J[2][0]*gauss_point_grad_basis(k,0) + inv_J[2][1]*gauss_point_grad_basis(k,1) + inv_J[2][2]*gauss_point_grad_basis(k,2); + + for (int j = 0; j < 3; j++) { + const double u_total = displacement(node_gid, j); + grad_u[j][0] += u_total * dpsig_k0; + grad_u[j][1] += u_total * dpsig_k1; + grad_u[j][2] += u_total * dpsig_k2; + } + } + + // get second PK stress and Green-Lagrange strain of current configuration + // *************************************************** + // WARNING: CURRENTLY ASSUMES ISOTROPIC LINEAR ELASTIC + // WARNING: NEED TO PUT THIS INTO A SEPARATE FUNCTION + // WARNING: POINTER CALLED CALC QUASI_STATIC STRESS + // *************************************************** + strain(0,0) = grad_u[0][0] + 0.5 * (grad_u[0][0]*grad_u[0][0] + grad_u[1][0]*grad_u[1][0] + grad_u[2][0]*grad_u[2][0]); + strain(1,1) = grad_u[1][1] + 0.5 * (grad_u[0][1]*grad_u[0][1] + grad_u[1][1]*grad_u[1][1] + grad_u[2][1]*grad_u[2][1]); + strain(2,2) = grad_u[2][2] + 0.5 * (grad_u[0][2]*grad_u[0][2] + grad_u[1][2]*grad_u[1][2] + grad_u[2][2]*grad_u[2][2]); + strain(1,2) = 0.5 * (grad_u[1][2] + grad_u[2][1] + grad_u[0][1]*grad_u[0][2] + grad_u[1][1]*grad_u[1][2] + grad_u[2][1]*grad_u[2][2]); + strain(2,1) = strain(1,2); + strain(0,2) = 0.5 * (grad_u[0][2] + grad_u[2][0] + grad_u[0][0]*grad_u[0][2] + grad_u[1][0]*grad_u[1][2] + grad_u[2][0]*grad_u[2][2]); + strain(2,0) = strain(0,2); + strain(0,1) = 0.5 * (grad_u[0][1] + grad_u[1][0] + grad_u[0][0]*grad_u[0][1] + grad_u[1][0]*grad_u[1][1] + grad_u[2][0]*grad_u[2][1]); + strain(1,0) = strain(0,1); + + /* for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + strain(i,j) = grad_u[i][j]; + } + } */ + + // Normal stresses + stress(0,0) = material_matrix[0][0] * strain(0,0) + material_matrix[0][1] * strain(1,1) + material_matrix[0][2] * strain(2,2); // Sxx + stress(1,1) = material_matrix[1][0] * strain(0,0) + material_matrix[1][1] * strain(1,1) + material_matrix[1][2] * strain(2,2); // Syy + stress(2,2) = material_matrix[2][0] * strain(0,0) + material_matrix[2][1] * strain(1,1) + material_matrix[2][2] * strain(2,2); // Szz + + // Shear stresses + stress(1,2) = material_matrix[3][3] * strain(1,2); // Syz + stress(2,1) = material_matrix[3][3] * strain(1,2); // Syz + + stress(0,2) = material_matrix[4][4] * strain(0,2); // Sxz + stress(2,0) = material_matrix[4][4] * strain(0,2); // Sxz + + stress(0,1) = material_matrix[5][5] * strain(0,1); // Sxy + stress(1,0) = material_matrix[5][5] * strain(0,1); // Sxy + +} // end post_process \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/qr_solver.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/qr_solver.cpp new file mode 100644 index 000000000..a9a1fcdbe --- /dev/null +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/qr_solver.cpp @@ -0,0 +1,291 @@ +#ifndef QR_H +#define QR_H +/********************************************************************************************** + © 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 +#include +#include + +#include "tlqs_solver_3D.hpp" + +#include "simulation_parameters.hpp" +#include "material.hpp" +#include "boundary_conditions.hpp" +#include "state.hpp" +#include "geometry_new.hpp" +#include "mesh_io.hpp" + +#include "matar.h" +using namespace mtr; + + +// Back substitution to solve Rx = y +void TLQS3D::QR_backsub(const CArrayKokkos &R, + const CArrayKokkos &y, + DCArrayKokkos &x) { + + size_t n = R.dims(0); + + for (int i = n - 1; i >= 0; --i) { + + RUN({ + x(i) = y(i); + }); + + double sum = 0.0; + double sum_lcl = 0.0; + + FOR_REDUCE_SUM(j, i + 1, n, + sum_lcl, { + sum_lcl -= R(i,j) * x(j); + }, sum); + + RUN({ + x(i) += sum; + x(i) /= R(i,i); + }); + + } // end for i + + return; + +} // end function + +// QR Decomposition using Modified Gram-Schmidt +void TLQS3D::QR_decompose(const CArrayKokkos &A, + FArrayKokkos &Q, + CArrayKokkos &R) { + + + const size_t m = A.dims(0); + const size_t n = A.dims(1); + + Q.set_values(0.0); + R.set_values(0.0); + + DCArrayKokkos v(n,m,"v"); + + // Copy columns of A to v, and taking transpose + FOR_ALL(i, 0, m, + j, 0, n, { + v(j, i) = A(i,j); + }); + Kokkos::fence(); + + + for (size_t i = 0; i < n; ++i) { + + // find the norm of a column in matrix v for row i + double tally = 0.0; + double tally_lcl = 0.0; + + + FOR_REDUCE_SUM(j, 0, m, + tally_lcl, { + tally_lcl += v(i,j) * v(i,j); + }, tally); + + Kokkos::fence(); + + RUN({ + R(i,i) = sqrt(tally); // row i norm + }); + // done with norm calc + + FOR_ALL(j, 0, m, { + Q(j,i) = v(i,j)/R(i,i); + }); + + Kokkos::fence(); + +// single parallelism +/* + FOR_ALL(jj, i+1, n, { + + R(i,jj) = 0.0; + + double sum=0; + + for(size_t k=0; k &Q, + const CArrayKokkos &R) +{ + //const size_t m = Q.dims(0); + const size_t n = Q.dims(1); + // Q(m,n,"Q"); + // R(n,n,"R"); + + double detR = 1.0; + int signQ = 1; + + // Accumulate det(R) and track sign adjustments + // for (size_t i = 0; i < n; ++i) { + // if (R(i,i) < 0.0) { + // signQ *= -1; // negative diagonal flips Q's sign + // detR *= -R(i,i); // store magnitude in detR + // } else { + // detR *= R(i,i); + // } + // } + + double prod_tally; + double prod_lcl = 1.0; + FOR_REDUCE_PRODUCT(i, 0, n, + prod_lcl, { + + if (R(i,i) < 0.0) { + prod_lcl *= -R(i,i); // store magnitude in detR + } else { + prod_lcl *= R(i,i); + } + }, prod_tally); // end j + Kokkos::fence(); + + + detR = prod_tally; + + + + prod_lcl = 1.0; + FOR_REDUCE_PRODUCT(i, 0, n, + prod_lcl, { + + if (R(i,i) < 0.0) { + prod_lcl *= -1; // negative diagonal flips Q's sign + } + }, prod_tally); // end j + + signQ = prod_tally; + + Kokkos::fence(); + + // Compute sign of det(Q) directly from Q if needed + // Here we infer sign from R adjustments only. + double detA = signQ * detR; + return detA; +} + + +// Solve for x in Ax = b using QR +// A[m,n] +// x[n] +// b[m] +void TLQS3D::QR_solver(const CArrayKokkos &A, + const CArrayKokkos &b, + DCArrayKokkos &x) { + + const size_t m = A.dims(0); + const size_t n = A.dims(1); + + FArrayKokkos Q(m,n,"Q"); + CArrayKokkos R(n,n,"R"); + CArrayKokkos y(n,"y"); + + QR_decompose(A, Q, R); + + + // Compute Q^t * b + FOR_FIRST(i, 0, n, { + + double sum = 0.0; + double sum_lcl = 0.0; + + // Q[m,n] so Q^t[n,m] b[m] + FOR_REDUCE_SUM_SECOND(j, 0, m, + sum_lcl, { + sum_lcl += Q(j,i) * b(j); + }, sum); + y(i) = sum; + + }); // end parallel i + + Kokkos::fence(); + + // Solve R x = y + QR_backsub(R, y, x); +} + +////////////////////////// + + +#endif // QR \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp index 151404f33..fe9fa318d 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_execute.cpp @@ -52,7 +52,8 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, swage::Mesh& mesh, - State_t& State) + State_t& State, + elements::fe_ref_elem_t& ref_elem) { // Conveinent local variables double fuzz = SimulationParamaters.DynamicOptions.fuzz; @@ -63,43 +64,117 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, int graphics_cyc_ival = SimulationParamaters.OutputOptions.graphics_iteration_step; // double time_initial = SimulationParamaters.DynamicOptions.time_initial; - double time_final = this->time_end; //SimulationParamaters.DynamicOptions.time_final; - double dt_min = SimulationParamaters.DynamicOptions.dt_min; - double dt_max = SimulationParamaters.DynamicOptions.dt_max; + double time_final = this->time_end; double dt_start = SimulationParamaters.DynamicOptions.dt_start; - double dt_cfl = SimulationParamaters.DynamicOptions.dt_cfl; - int rk_num_stages = SimulationParamaters.DynamicOptions.rk_num_stages; int cycle_stop = SimulationParamaters.DynamicOptions.cycle_stop; // initialize time, time_step, and cycles - double time_value = this->time_start; // was 0.0 + double time_value = this->time_start; double dt = dt_start; + // defining displacement at start of tlqs solve + FOR_ALL(i, 0, static_cast(mesh.num_nodes), + j, 0, 3, { + State.node.displacement(i,j) = State.node.coords(i,j) - State.node.coords_t0(i,j); + }); + + // ******************************* + // local variables for this solver + // ******************************* + + // element stiffness and force arrays + CArrayKokkos K_elem(mesh.num_elems,3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); /// K1 + K2 + CArrayKokkos F_elem(mesh.num_elems,3*mesh.num_nodes_in_elem); /// F02 - F01 + + // additive schwarz preconditioning variables + /* CArrayKokkos K_elem_inv(mesh.num_elems,3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); + CArrayKokkos intermediate_K_elem_inv(mesh.num_elems,3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); + CArrayKokkos temporary(mesh.num_elems, 3*mesh.num_nodes_in_elem); + CArrayKokkos perm(mesh.num_elems, 3*mesh.num_nodes_in_elem); + CArrayKokkos zk(3 * mesh.num_nodes); + CArrayKokkos zkp1(3 * mesh.num_nodes); */ + + // conjugate gradient method vectors + CArrayKokkos p(3*mesh.num_nodes); + CArrayKokkos rk(3*mesh.num_nodes); + CArrayKokkos rkp1(3*mesh.num_nodes); + + // Anderson acceleration variables + size_t window_size = 5; + DCArrayKokkos anderson_weights(window_size); + CArrayKokkos curr_anderson_residual(3*mesh.num_nodes); + CArrayKokkos hist_anderson_residual(3*mesh.num_nodes,window_size); + hist_anderson_residual.set_values(0); + CArrayKokkos hist_displacement_iter(3*mesh.num_nodes,window_size); + hist_displacement_iter.set_values(0); + + // Picard iteration vectors + // + // displacement_step : current best estimate of the total displacement increment + // for this load step. Passed to get_gradients so that K/F + // are evaluated at the current linearisation point. + // Set (not accumulated) to displacement_iter_kp1 after each + // Anderson-accelerated Picard iterate. + // + // displacement_iter_k : Picard iterate coming *into* the current iteration (x_k). + // Initialised to zero at the start of every load step; + // updated to displacement_iter_kp1 at the end of each iter. + // + // displacement_iter_kp1: result of the CG solve for this Picard iteration G(x_k), + // then overwritten with the Anderson-accelerated update x_{k+1}. + // Reset to zero before every CG solve. + CArrayKokkos displacement_step(3*mesh.num_nodes); /// current load-step displacement estimate + CArrayKokkos displacement_iter_k(3*mesh.num_nodes); /// x_k (Picard iterate in) + CArrayKokkos displacement_iter_kp1(3*mesh.num_nodes); /// G(x_k) then x_{k+1} + + // variables for chebyshev smoothing + CArrayKokkos D_inv(3 * mesh.num_nodes); + CArrayKokkos zk(3 * mesh.num_nodes); + CArrayKokkos zkp1(3 * mesh.num_nodes); + CArrayKokkos delta_z(3 * mesh.num_nodes); + CArrayKokkos temporary(3 * mesh.num_nodes); + + + // Algebraic Multigrid variables + /* int max_layers = 10; + int dof_cutoff = 3000; // NOTE: LOOK INTO THIS VALUE MORE + double theta = 0.25; + double omega = 2/3; + CArrayKokkos A1; + CArrayKokkos A2; + CArrayKokkos A3; + CArrayKokkos A4; + CArrayKokkos A5; + CArrayKokkos A6; + CArrayKokkos A7; + CArrayKokkos A8; + CArrayKokkos A9; + CArrayKokkos A10; + CArrayKokkos P1; + CArrayKokkos P2; + CArrayKokkos P3; + CArrayKokkos P4; + CArrayKokkos P5; + CArrayKokkos P6; + CArrayKokkos P7; + CArrayKokkos P8; + CArrayKokkos P9; + CArrayKokkos P10; */ + // Create mesh writer - MeshWriter mesh_writer; // Note: Pull to driver after refactoring evolution + MeshWriter mesh_writer; // --- Graphics vars ---- CArray graphics_times = CArray(20000); - graphics_times(0) = this->time_start; // was zero - double graphics_time = this->time_start; // the times for writing graphics dump, was started at 0.0 - - - /// WARNING WARNING WARNING: REMOVE BEFORE BUILDING SOLVER, THIS IS A PLACEHOLDER FOR THE TLQS SOLVER - std::cout << "Sucessfully called the TLQS execute function" << std::endl; - return; - - - // Apply initial boundary conditions - + graphics_times(0) = this->time_start; + double graphics_time = this->time_start; double cached_pregraphics_dt = fuzz; // the number of materials specified by the user input const size_t num_mats = Materials.num_mats; - - // a flag to exit the calculation size_t stop_calc = 0; @@ -107,23 +182,26 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // Write initial state at t=0 printf("Writing outputs to file at %f \n", graphics_time); - mesh_writer.write_mesh( - mesh, - State, - SimulationParamaters, - dt, - time_value, - graphics_times, - TLQS3D_State::required_node_state, - TLQS3D_State::required_gauss_pt_state, - TLQS3D_State::required_material_pt_state, - this->solver_id); + mesh_writer.write_mesh_Pn(mesh, + State, + SimulationParamaters, + dt, + time_value, + graphics_times, + TLQS3D_State::required_node_state, + TLQS3D_State::required_gauss_pt_state, + TLQS3D_State::required_material_pt_state, + this->solver_id, + ref_elem); - - graphics_time = time_value + graphics_dt_ival; - // loop over the max number of time integration cycles + // ****************************************************************************************************** + // setting max_iter, need to figure out either a good general number or make it a user set with a default + // ****************************************************************************************************** + int max_iter = 5000; + + // loop over the max number of load steps for (size_t cycle = 0; cycle < cycle_stop; cycle++) { // stop calculation if flag if (stop_calc == 1) { @@ -131,49 +209,456 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, } cached_pregraphics_dt = dt; - - // the smallest time step across all materials - double min_dt_calc = dt_max; - - // calculating time step per material - for(size_t mat_id = 0; mat_id < num_mats; mat_id++){ - - // initialize the material dt - double dt_mat = dt; - - // Initialize timestep - timestep_init(State.node.coords, - State.node.coords_n0, - State.node.vel, - State.node.vel_n0, - State.MaterialPoints.sie, - State.MaterialPoints.sie_n0, - State.MaterialPoints.stress, - State.MaterialPoints.stress_n0, - mesh.num_dims, - mesh.num_elems, - mesh.num_nodes, - mat_id); - - - } // end for loop over all mats - // Print the initial time step and time value if (cycle == 0) { printf("cycle = %lu, time = %f, time step = %f \n", cycle, time_value, dt); } - // print time step every 10 cycles + // print time step every 20 cycles else if (cycle % 20 == 0) { printf("cycle = %lu, time = %f, time step = %f \n", cycle, time_value, dt); - } // end if - - - // --------------------------------------------------------------------- - // integrate the solution forward - // --------------------------------------------------------------------- - - + } + + // ----------------------------------------------------------- + // Reset all Picard / Anderson state for this new load step. + // displacement_step is the running estimate of the load-step + // displacement, so it also starts from zero. + // ----------------------------------------------------------- + displacement_step.set_values(0); + displacement_iter_k.set_values(0); + displacement_iter_kp1.set_values(0); + hist_anderson_residual.set_values(0); + hist_displacement_iter.set_values(0); + Kokkos::fence(); + + // start Picard iteration loop + for (int iter = 0; iter < max_iter; iter++) { + //std::cout << "ITER: " << iter << std::endl; + const auto start_time = std::chrono::steady_clock::now(); + + // dirichlet (displacement) type + boundary_displacement(mesh, BoundaryConditions, K_elem, F_elem, displacement_step, dt, time_value, time_start, time_end); + + auto point_A = std::chrono::steady_clock::now(); + auto elapsed_A = std::chrono::duration_cast(point_A - start_time).count(); + //std::cout << "Time elapsed for first bc: " << elapsed_A << " ms\n"; + // *************************************************** + // get element arrays + // *************************************************** + K_elem.set_values(0); + F_elem.set_values(0); + Kokkos::fence(); + + // looping through materials + for (int mat_id = 0; mat_id < num_mats; mat_id++) { + + // parallel loop through elements + FOR_ALL(elem, 0, State.MaterialToMeshMaps.num_mat_elems.host(mat_id), { + + // setting up views and temp memory + const size_t elem_id = State.MaterialToMeshMaps.elem_in_mat_elem(mat_id, elem); + ViewCArrayKokkos nodes_in_curr_elem(&mesh.nodes_in_elem(elem_id,0),mesh.num_nodes_in_elem); + double grad_u[3][3]; + double inv_J[3][3]; + double det_J; + double PK2_curr_config[6]; + double material_matrix[6][6]; + ViewCArrayKokkos curr_K_elem(&K_elem(elem_id,0,0),3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); + ViewCArrayKokkos curr_F_elem(&F_elem(elem_id,0),3*mesh.num_nodes_in_elem); + + // looping through material points + for (int mat_pt = 0; mat_pt < ref_elem.gauss_point_grad_basis.dims(0); mat_pt++) { + // setting up view and getting material matrix + ViewCArrayKokkos curr_grad_basis(&ref_elem.gauss_point_grad_basis(mat_pt,0,0),ref_elem.num_basis, mesh.num_dims); + Materials.MaterialFunctions(mat_id).fill_C_matrix(Materials.strength_global_vars, material_matrix, mat_id); + + // tallying to element array + TLQS3D::get_gradients(material_matrix, nodes_in_curr_elem, State.node.coords_t0, State.node.displacement, displacement_step, curr_grad_basis, grad_u, inv_J, det_J, PK2_curr_config); + double weight = ref_elem.gauss_point_weights(mat_pt)*det_J; + double local_mat_vol_frac = State.MaterialPoints.mat_volfrac(mat_id, elem); + + TLQS3D::tally_elem_arrays(material_matrix, grad_u, inv_J, curr_grad_basis, weight, PK2_curr_config, curr_K_elem, curr_F_elem, local_mat_vol_frac); + } // end mat_pt + + }); // end elem + + Kokkos::fence(); + + } // end mat_id + auto point_B = std::chrono::steady_clock::now(); + auto elapsed_B = std::chrono::duration_cast(point_B - point_A).count(); + //std::cout << "Time elapsed to build element arrays: " << elapsed_B << " ms\n"; + + // *************************************************** + // end element arrays + // *************************************************** + + // *************************************************** + // apply boundary conditions + // *************************************************** + + // neumann (traction) type + + // dirichlet (displacement) type + boundary_displacement(mesh, BoundaryConditions, K_elem, F_elem, displacement_step, dt, time_value, time_start, time_end); + auto point_C = std::chrono::steady_clock::now(); + auto elapsed_C = std::chrono::duration_cast(point_C - point_B).count(); + //std::cout << "Time elapsed for second bc: " << elapsed_C << " ms\n"; + // *************************************************** + // end boundary conditions + // *************************************************** + + /* // get inverse of K_elem arrays for each element + // initialize K_elem_inv to K_elem since LU functions are destructive + FOR_ALL(i, 0, mesh.num_elems, + j, 0, 3*mesh.num_nodes_in_elem, + k, 0, 3*mesh.num_nodes_in_elem,{ + intermediate_K_elem_inv(i,j,k) = K_elem(i,j,k); + }); + + // Regularize: K_e <- K_e + eps * (trace(K_e)/n_local) * I + // This lifts the null space without significantly affecting the deformation modes + FOR_ALL(elem_gid, 0, mesh.num_elems, { + const size_t n_local = 3 * mesh.num_nodes_in_elem; + double trace = 0.0; + for (size_t d = 0; d < n_local; d++) { + trace += intermediate_K_elem_inv(elem_gid, d, d); + } + const double reg = 1e-6 * trace / (double)n_local; + for (size_t d = 0; d < n_local; d++) { + intermediate_K_elem_inv(elem_gid, d, d) += reg; + } + }); + Kokkos::fence(); + + FOR_ALL(i, 0, mesh.num_elems,{ + ViewCArrayKokkos curr_intermediate_K_elem_inv(&intermediate_K_elem_inv(i,0,0),3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); + ViewCArrayKokkos curr_temporary(&temporary(i, 0), 3*mesh.num_nodes_in_elem); + ViewCArrayKokkos curr_perm(&perm(i, 0), 3*mesh.num_nodes_in_elem); + ViewCArrayKokkos curr_K_elem_inv(&K_elem_inv(i,0,0),3*mesh.num_nodes_in_elem,3*mesh.num_nodes_in_elem); + int parity; + LU_decompose(curr_intermediate_K_elem_inv, curr_perm, curr_temporary, parity); + LU_invert(curr_intermediate_K_elem_inv, curr_perm, curr_K_elem_inv, curr_temporary); + }); */ + + // *************************************************** + // begin conjugate gradient solve + // *************************************************** + + // Reset the CG solution vector each Picard iteration + displacement_iter_kp1.set_values(0); + rk.set_values(0); + Kokkos::fence(); + + // getting inverse of the diagonal like diagonal jacobi + get_diagonal_inverse(D_inv, K_elem, mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem); + + // getting spectral bounds for chebyshev smoothing + // We pass zk and temporary here safely because they are currently uninitialized + // scratch space before the CG loop kicks off. + double alpha = 0.0; + double beta = 0.0; + const int cheb_degree = 3; // Choose your Chebyshev polynomial degree (typically 2 to 5) + + get_chebyshev_bounds(alpha, beta, D_inv, K_elem, mesh.num_nodes, + mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, + zk, temporary, 15); // Running 15 power iterations + + // getting r0 = (02F - 01F) - K * displacement_iter_k + get_r0(mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, F_elem, K_elem, displacement_iter_kp1, rk); + + // smoothing with chebyshev polynomial + apply_chebyshev_preconditioner(rk, zk, D_inv, zk, delta_z, temporary, K_elem, + mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, + alpha, beta, cheb_degree); + + // z0 = M_inv * r0, p0 = z0 + //get_z0(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.nodes_in_elem, mesh.elems_in_node, K_elem, rk, zk); + FOR_ALL(i, 0, 3 * mesh.num_nodes, { + //zk(i) = M_inv(i) * rk(i); + p(i) = zk(i); + }); + Kokkos::fence(); + + // start of PCG iteration loop + for (int cgm_iter = 0; cgm_iter < max_iter; cgm_iter++) { + + // r_k^T * z_k + double rktzk = 0.0; + double loc_rktzk = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rktzk, { + loc_rktzk += rk(i) * zk(i); + }, rktzk); + Kokkos::fence(); + + // alpha_k = (r_k^T * z_k) / (p_k^T * K * p_k) + double alpha_k = get_alpha(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.nodes_in_elem, K_elem, rktzk, p); + + // displacement_iter_kp1 = displacement_iter_kp1 + alpha_k * p_k + FOR_ALL(i, 0, 3*mesh.num_nodes, { + displacement_iter_kp1(i) += alpha_k * p(i); + }); + Kokkos::fence(); + + // r_{k+1} = r_k - alpha_k * K * p_k + get_rkp1(mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, K_elem, rk, p, alpha_k, rkp1); + + // smoothing with chebyshev polynomial + apply_chebyshev_preconditioner(rkp1, zkp1, D_inv, zk, delta_z, temporary, K_elem, + mesh.num_nodes, mesh.elems_in_node, mesh.num_nodes_in_elem, mesh.nodes_in_elem, + alpha, beta, cheb_degree); + + // z_{k+1} = M_inv * r_{k+1} + //get_zkp1(mesh.num_nodes, mesh.num_nodes_in_elem, mesh.nodes_in_elem, mesh.elems_in_node, K_elem, rkp1, zkp1); + /* FOR_ALL(i, 0, 3*mesh.num_nodes, { + zkp1(i) = M_inv(i) * rkp1(i); + }); + Kokkos::fence(); */ + + // check convergence on true residual norm + double rkp1trkp1 = 0.0; + double loc_rkp1trkp1 = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rkp1trkp1, { + loc_rkp1trkp1 += rkp1(i) * rkp1(i); + }, rkp1trkp1); + Kokkos::fence(); + double norm = sqrt(rkp1trkp1); + //std::cout << "CGM iter " << cgm_iter << " residual norm: " << norm << "\n"; + if (norm < 1.0/*1E-10*/) { + break; + } + + // r_{k+1}^T * z_{k+1} + double rkp1tzkp1 = 0.0; + double loc_rkp1tzkp1 = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_rkp1tzkp1, { + loc_rkp1tzkp1 += rkp1(i) * zkp1(i); + }, rkp1tzkp1); + Kokkos::fence(); + + // beta_k = (r_{k+1}^T * z_{k+1}) / (r_k^T * z_k) + double beta_k = rkp1tzkp1 / (rktzk + 1e-16); + + // p_{k+1} = z_{k+1} + beta_k * p_k + FOR_ALL(i, 0, 3*mesh.num_nodes, { + p(i) = zkp1(i) + beta_k * p(i); + }); + + // update rk, zk for next iteration + FOR_ALL(i, 0, 3*mesh.num_nodes, { + rk(i) = rkp1(i); + zk(i) = zkp1(i); + }); + Kokkos::fence(); + + } // end PCG iteration loop + auto point_D = std::chrono::steady_clock::now(); + auto elapsed_D = std::chrono::duration_cast(point_D - point_C).count(); + //std::cout << "Time elapsed for CGM: " << elapsed_D << " ms\n"; + // *************************************************** + // end conjugate gradient solve + // displacement_iter_kp1 now holds G(x_k): the Picard + // fixed-point map applied to the current iterate x_k. + // *************************************************** + + // *************************************************** + // Anderson acceleration of Picard solve + // + // Standard AA-m formulation: + // f_k = G(x_k) - x_k (residual of fixed-point eq.) + // x_{k+1} = sum_i( theta_i * G(x_{hist_i}) ) + // where theta minimises ||F * theta|| s.t. sum(theta) = 1, + // and F = [f_{oldest} ... f_k] (columns, m_anderson of them). + // + // The history is stored in a circular buffer of depth window_size. + // Column overwritten each iteration = iter % window_size. + // A temporary copy of the residual history is passed to qr_solve() + // because QR factorisation is destructive. + // *************************************************** + + // --- Step 1: compute Anderson residual f_k = G(x_k) - x_k --- + FOR_ALL(i, 0, 3*mesh.num_nodes, { + curr_anderson_residual(i) = displacement_iter_kp1(i) - displacement_iter_k(i); + }); + Kokkos::fence(); + + // --- Step 2: store f_k and G(x_k) in circular buffer --- + // Column to overwrite is determined by the modulo of the current + // Picard iteration index against the window size. + int anderson_col = iter % (int)window_size; + + FOR_ALL(i, 0, 3*mesh.num_nodes, { + hist_anderson_residual(i, anderson_col) = curr_anderson_residual(i); + hist_displacement_iter(i, anderson_col) = displacement_iter_kp1(i); + }); + Kokkos::fence(); + + // --- Step 3: determine how many history columns to use --- + // On early iterations the buffer is only partially filled. + int m_anderson = std::min(iter + 1, (int)window_size); + + // --- Step 4: build a temporary copy of the residual history --- + // The copy is ordered oldest -> newest so column 0 of temp is the + // oldest retained residual and column m_anderson-1 is f_k. + // + // Mapping from window slot w (0 = oldest, m_anderson-1 = newest) + // to circular buffer column: + // if buffer not yet full -> hist_col = w + // if buffer full -> hist_col = (anderson_col - m_anderson + 1 + w + // + window_size) % window_size + // Both branches reduce to the same expression when the buffer is not + // yet full (anderson_col == iter, m_anderson == iter+1), but we keep + // them explicit for clarity. + // + // We allocate a fresh CArrayKokkos each iteration so the size matches + // m_anderson exactly. The destructor frees the memory when we leave + // this scope. + CArrayKokkos temp_residual_hist(3*mesh.num_nodes, m_anderson); + + for (int w = 0; w < m_anderson; w++) { + int hist_col; + if (iter + 1 <= (int)window_size) { + // Buffer not yet full: columns 0..iter are valid, in order. + hist_col = w; + } else { + // Buffer full: walk forward from the slot after anderson_col + // (= oldest) to anderson_col itself (= newest). + hist_col = (anderson_col - m_anderson + 1 + w + (int)window_size) % (int)window_size; + } + + FOR_ALL(i, 0, 3*mesh.num_nodes, { + temp_residual_hist(i, w) = hist_anderson_residual(i, hist_col); + }); + Kokkos::fence(); + } + + // --- Step 5: solve for Anderson mixing weights via QR --- + // qr_solve() finds theta such that: + // min_theta || temp_residual_hist * theta ||_2 + // subject to sum(theta_i) = 1 + // and stores the result in anderson_weights(0..m_anderson-1). + // + // The call is destructive with respect to temp_residual_hist; the + // circular buffer hist_anderson_residual is untouched because we + // passed a copy. + anderson_weights.set_values(0); + QR_solver(temp_residual_hist, curr_anderson_residual, anderson_weights); + anderson_weights.update_host(); + + // --- Step 6: form the Anderson-accelerated iterate --- + // x_{k+1} = sum_{w=0}^{m_anderson-1} theta_w * G(x_{hist_col_w}) + // + // anderson_weights(w) are assumed to be accessible on the host after + // qr_solve() returns (sync/copy is the responsibility of qr_solve). + displacement_iter_kp1.set_values(0); + Kokkos::fence(); + + for (int w = 0; w < m_anderson; w++) { + int hist_col; + if (iter + 1 <= (int)window_size) { + hist_col = w; + } else { + hist_col = (anderson_col - m_anderson + 1 + w + (int)window_size) % (int)window_size; + } + + // Read weight on host; capture by value into the kernel. + // (Requires qr_solve to have made anderson_weights host-accessible.) + double theta_w = anderson_weights.host(w); + + FOR_ALL(i, 0, 3*mesh.num_nodes, { + displacement_iter_kp1(i) += theta_w * hist_displacement_iter(i, hist_col); + }); + Kokkos::fence(); + } + auto point_E = std::chrono::steady_clock::now(); + auto elapsed_E = std::chrono::duration_cast(point_E - point_D).count(); + //std::cout << "Time elapsed for anderson: " << elapsed_E << " ms\n"; + // *************************************************** + // end Anderson acceleration + // *************************************************** + + // update displacement step vector for convergence check and next iteration + FOR_ALL(i, 0, 3*mesh.num_nodes, { + displacement_step(i) += displacement_iter_kp1(i); + }); + Kokkos::fence(); + + // convergence check + /* double norm_num = 0.0; + double loc_norm_num = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_norm_num, { + loc_norm_num += curr_anderson_residual(i) * curr_anderson_residual(i); + }, norm_num); + + double norm_den = 0.0; + double loc_norm_den = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_norm_den, { + loc_norm_den += displacement_step(i) * displacement_step(i); + }, norm_den); + + double norm = sqrt(norm_num / norm_den); */ + double norm = 0.0; + double loc_norm = 0.0; + FOR_REDUCE_SUM(i, 0, 3*mesh.num_nodes, loc_norm, { + loc_norm += curr_anderson_residual(i) * curr_anderson_residual(i); + }, norm); + + std::cout << "ITER: " << iter << " ANDERSON RESIDUAL NORM: " << norm << std::endl; + if (norm < 1E-12 && iter > 1) { + std::cout << "PICARD CONVERGED AT ITER: " << iter+1 << std::endl; + break; + } + + // Update x_k <- x_{k+1} for the next Picard iteration. + FOR_ALL(i, 0, 3*mesh.num_nodes, { + displacement_iter_k(i) = displacement_iter_kp1(i); + }); + Kokkos::fence(); + auto point_F = std::chrono::steady_clock::now(); + auto elapsed_F = std::chrono::duration_cast(point_F - point_E).count(); + //std::cout << "Time elapsed for picard convergence and update: " << elapsed_F << " ms\n"; + auto elapsed_iter = std::chrono::duration_cast(point_F - start_time).count(); + //std::cout << "Time elapsed for this full iteration: " << elapsed_iter << " ms\n"; + } // end Picard iteration loop + + // updating total displacement for next load step + FOR_ALL(i, 0, (int)mesh.num_nodes, + j, 0, 3, { + State.node.displacement(i,j) += displacement_step(3*i + j); + State.node.coords(i,j) = State.node.coords_t0(i,j) + State.node.displacement(i,j); + }); + Kokkos::fence(); + + // filling in stress and strain for output + for (int mat_id = 0; mat_id < num_mats; mat_id++) { + + FOR_ALL(elem, 0, State.MaterialToMeshMaps.num_mat_elems.host(mat_id), { + + // setting up views and temp memory + const size_t elem_id = State.MaterialToMeshMaps.elem_in_mat_elem(mat_id, elem); + + ViewCArrayKokkos nodes_in_curr_elem(&mesh.nodes_in_elem(elem_id,0),mesh.num_nodes_in_elem); + double material_matrix[6][6]; + + // looping through material points + for (int mat_pt = 0; mat_pt < ref_elem.gauss_point_grad_basis.dims(0); mat_pt++) { + // setting up view and getting material matrix + ViewCArrayKokkos curr_grad_basis(&ref_elem.gauss_point_grad_basis(mat_pt,0,0),ref_elem.num_basis, mesh.num_dims); + Materials.MaterialFunctions(mat_id).fill_C_matrix(Materials.strength_global_vars, material_matrix, mat_id); + + // views into stress and strain + ViewCArrayKokkos stress_view(&State.MaterialPoints.stress(mat_id, State.points_in_mat_elem(elem,mat_pt),0,0),3,3); + ViewCArrayKokkos strain_view(&State.MaterialPoints.strain(mat_id, State.points_in_mat_elem(elem,mat_pt),0,0),3,3); + + // tallying to element array + TLQS3D::post_process(material_matrix, nodes_in_curr_elem, State.node.coords_t0, State.node.displacement, curr_grad_basis, stress_view, strain_view); + + } // end mat_pt + + }); // end elem + Kokkos::fence(); + + } // end mat_id // increment the time time_value += dt; @@ -196,7 +681,7 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, // write outputs if (write == 1) { printf("Writing outputs to file at %f \n", graphics_time); - mesh_writer.write_mesh(mesh, + mesh_writer.write_mesh_Pn(mesh, State, SimulationParamaters, dt, @@ -205,7 +690,8 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, TLQS3D_State::required_node_state, TLQS3D_State::required_gauss_pt_state, TLQS3D_State::required_material_pt_state, - this->solver_id); + this->solver_id, + ref_elem); graphics_time = time_value + graphics_dt_ival; @@ -223,6 +709,4 @@ void TLQS3D::execute(SimulationParameters_t& SimulationParamaters, printf("\nCalculation time in seconds: %f \n", calc_time * 1e-9); -} // end of SGH execute - - +} // end of TLQS3D execute \ No newline at end of file diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_initialize.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_initialize.cpp index 8b0d48639..3cf30edd4 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_initialize.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_initialize.cpp @@ -44,7 +44,7 @@ void TLQS3D::initialize(SimulationParameters_t& SimulationParamaters, State_t& State) const { const size_t num_nodes = mesh.num_nodes; - const size_t num_gauss_pts = mesh.num_elems; + const size_t num_gauss_pts = mesh.num_gauss_in_elem*mesh.num_elems; const size_t num_corners = mesh.num_corners; const size_t num_dims = mesh.num_dims; @@ -59,15 +59,15 @@ void TLQS3D::initialize(SimulationParameters_t& SimulationParamaters, State.corner.initialize(num_corners, num_dims, TLQS3D_State::required_corner_state); // check that the fills specify the required nodal fields - bool filled_nodal_state = + /* bool filled_nodal_state = check_fill_node_states(TLQS3D_State::required_fill_node_state, - SimulationParamaters.InitialConditionSetup.fill_node_states); + SimulationParamaters.region_setups.fill_node_states); */ - if (filled_nodal_state == false){ + /* if (filled_nodal_state == false){ std::cout <<" Missing required nodal state in the fill instructions for the TLQS solver \n"; - std::cout <<" The nodal velocity must be specified. \n" << std::endl; + std::cout <<" The nodal displacement must be specified. \n" << std::endl; throw std::runtime_error("**** Provide fill instructions for all required nodal variables ****"); - } + } */ std::cout << "TLQS solver initialized \n"; diff --git a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_setup.cpp b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_setup.cpp index 2a3d21b4d..c2434d234 100644 --- a/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_setup.cpp +++ b/apps/multiphysics/src/Solvers/TLQS_solver_3D/src/tlqs_setup.cpp @@ -60,6 +60,19 @@ void TLQS3D::setup(SimulationParameters_t& SimulationParamaters, const size_t num_mats = Materials.num_mats; // the number of materials on the mesh + // calculate pressure, sound speed, and stress for each material + for (int mat_id = 0; mat_id < num_mats; mat_id++) { + + // call the initialization function for state vars + init_state_vars(Materials, + mesh, + State.MaterialPoints.eos_state_vars, + State.MaterialPoints.strength_state_vars, + State.MaterialToMeshMaps.elem_in_mat_elem, + State.MaterialPoints.num_material_points.host(mat_id), + mat_id); + + } // for loop over mat_id std::cout << "TLQS solver setup \n"; return; diff --git a/apps/multiphysics/src/Solvers/level_set_solver/include/level_set_solver.hpp b/apps/multiphysics/src/Solvers/level_set_solver/include/level_set_solver.hpp index 6fff7601e..a68512add 100644 --- a/apps/multiphysics/src/Solvers/level_set_solver/include/level_set_solver.hpp +++ b/apps/multiphysics/src/Solvers/level_set_solver/include/level_set_solver.hpp @@ -149,7 +149,8 @@ class LevelSet : public Solver Material_t& Materials, BoundaryCondition_t& Boundary, swage::Mesh& mesh, - State_t& State) override; + State_t& State, + elements::fe_ref_elem_t& ref_elem) override; ///////////////////////////////////////////////////////////////////////////// /// diff --git a/apps/multiphysics/src/Solvers/level_set_solver/src/level_set_execute.cpp b/apps/multiphysics/src/Solvers/level_set_solver/src/level_set_execute.cpp index acda5f0cc..0ef650dd3 100644 --- a/apps/multiphysics/src/Solvers/level_set_solver/src/level_set_execute.cpp +++ b/apps/multiphysics/src/Solvers/level_set_solver/src/level_set_execute.cpp @@ -53,7 +53,8 @@ void LevelSet::execute(SimulationParameters_t& SimulationParamaters, Material_t& Materials, BoundaryCondition_t& BoundaryConditions, swage::Mesh& mesh, - State_t& State) + State_t& State, + elements::fe_ref_elem_t& ref_elem) { // arrays local to this solver diff --git a/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp new file mode 100644 index 000000000..d9441399d --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/cyclic_displacement_bc.hpp @@ -0,0 +1,122 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_CYCLIC_H +#define BOUNDARY_DISP_CYCLIC_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace CyclicDisplacementBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn displacement +/// +/// \brief This is a function to set the displacement along a symmetry plane or +/// a wall. This fcn imposes a normal displacement, in the specified +/// direction, to be equal to zero +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + // disp bc global vars col 0: amplitude + // disp bc global vars col 1: frequency in Hz + // disp bc global vars col 2: direction + + const size_t constrained_dir = static_cast(disp_bc_global_vars(bdy_set, 2)); + const size_t num_elems_in_node = mesh.elems_in_node.stride(bdy_node_gid); + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const size_t num_dof_in_elem = 3 * num_nodes_in_elem; + + double ang_freq = 2*3.14159265358979323846*disp_bc_global_vars(bdy_set, 1); + double disp_curr = disp_bc_global_vars(bdy_set, 0)*sin((time_value-time_start)*ang_freq); + double disp_next = disp_bc_global_vars(bdy_set, 0)*sin((time_value+dt-time_start)*ang_freq); + const double total_step_target = disp_next - disp_curr; + + displacement_step(3 * bdy_node_gid + constrained_dir) = total_step_target; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); + + // Find the local node index for bdy_node_gid within this element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (mesh.nodes_in_elem(elem_gid, a) == bdy_node_gid) { + local_node_lid = a; + break; + } + } + + const size_t constrained_dof = 3 * local_node_lid + constrained_dir; + + for (size_t col = 0; col < num_dof_in_elem; col++) + K_elem(elem_gid, constrained_dof, col) = 0.0; + + for (size_t row = 0; row < num_dof_in_elem; row++) + K_elem(elem_gid, row, constrained_dof) = 0.0; + + K_elem(elem_gid, constrained_dof, constrained_dof) = 1.0 / (double)num_elems_in_node; + F_elem(elem_gid, constrained_dof) = 0.0; + } + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/fixed_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/fixed_displacement_bc.hpp new file mode 100644 index 000000000..42442aaeb --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/fixed_displacement_bc.hpp @@ -0,0 +1,112 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_FIXED_H +#define BOUNDARY_DISP_FIXED_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary displacement is set to zero in all directions +/// +/// \brief This is a function to set the displacement all directions to a value +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +namespace FixedDisplacementBC +{ +// add an enum for boundary statevars and global vars + +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + const size_t num_elems_in_node = mesh.elems_in_node.stride(bdy_node_gid); + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const size_t num_dof_in_elem = 3 * num_nodes_in_elem; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); + + // Find the local node index for bdy_node_gid within this element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (mesh.nodes_in_elem(elem_gid, a) == bdy_node_gid) { + local_node_lid = a; + break; + } + } + + // Zero rows and columns for all 3 constrained DOFs of this node + for (size_t p = 0; p < 3; p++) { + const size_t constrained_dof = 3 * local_node_lid + p; + + for (size_t col = 0; col < num_dof_in_elem; col++) + K_elem(elem_gid, constrained_dof, col) = 0.0; + + for (size_t row = 0; row < num_dof_in_elem; row++) + K_elem(elem_gid, row, constrained_dof) = 0.0; + + K_elem(elem_gid, constrained_dof, constrained_dof) = 1.0 / (double)num_elems_in_node; + F_elem(elem_gid, constrained_dof) = 0; + } + } + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/no_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/no_displacement_bc.hpp new file mode 100644 index 000000000..a192c9729 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/no_displacement_bc.hpp @@ -0,0 +1,80 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_NONE_H +#define BOUNDARY_DISP_NONE_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace NoDisplacementBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary displacement is never set +/// +/// \brief This is a function to set the displacement all directions to a value +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + // this is a blank function by design + return; +} // end displacement +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp new file mode 100644 index 000000000..0589351c4 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/piston_displacement_bc.hpp @@ -0,0 +1,119 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_PISTON_H +#define BOUNDARY_DISP_PISTON_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace PistonDisplacementBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn displacement +/// +/// \brief This is a function to set the displacement along a symmetry plane or +/// a wall. This fcn imposes a normal displacement, in the specified +/// direction, to be equal to zero +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + // disp bc global vars col 0: final u + // disp bc global vars col 1: direction + + const size_t constrained_dir = static_cast(disp_bc_global_vars(bdy_set, 1)); + const size_t num_elems_in_node = mesh.elems_in_node.stride(bdy_node_gid); + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const size_t num_dof_in_elem = 3 * num_nodes_in_elem; + + const double total_step_target = disp_bc_global_vars(bdy_set, 0) * (dt / (time_end - time_start)); + + displacement_step(3 * bdy_node_gid + constrained_dir) = total_step_target; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); + + // Find the local node index for bdy_node_gid within this element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (mesh.nodes_in_elem(elem_gid, a) == bdy_node_gid) { + local_node_lid = a; + break; + } + } + + const size_t constrained_dof = 3 * local_node_lid + constrained_dir; + + for (size_t col = 0; col < num_dof_in_elem; col++) + K_elem(elem_gid, constrained_dof, col) = 0.0; + + for (size_t row = 0; row < num_dof_in_elem; row++) + K_elem(elem_gid, row, constrained_dof) = 0.0; + + K_elem(elem_gid, constrained_dof, constrained_dof) = 1.0 / (double)num_elems_in_node; + F_elem(elem_gid, constrained_dof) = 0.0; + + } + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/reflected_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/reflected_displacement_bc.hpp new file mode 100644 index 000000000..eb7fde097 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/reflected_displacement_bc.hpp @@ -0,0 +1,113 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_REFLECTED_H +#define BOUNDARY_DISP_REFLECTED_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace ReflectedDisplacementBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn displacement +/// +/// \brief This is a function to set the displacement along a symmetry plane or +/// a wall. This fcn imposes a normal displacement, in the specified +/// direction, to be equal to zero +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + // disp bc global vars col 0: direction + + const size_t constrained_dir = static_cast(disp_bc_global_vars(bdy_set, 0)); + const size_t num_elems_in_node = mesh.elems_in_node.stride(bdy_node_gid); + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const size_t num_dof_in_elem = 3 * num_nodes_in_elem; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); + + // Find the local node index for bdy_node_gid within this element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (mesh.nodes_in_elem(elem_gid, a) == bdy_node_gid) { + local_node_lid = a; + break; + } + } + + const size_t constrained_dof = 3 * local_node_lid + constrained_dir; + + for (size_t col = 0; col < num_dof_in_elem; col++) + K_elem(elem_gid, constrained_dof, col) = 0.0; + + for (size_t row = 0; row < num_dof_in_elem; row++) + K_elem(elem_gid, row, constrained_dof) = 0.0; + + K_elem(elem_gid, constrained_dof, constrained_dof) = 1.0 / (double)num_elems_in_node; + F_elem(elem_gid, constrained_dof) = 0; + } + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp new file mode 100644 index 000000000..40884cdf9 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/total_displacement_bc.hpp @@ -0,0 +1,124 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_TOTAL_H +#define BOUNDARY_DISP_TOTAL_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary displacement is set to zero in all directions +/// +/// \brief This is a function to set the displacement all directions to a value +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +namespace TotalDisplacementBC +{ +// add an enum for boundary statevars and global vars + +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + // disp bc global vars col 0: final ux + // disp bc global vars col 1: final uy + // disp bc global vars col 2: final uz + + const size_t num_elems_in_node = mesh.elems_in_node.stride(bdy_node_gid); + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const size_t num_dof_in_elem = 3 * num_nodes_in_elem; + + const double total_step_target_0 = disp_bc_global_vars(bdy_set, 0) * (dt / (time_end - time_start)); + const double total_step_target_1 = disp_bc_global_vars(bdy_set, 1) * (dt / (time_end - time_start)); + const double total_step_target_2 = disp_bc_global_vars(bdy_set, 2) * (dt / (time_end - time_start)); + + displacement_step(3 * bdy_node_gid) = total_step_target_0; + displacement_step(3 * bdy_node_gid + 1) = total_step_target_1; + displacement_step(3 * bdy_node_gid + 2) = total_step_target_2; + + for (size_t elem_lid = 0; elem_lid < num_elems_in_node; elem_lid++) { + const size_t elem_gid = mesh.elems_in_node(bdy_node_gid, elem_lid); + + // Find the local node index for bdy_node_gid within this element + size_t local_node_lid = num_nodes_in_elem; // sentinel + for (size_t a = 0; a < num_nodes_in_elem; a++) { + if (mesh.nodes_in_elem(elem_gid, a) == bdy_node_gid) { + local_node_lid = a; + break; + } + } + + // Zero rows and columns for all 3 constrained DOFs of this node + for (size_t p = 0; p < 3; p++) { + const size_t constrained_dof = 3 * local_node_lid + p; + + for (size_t col = 0; col < num_dof_in_elem; col++) + K_elem(elem_gid, constrained_dof, col) = 0.0; + + for (size_t row = 0; row < num_dof_in_elem; row++) + K_elem(elem_gid, row, constrained_dof) = 0.0; + + K_elem(elem_gid, constrained_dof, constrained_dof) = 1.0 / (double)num_elems_in_node; + F_elem(elem_gid, constrained_dof) = 0.0; + } + } + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/displacement/user_defined_displacement_bc.hpp b/apps/multiphysics/src/boundary_conditions/displacement/user_defined_displacement_bc.hpp new file mode 100644 index 000000000..f5dde6946 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/displacement/user_defined_displacement_bc.hpp @@ -0,0 +1,83 @@ +/********************************************************************************************** +� 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 BOUNDARY_DISP_USER_DEFINED_H +#define BOUNDARY_DISP_USER_DEFINED_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary displacement is user defined +/// +/// \brief This is a function to set the displacement based on user implementation +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node displacement +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +namespace UserDefinedDisplacementBC +{ +// add an enum for boundary statevars and global vars + +KOKKOS_FUNCTION +static void displacement(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + // add user coding here + + return; +} // end func +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/quasi_static_stress/cyclic_quasi_static_stress_bc.hpp b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/cyclic_quasi_static_stress_bc.hpp new file mode 100644 index 000000000..804a196f5 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/cyclic_quasi_static_stress_bc.hpp @@ -0,0 +1,77 @@ +/********************************************************************************************** +� 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 BOUNDARY_QSTATX_STRESS_CYCLIC_H +#define BOUNDARY_QSTATX_STRESS_CYCLIC_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace CyclicQstatxStressBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary stress does not exist, its a free surface +/// +/// \brief This is a function for a free surface, the default case +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node force +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void qstatx_stress(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + + return; +} // end qstatx stress +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/quasi_static_stress/linear_quasi_static_stress_bc.hpp b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/linear_quasi_static_stress_bc.hpp new file mode 100644 index 000000000..c9cd34db0 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/linear_quasi_static_stress_bc.hpp @@ -0,0 +1,77 @@ +/********************************************************************************************** +� 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 BOUNDARY_QSTATX_STRESS_LINEAR_H +#define BOUNDARY_QSTATX_STRESS_LINEAR_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace LinearQstatxStressBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary stress does not exist, its a free surface +/// +/// \brief This is a function for a free surface, the default case +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node force +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void qstatx_stress(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + + return; +} // end qstatx stress +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/quasi_static_stress/no_quasi_static_stress_bc.hpp b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/no_quasi_static_stress_bc.hpp new file mode 100644 index 000000000..1240160fd --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/no_quasi_static_stress_bc.hpp @@ -0,0 +1,77 @@ +/********************************************************************************************** +� 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 BOUNDARY_QSTATX_STRESS_NONE_H +#define BOUNDARY_QSTATX_STRESS_NONE_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace NoQstatxStressBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary stress does not exist, its a free surface +/// +/// \brief This is a function for a free surface, the default case +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node force +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void qstatx_stress(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + + return; +} // end qstatx stress +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/quasi_static_stress/uniform_quasi_static_stress_bc.hpp b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/uniform_quasi_static_stress_bc.hpp new file mode 100644 index 000000000..acb00b82a --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/uniform_quasi_static_stress_bc.hpp @@ -0,0 +1,77 @@ +/********************************************************************************************** +� 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 BOUNDARY_QSTATX_STRESS_UNIFORM_H +#define BOUNDARY_QSTATX_STRESS_UNIFORM_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace UniformQstatxStressBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary stress does not exist, its a free surface +/// +/// \brief This is a function for a free surface, the default case +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node force +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void qstatx_stress(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + + return; +} // end qstatx stress +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/boundary_conditions/quasi_static_stress/user_defined_quasi_static_stress_bc.hpp b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/user_defined_quasi_static_stress_bc.hpp new file mode 100644 index 000000000..349cfb4b8 --- /dev/null +++ b/apps/multiphysics/src/boundary_conditions/quasi_static_stress/user_defined_quasi_static_stress_bc.hpp @@ -0,0 +1,77 @@ +/********************************************************************************************** +� 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 BOUNDARY_QSTATX_STRESS_USER_DEFINED_H +#define BOUNDARY_QSTATX_STRESS_USER_DEFINED_H + +#include "boundary_conditions.hpp" + +struct BoundaryConditionEnums_t; + +namespace UserDefinedQstatxStressBC +{ +///////////////////////////////////////////////////////////////////////////// +/// +/// \fn Boundary stress does not exist, its a free surface +/// +/// \brief This is a function for a free surface, the default case +/// +/// \param Mesh object +/// \param Boundary condition enums to select options +/// \param Boundary condition global variables array +/// \param Boundary condition state variables array +/// \param Node force +/// \param Time of the simulation +/// \param Boundary global index for the surface node +/// \param Boundary set local id +/// +///////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +static void qstatx_stress(const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) +{ + + + return; +} // end qstatx stress +} // end namespace + +#endif // end Header Guard \ No newline at end of file diff --git a/apps/multiphysics/src/common/include/boundary_conditions.hpp b/apps/multiphysics/src/common/include/boundary_conditions.hpp index ed46d3030..d3939d6d2 100644 --- a/apps/multiphysics/src/common/include/boundary_conditions.hpp +++ b/apps/multiphysics/src/common/include/boundary_conditions.hpp @@ -66,7 +66,20 @@ enum BCVelocityModels reflectedVelocityBC = 3, zeroVelocityBC = 4, userDefinedVelocityBC = 5, - pistonVelocityBC = 6 + pistonVelocityBC = 6, + rollerVelocityBC = 7 +}; + +enum BCDisplacementModels +{ + noDisplacementBC = 0, + totalDisplacementBC = 1, + reflectedDisplacementBC = 2, + fixedDisplacementBC = 3, + userDefinedDisplacementBC = 4, + pistonDisplacementBC = 5, + rollerDisplacementBC = 6, + cyclicDisplacementBC = 7 }; // types of temperature boundary conditions @@ -91,12 +104,15 @@ enum BCStressModels preloadContact = 5, fractureStressBC = 6 }; -// future model options: -// displacementBC -// temperatureBC -// contactBC - +enum BCQstatxStressModels +{ + noQstatxStressBC = 0, + uniformQstatxStressBC = 1, + linearQstatxStressBC = 2, + cyclicQstatxStressBC = 3, + userDefinedQstatxStressBC = 4, +}; enum BCFcnLocation { @@ -129,7 +145,22 @@ static std::map bc_velocity_ { "reflected", boundary_conditions::reflectedVelocityBC }, { "fixed", boundary_conditions::zeroVelocityBC }, { "user_defined", boundary_conditions::userDefinedVelocityBC }, - { "piston", boundary_conditions::pistonVelocityBC } + { "piston", boundary_conditions::pistonVelocityBC }, + { "roller", boundary_conditions::rollerVelocityBC } +}; + + +// Displacement models +static std::map bc_displacement_model_map +{ + { "none", boundary_conditions::noDisplacementBC }, + { "total", boundary_conditions::totalDisplacementBC }, + { "reflected", boundary_conditions::reflectedDisplacementBC }, + { "fixed", boundary_conditions::fixedDisplacementBC }, + { "user_defined", boundary_conditions::userDefinedDisplacementBC }, + { "piston", boundary_conditions::pistonDisplacementBC }, + { "roller", boundary_conditions::rollerDisplacementBC }, + { "cyclic", boundary_conditions::cyclicDisplacementBC } }; @@ -156,6 +187,14 @@ static std::map bc_stress_mode { "fracture", boundary_conditions::fractureStressBC } }; +static std::map bc_qstatx_stress_model_map +{ + { "none", boundary_conditions::noQstatxStressBC }, + { "uniform", boundary_conditions::uniformQstatxStressBC }, + { "linear", boundary_conditions::linearQstatxStressBC }, + { "cyclic", boundary_conditions::cyclicQstatxStressBC }, + { "user_defined", boundary_conditions::userDefinedQstatxStressBC }, +}; static std::map bc_location_map @@ -193,11 +232,17 @@ struct BoundaryConditionEnums_t // BC model for velocity boundary_conditions::BCVelocityModels BCVelocityModel = boundary_conditions::noVelocityBC; ///< Type of velocity boundary condition + // BC model for Displacement + boundary_conditions::BCDisplacementModels BCDisplacementModel = boundary_conditions::noDisplacementBC; ///< Type of displacement boundary condition + // BC model for temperature boundary_conditions::BCTemperatureModels BCTemperatureModel = boundary_conditions::noTemperatureBC; ///< Type of temperature boundary condition // BC model for stress boundary_conditions::BCStressModels BCStressModel = boundary_conditions::noStressBC; ///< Type of stress boundary condition + + // BC model for quasi static stress + boundary_conditions::BCQstatxStressModels BCQstatxStressModel = boundary_conditions::noQstatxStressBC; ///< Type of stress boundary condition boundary_conditions::BCFcnLocation Location = boundary_conditions::device; // host or device BC function }; // end boundary condition enums @@ -222,6 +267,21 @@ struct BoundaryConditionFunctions_t const size_t bdy_node_gid, const size_t bdy_set) = NULL; + // function pointer for displacement BC's + void (*displacement) (const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& disp_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const CArrayKokkos& K_elem, + const CArrayKokkos& F_elem, + const CArrayKokkos& displacement_step, + const double dt, + const double time_value, + const double time_start, + const double time_end, + const size_t bdy_node_gid, + const size_t bdy_set) = NULL; + // function pointer for temperature BC's void (*temperature) (const swage::Mesh& mesh, const DCArrayKokkos& BoundaryConditionEnums, @@ -248,7 +308,7 @@ struct BoundaryConditionFunctions_t // function pointer for stress BC's void (*stress) (const swage::Mesh& mesh, const DCArrayKokkos& BoundaryConditionEnums, - const RaggedRightArrayKokkos& vel_bc_global_vars, + const RaggedRightArrayKokkos& stress_bc_global_vars, const DCArrayKokkos& bc_state_vars, const ViewCArrayKokkos & corner_surf_force, const ViewCArrayKokkos & corner_surf_normal, @@ -257,6 +317,17 @@ struct BoundaryConditionFunctions_t const size_t bdy_node_gid, const size_t bdy_set) = NULL; + // function pointer for stress BC's + void (*qstatx_stress) (const swage::Mesh& mesh, + const DCArrayKokkos& BoundaryConditionEnums, + const RaggedRightArrayKokkos& qstatx_stress_bc_global_vars, + const DCArrayKokkos& bc_state_vars, + const ViewCArrayKokkos & corner_surf_force, + const ViewCArrayKokkos & corner_surf_normal, + const double time_value, + const size_t bdy_node_gid, + const size_t bdy_set) = NULL; + }; // end boundary condition fcns @@ -286,9 +357,16 @@ struct BoundaryCondition_t DCArrayKokkos vel_bdy_sets_in_solver; // (solver_id, bc_lid) DCArrayKokkos num_vel_bdy_sets_in_solver; // (solver_id) + // making a psuedo dual ragged right + DCArrayKokkos disp_bdy_sets_in_solver; // (solver_id, bc_lid) + DCArrayKokkos num_disp_bdy_sets_in_solver; // (solver_id) + DCArrayKokkos stress_bdy_sets_in_solver; // (solver_id, bc_lid) DCArrayKokkos num_stress_bdy_sets_in_solver; // (solver_id) + DCArrayKokkos qstatx_stress_bdy_sets_in_solver; // (solver_id, bc_lid) + DCArrayKokkos num_qstatx_stress_bdy_sets_in_solver; // (solver_id) + // keep adding ragged storage for the other BC models -- temp, displacement, etc. DCArrayKokkos temperature_bdy_sets_in_solver; // (solver_id, lids) DCArrayKokkos num_temperature_bdy_sets_in_solver; // (solver_id) @@ -307,12 +385,20 @@ struct BoundaryCondition_t // global variables for velocity boundary condition models RaggedRightArrayKokkos velocity_bc_global_vars; // (bc_id, vars...) - CArrayKokkos num_velocity_bc_global_vars; + CArrayKokkos num_velocity_bc_global_vars; + + // global variables for displacement boundary condition models + RaggedRightArrayKokkos displacement_bc_global_vars; // (bc_id, vars...) + CArrayKokkos num_displacement_bc_global_vars; // global variables for stress boundary condition models RaggedRightArrayKokkos stress_bc_global_vars; // (bc_id, vars...) CArrayKokkos num_stress_bc_global_vars; + // global variables for qstatx stress boundary condition models + RaggedRightArrayKokkos qstatx_stress_bc_global_vars; // (bc_id, vars...) + CArrayKokkos num_qstatx_stress_bc_global_vars; + // global variables for temperature boundary condition models RaggedRightArrayKokkos temperature_bc_global_vars; CArrayKokkos num_temperature_bc_global_vars; @@ -333,8 +419,12 @@ static std::vector str_bc_inps "surface", "velocity_model", "velocity_bc_global_vars", + "displacement_model", + "displacement_bc_global_vars", "stress_model", "stress_bc_global_vars", + "qstatx_stress_model", + "qstatx_stress_bc_global_vars", "contact_bc_vars", "temperature_model", "temperature_bc_global_vars", diff --git a/apps/multiphysics/src/common/include/material.hpp b/apps/multiphysics/src/common/include/material.hpp index c87a4f513..edc396ed6 100644 --- a/apps/multiphysics/src/common/include/material.hpp +++ b/apps/multiphysics/src/common/include/material.hpp @@ -61,6 +61,7 @@ namespace model hostANNStrength = 4, hypoElasticPlasticStrength = 5, hypoElasticPlasticStrengthRZ = 6, + QSIsotropicLinearElastic = 7, }; // EOS model types @@ -173,6 +174,7 @@ static std::map strength_models_map { "hypo_elastic_plastic_strength", model::hypoElasticPlasticStrength }, { "hypo_elastic_plastic_strength_rz", model::hypoElasticPlasticStrengthRZ }, { "host_ann_strength", model::hostANNStrength }, + { "qstatx_isotropic_linear_elastic", model::QSIsotropicLinearElastic}, }; @@ -389,6 +391,10 @@ struct MaterialFunctions_t const size_t num_material_points, const size_t mat_id) = NULL; + void (*fill_C_matrix) ( + const RaggedRightArrayKokkos &strength_global_vars, + double C[6][6], + const size_t mat_id) = NULL; // -- Erosion -- diff --git a/apps/multiphysics/src/common/include/mesh_io.hpp b/apps/multiphysics/src/common/include/mesh_io.hpp index cd172f80c..2d9ed2741 100644 --- a/apps/multiphysics/src/common/include/mesh_io.hpp +++ b/apps/multiphysics/src/common/include/mesh_io.hpp @@ -352,7 +352,8 @@ class MeshReader void read_mesh(swage::Mesh& mesh, MPICArrayKokkos& node_coords, MeshInput_t& mesh_inps, - int num_dims) + int num_dims, + bool HexN) { if (mesh_file_ == NULL) { throw std::runtime_error("**** No mesh path given for read_mesh ****"); @@ -387,10 +388,15 @@ class MeshReader read_Abaqus_mesh(mesh, node_coords, mesh_inps, num_dims); } else if(extension == "vtk"){ // vtk file format - read_vtk_mesh(mesh, node_coords, mesh_inps, num_dims); + if (!HexN) { + read_vtk_mesh(mesh, node_coords, mesh_inps, num_dims); + } + else { + read_vtk_hexN_mesh(mesh, node_coords, mesh_inps, num_dims); + } } else if(extension == "vtu"){ // vtu file format - read_vtu_mesh(mesh, node_coords, mesh_inps, num_dims); + read_vtu_mesh(mesh, node_coords, mesh_inps, num_dims, HexN); } else{ throw std::runtime_error("**** Mesh file extension not understood ****"); @@ -916,6 +922,185 @@ class MeshReader } // end of VTKread function + ///////////////////////////////////////////////////////////////////////////// + /// + /// \fn read_vtk_hexN_mesh + /// + /// \brief Read ASCII .vtk mesh file for arbitrary order hex elements (HexN) + /// using Matar CArrays for internal mapping. + /// + /// \param Simulation mesh + /// \param Simulation state + /// \param Node state struct + /// \param Number of dimensions + /// + ///////////////////////////////////////////////////////////////////////////// + void read_vtk_hexN_mesh(swage::Mesh& mesh, + MPICArrayKokkos& node_coords, + MeshInput_t& mesh_inps, + int num_dims) + { + + std::cout << "Reading VTK HexN mesh" << std::endl; + + int i; // used for writing information to file + int node_gid; // the global id for the point + int elem_gid; // the global id for the elem + + size_t num_nodes_in_elem = 1; + for (int dim = 0; dim < num_dims; dim++) { + num_nodes_in_elem *= (mesh_inps.p_order+1); + } + std::string token; + bool found = false; + + std::ifstream in; + in.open(mesh_file_); // Uses mesh_file_ from the class/scope + + + // --- 1. Find and Read POINTS --- + i = 0; + while (found==false) { + std::string str; + std::string delimiter = " "; + std::getline(in, str); + std::vector v = split (str, delimiter); + + if(v[0] == "POINTS"){ + size_t num_nodes = std::stoi(v[1]); + printf("Number of nodes read in %zu\n", num_nodes); + + mesh.initialize_nodes(num_nodes); + //std::vector required_node_state = { node_state::coords }; + //node.initialize(num_nodes, num_dims, required_node_state); + + found=true; + } + + if (i > 1000){ + std::cerr << "ERROR: Failed to find POINTS in file" << std::endl; + break; + } + i++; + } + + node_coords = MPICArrayKokkos(mesh.num_nodes, num_dims, "Node_coordinates_in_mesh_io"); + + // Read node coordinates and apply user scaling + for (node_gid=0; node_gid < mesh.num_nodes; node_gid++){ + std::string str; + std::getline(in, str); + std::vector v = split (str, " "); + + node_coords.host(node_gid, 0) = mesh_inps.scale_x * std::stod(v[0]); + node_coords.host(node_gid, 1) = mesh_inps.scale_y * std::stod(v[1]); + if(num_dims == 3){ + node_coords.host(node_gid, 2) = mesh_inps.scale_z * std::stod(v[2]); + } + } + node_coords.update_device(); + + + // --- 2. Find and Read CELLS --- + found = false; + i = 0; + size_t num_elem = 0; + while (found==false) { + std::string str; + std::getline(in, str); + std::vector v = split (str, " "); + + if(v[0] == "CELLS"){ + num_elem = std::stoi(v[1]); + printf("Number of elements read in %zu\n", num_elem); + mesh.initialize_elems_Pn(num_elem, num_dims, mesh_inps.p_order); + found = true; + } + + if (i > 1000){ + printf("ERROR: Failed to find CELLS \n"); + break; + } + i++; + } + + // --- 3. Connectivity and Reordering --- + CArray convert_vtk_to_fierro(num_nodes_in_elem); + std::cout << "HEREHEREHERE " << num_nodes_in_elem << std::endl << std::endl; + bool map_built = false; + + for (elem_gid=0; elem_gid < num_elem; elem_gid++) { + std::string str; + std::getline(in, str); + std::vector v = split (str, " "); + num_nodes_in_elem = std::stoi(v[0]); + const int num_1D_points = std::round(std::cbrt(num_nodes_in_elem)); + const int Pn_order = num_1D_points - 1; + + int this_point = 0; + for (int k=0; k <= Pn_order; k++){ + for (int j=0; j <= Pn_order; j++){ + for (int i_idx=0; i_idx <= Pn_order; i_idx++){ + + int order[3] = {Pn_order, Pn_order, Pn_order}; + int this_index = PointIndexFromIJK(i_idx, j, k, order); + + convert_vtk_to_fierro(this_point) = this_index; + this_point++; + } + } + } + + // Map connectivity from VTK to Fierro/Swage mesh structure + for (size_t node_lid=0; node_lid < num_nodes_in_elem; node_lid++){ + int vtk_index = convert_vtk_to_fierro(node_lid); + mesh.nodes_in_elem.host(elem_gid, node_lid) = (size_t)std::stod(v[vtk_index+1]); + } + } + + mesh.nodes_in_elem.update_device(); + + // Initialize corners based on dynamic nodes-per-element count + size_t num_corners = num_elem * num_nodes_in_elem; + mesh.initialize_corners(num_corners); + + // Build connectivity (Faces, etc.) + mesh.build_connectivity(); + + + // --- 4. Validate CELL_TYPES --- + found = false; + i = 0; + size_t elem_type = 0; + while (found==false) { + std::string str; + std::getline(in, str); + std::vector v = split (str, " "); + + if(v[0] == "CELL_TYPES"){ + std::getline(in, str); + elem_type = std::stoi(str); + found = true; + } + + if (i > 1000){ + printf("ERROR: Failed to find CELL_TYPES \n"); + break; + } + i++; + } + + printf("Element type read = %zu \n", elem_type); + + // 12 = Linear Hex, 72 = Lagrange Hex (High Order) + if(elem_type != 12 && elem_type != 72) { + std::cerr << "WARNING: element type " << elem_type << " may not be a supported Hex type." << std::endl; + } + + in.close(); + + } // end of read_vtk_hexN_mesh + ///////////////////////////////////////////////////////////////////////////// /// @@ -932,7 +1117,8 @@ class MeshReader void read_vtu_mesh(swage::Mesh& mesh, MPICArrayKokkos& node_coords, MeshInput_t& mesh_inps, - int num_dims) + int num_dims, + bool HexN) { std::cout<<"Reading VTU file in a multiblock VTK mesh"<& node_coords, - SimulationParameters_t& SimulationParamaters) + SimulationParameters_t& SimulationParamaters, + bool HexN) { if (SimulationParamaters.MeshInput.num_dims == 2) { if (SimulationParamaters.MeshInput.type == mesh_input::Polar) { @@ -1277,7 +1469,12 @@ class MeshBuilder } } else if (SimulationParamaters.MeshInput.num_dims == 3) { - build_3d_box(mesh, node_coords, SimulationParamaters); + if (!HexN) { + build_3d_box(mesh, node_coords, SimulationParamaters); + } + else { + build_3d_HexN_box(mesh, node_coords, SimulationParamaters); + } } else{ throw std::runtime_error("**** ONLY 2D RZ OR 3D MESHES ARE SUPPORTED ****"); @@ -1661,7 +1858,6 @@ class MeshBuilder MPICArrayKokkos& node_coords, SimulationParameters_t& SimulationParamaters) const { - printf(" ***** WARNING:: build_3d_HexN_box not yet implemented\n"); const int num_dim = 3; // SimulationParamaters.MeshInput.length.update_host(); @@ -1749,7 +1945,7 @@ class MeshBuilder // initialize elem variables - mesh.initialize_elems(num_elems, num_dim); + mesh.initialize_elems_Pn(num_elems, num_dim, SimulationParamaters.MeshInput.p_order); // --- Build elems --- @@ -1892,6 +2088,9 @@ class MeshWriter case material_pt_state::stress: State.MaterialPoints.stress.update_host(); break; + case material_pt_state::strain: + State.MaterialPoints.strain.update_host(); + break; // additional vars for thermal-mechanical solver case material_pt_state::thermal_conductivity: @@ -1902,6 +2101,10 @@ class MeshWriter State.MaterialPoints.specific_heat.update_host(); break; + case material_pt_state::heat_flux: + State.MaterialPoints.q_flux.update_host(); + break; + // add other variables here // not used @@ -1911,8 +2114,7 @@ class MeshWriter break; case material_pt_state::poisson_ratios: break; - case material_pt_state::heat_flux: - break; + default: std::cout<<"Desired material point state not understood in outputs"<1){ - for (int k = 0; k <= Pn_order_z; k++) { + CArray convert_fierro_to_vtk((Pn_order+1)*(Pn_order+1)*(Pn_order+1)); + int this_node_lid = 0; + for (int k = 0; k <= Pn_order; k++) { for (int j = 0; j <= Pn_order; j++) { for (int i = 0; i <= Pn_order; i++) { - size_t node_lid = PointIndexFromIJK(i, j, k, order); - fprintf(out[0], "%lu ", nodes_in_elem_host(elem_gid, node_lid)); + // convert this_point index to the FE index convention + size_t vtk_index = PointIndexFromIJK(i, j, k, order); + + // store the points in this elem according the the finite + // element numbering convention + convert_fierro_to_vtk(vtk_index) = this_node_lid; + // increment the point counting index + this_node_lid ++; + + } // end for icount + } // end for jcount + } // end for kcount + // Write connectivity: all node IDs for all elements, space-separated + for (size_t elem_gid = 0; elem_gid < num_elems; elem_gid++) { + + int this_node_lid = 0; + for (int k = 0; k <= Pn_order; k++) { + for (int j = 0; j <= Pn_order; j++) { + for (int i = 0; i <= Pn_order; i++) { + // size_t node_lid = PointIndexFromIJK(i, j, k, order); + int node_lid = convert_fierro_to_vtk(this_node_lid); + fprintf(out[0], "%lu ", nodes_in_elem_host(elem_gid, node_lid)); + this_node_lid ++; + } } } } // end for @@ -4501,29 +4871,1761 @@ class MeshWriter } // end write vtu - ///////////////////////////////////////////////////////////////////////////// /// - /// \fn write_pvd - /// - /// \brief Writes a pvd ASCII output file for the element and nodal fields + /// \fn writes mesh with the format given in the input.yaml file /// - /// \param Vector of all graphics output times - /// \param element average field names - /// \param current time value - /// \param graphics index + /// \param Simulation mesh + /// \param Element related state + /// \param Node related state + /// \param Corner related state + /// \param Simulation input parameters /// ///////////////////////////////////////////////////////////////////////////// - void write_pvd(CArray& graphics_times, - double time_value, - int graphics_id, - const size_t solver_id, - int mpi_rank) - { - if (mpi_rank != 0) { - return; - } - + void write_mesh_Pn(swage::Mesh& mesh, + State_t& State, + SimulationParameters_t& SimulationParamaters, + double dt, + double time_value, + CArray graphics_times, + std::vector node_states, + std::vector gauss_pt_states, + std::vector material_pt_states, + const size_t solver_id, + elements::fe_ref_elem_t& ref_elem) + { + + + // node_state is an enum for possible fields (e.g., coords, velocity, etc.), see state.h + // gauss_pt_state is an enum for possible fields (e.g., vol, divergence, etc.) + // material_pt_state is an enum for possible fields (e.g., den, pres, etc.) + + + // ******************* + // Update host + // ******************* + + const size_t num_mats = State.MaterialPoints.num_material_points.size(); + + // material point values + + // Update host data for mat_pt state + for (auto field : material_pt_states){ + switch(field){ + // scalar vars to write out + case material_pt_state::density: + State.MaterialPoints.den.update_host(); + break; + case material_pt_state::pressure: + State.MaterialPoints.pres.update_host(); + break; + case material_pt_state::specific_internal_energy: + State.MaterialPoints.sie.update_host(); + break; + case material_pt_state::sound_speed: + State.MaterialPoints.sspd.update_host(); + break; + case material_pt_state::mass: + State.MaterialPoints.mass.update_host(); + break; + case material_pt_state::volume_fraction: + State.MaterialPoints.mat_volfrac.update_host(); + State.MaterialPoints.geo_volfrac.update_host(); + break; + case material_pt_state::eroded_flag: + State.MaterialPoints.eroded.update_host(); + break; + // tensor vars to write out + case material_pt_state::stress: + State.MaterialPoints.stress.update_host(); + break; + + case material_pt_state::strain: + State.MaterialPoints.strain.update_host(); + break; + + // additional vars for thermal-mechanical solver + case material_pt_state::thermal_conductivity: + State.MaterialPoints.conductivity.update_host(); + break; + + case material_pt_state::specific_heat: + State.MaterialPoints.specific_heat.update_host(); + break; + + case material_pt_state::heat_flux: + State.MaterialPoints.q_flux.update_host(); + break; + + // add other variables here + + // not used + case material_pt_state::elastic_modulii: + break; + case material_pt_state::shear_modulii: + break; + case material_pt_state::poisson_ratios: + break; + + default: + std::cout<<"Desired material point state not understood in outputs"< elem_scalar_var_names(num_elem_scalar_vars); + std::vector elem_tensor_var_names(num_elem_tensor_vars); + + // Scalar, vector, and tensor values associated with a material in part elems + std::vector mat_elem_scalar_var_names(num_mat_pt_scalar_vars); + std::vector mat_elem_tensor_var_names(num_mat_pt_tensor_vars); + + + // the ids to access a variable in the mat_scalar_var_name or tensor list + int mat_den_id = -1; + int mat_pres_id = -1; + int mat_sie_id = -1; + int mat_sspd_id = -1; + int mat_mass_id = -1; + int mat_mat_volfrac_id = -1; + int mat_geo_volfrac_id = -1; // geometric volume fraction of part + int mat_eroded_id = -1; + int mat_stress_id = -1; + int mat_strain_id = -1; + int mat_heat_flux_id = -1; + + int mat_conductivity_id = -1; + int mat_specific_heat_id = -1; + + // the index for the scalar, vector, and tensor fields + size_t var = 0; + size_t vector_var = 0; + size_t tensor_var = 0; + + // material point state to output + for (auto field : SimulationParamaters.OutputOptions.output_mat_pt_state){ + switch(field){ + // scalar vars + case material_pt_state::density: + mat_elem_scalar_var_names[var] = "mat_den"; + mat_den_id = var; + var++; + break; + case material_pt_state::pressure: + mat_elem_scalar_var_names[var] = "mat_pres"; + mat_pres_id = var; + var++; + break; + case material_pt_state::specific_internal_energy: + mat_elem_scalar_var_names[var] = "mat_sie"; + mat_sie_id = var; + var++; + break; + case material_pt_state::sound_speed: + mat_elem_scalar_var_names[var] = "mat_sspd"; + mat_sspd_id = var; + var++; + break; + case material_pt_state::mass: + mat_elem_scalar_var_names[var] = "mat_mass"; + mat_mass_id = var; + var++; + break; + case material_pt_state::volume_fraction: + mat_elem_scalar_var_names[var] = "mat_volfrac"; + mat_mat_volfrac_id = var; + var++; + + mat_elem_scalar_var_names[var] = "mat_geo_volfrac"; + mat_geo_volfrac_id = var; + var++; + break; + case material_pt_state::eroded_flag: + mat_elem_scalar_var_names[var] = "mat_eroded"; + mat_eroded_id = var; + var++; + break; + // tensor vars + case material_pt_state::stress: + mat_elem_tensor_var_names[tensor_var] = "mat_stress"; + mat_stress_id = tensor_var; + tensor_var++; + break; + + case material_pt_state::strain: + mat_elem_tensor_var_names[tensor_var] = "mat_strain"; + mat_strain_id = tensor_var; + tensor_var++; + break; + + + // additional vars for thermal-mechanical solver + case material_pt_state::thermal_conductivity: + mat_elem_scalar_var_names[var] = "mat_thermal_K"; + mat_conductivity_id = var; + var++; + break; + + case material_pt_state::specific_heat: + mat_elem_scalar_var_names[var] = "mat_Cp"; + mat_specific_heat_id = var; + var++; + break; + + case material_pt_state::heat_flux: + mat_elem_scalar_var_names[var] = "mat_heat_flux"; + mat_heat_flux_id = var; + var++; + break; + + + // add other variables here + + // not used + case material_pt_state::elastic_modulii: + break; + case material_pt_state::shear_modulii: + break; + case material_pt_state::poisson_ratios: + break; + + } // end switch + } // end for over mat_pt_states + + + // element average fields to output + + // the ids to access a variable in the elem_scalar_var_name or tensor list + int den_id = -1; + int pres_id = -1; + int sie_id = -1; + int sspd_id = -1; + int mass_id = -1; + int stress_id = -1; + int strain_id = -1; + + int conductivity_id = -1; + int specific_heat_id = -1; + + // reset the counters + var = 0; + vector_var = 0; + tensor_var = 0; + + // element state to output + for (auto field : SimulationParamaters.OutputOptions.output_elem_state){ + switch(field){ + // scalar vars + case material_pt_state::density: + elem_scalar_var_names[var] = "den"; + den_id = var; + var++; + break; + case material_pt_state::pressure: + elem_scalar_var_names[var] = "pres"; + pres_id = var; + var++; + break; + case material_pt_state::specific_internal_energy: + elem_scalar_var_names[var] = "sie"; + sie_id = var; + var++; + break; + case material_pt_state::sound_speed: + elem_scalar_var_names[var] = "sspd"; + sspd_id = var; + var++; + break; + case material_pt_state::mass: + elem_scalar_var_names[var] = "mass"; + mass_id = var; + var++; + break; + // tensor vars + case material_pt_state::stress: + elem_tensor_var_names[tensor_var] = "stress"; + stress_id = tensor_var; + tensor_var++; + break; + + case material_pt_state::strain: + elem_tensor_var_names[tensor_var] = "strain"; + strain_id = tensor_var; + tensor_var++; + break; + + // heat transfer variables + case material_pt_state::thermal_conductivity: + elem_scalar_var_names[var] = "thermal_K"; + conductivity_id = var; + var++; + break; + + case material_pt_state::specific_heat: + elem_scalar_var_names[var] = "Cp"; + specific_heat_id = var; + var++; + break; + + // add other variables here + + // not used + case material_pt_state::volume_fraction: + break; + case material_pt_state::eroded_flag: + break; + case material_pt_state::elastic_modulii: + break; + case material_pt_state::shear_modulii: + break; + case material_pt_state::poisson_ratios: + break; + case material_pt_state::heat_flux: + break; + } // end switch + } // end for over mat_pt_states + + // append Gauss point vars to the element arrays + int vol_id = -1; + int div_id = -1; + int level_set_id = -1; + int vel_grad_id = -1; + int shock_detector_id = -1; + + + for (auto field : SimulationParamaters.OutputOptions.output_gauss_pt_state){ + switch(field){ + // scalars + case gauss_pt_state::volume: + elem_scalar_var_names[var] = "vol"; + vol_id = var; + var++; + break; + case gauss_pt_state::divergence_velocity: + elem_scalar_var_names[var] = "div"; + div_id = var; + var++; + break; + + case gauss_pt_state::shock_detector: + elem_scalar_var_names[var] = "shock_detector"; + shock_detector_id = var; + var++; + break; + + case gauss_pt_state::level_set: + elem_scalar_var_names[var] = "level_set"; + level_set_id = var; + var++; + break; + + // tensors + case gauss_pt_state::gradient_velocity: + elem_tensor_var_names[tensor_var] = "vel_grad"; + vel_grad_id = tensor_var; + tensor_var++; + break; + } // end switch + } // end loop over gauss_pt_states + + + // ******************* + // nodal values + // ******************* + + size_t num_node_scalar_vars = 0; + size_t num_node_vector_vars = 0; + + for (auto field : SimulationParamaters.OutputOptions.output_node_state){ + switch(field){ + // --- scalars + case node_state::mass: + num_node_scalar_vars ++; + break; + case node_state::temp: + num_node_scalar_vars ++; + break; + // -- vectors + case node_state::coords: + num_node_vector_vars ++; + break; + case node_state::velocity: + num_node_vector_vars ++; // for velocity + num_node_vector_vars ++; // for acceleration + break; + case node_state::gradient_level_set: + num_node_vector_vars ++; + break; + case node_state::displacement: + num_node_vector_vars ++; + break; + case node_state::coords_t0: + // not necessary to output as coords at t=0 will be coords_t0 + break; + case node_state::force: + break; + + // heat transer vars + case node_state::heat_transfer: + break; + } // end switch + } // end for over + Kokkos::fence(); + + + // Scalar and vector values associated with a node + std::vector node_scalar_var_names(num_node_scalar_vars); + std::vector node_vector_var_names(num_node_vector_vars); + + int node_mass_id = -1; + int node_vel_id = -1; + int node_accel_id = -1; + int node_coord_id = -1; + int node_temp_id = -1; + int node_grad_level_set_id = -1; + int node_disp_id = -1; + + // reset counters for node fields + var = 0; + vector_var = 0; + tensor_var = 0; + + for (auto field : SimulationParamaters.OutputOptions.output_node_state){ + switch(field){ + // scalars + case node_state::mass: + node_scalar_var_names[var] = "node_mass"; + node_mass_id = var; + var++; + break; + case node_state::temp: + node_scalar_var_names[var] = "node_temp"; + node_temp_id = var; + var++; + break; + + // vector fields + + case node_state::coords: + node_vector_var_names[vector_var] = "node_coords"; + node_coord_id = vector_var; + vector_var++; + break; + + case node_state::velocity: + node_vector_var_names[vector_var] = "node_vel"; + node_vel_id = vector_var; + vector_var++; + + node_vector_var_names[vector_var] = "node_accel"; + node_accel_id = vector_var; + vector_var++; + break; + + case node_state::gradient_level_set: + node_vector_var_names[vector_var] = "node_grad_lvlset"; + node_grad_level_set_id = vector_var; + vector_var++; + break; + + case node_state::displacement: + node_vector_var_names[vector_var] = "node_disp"; + node_disp_id = vector_var; + vector_var ++; + break; + + // -- not used vars + case node_state::coords_t0: + break; + + case node_state::force: + break; + + // heat transer vars + case node_state::heat_transfer: + break; + + // tensors + + } // end switch + } // end for over + + + // ************************************** + // build and save element average fields + // ************************************** + + // short hand + const size_t num_nodes = mesh.num_nodes; + const size_t num_elems = mesh.num_elems; + const size_t num_dims = mesh.num_dims; + const size_t num_nodes_in_elem = mesh.num_nodes_in_elem; + const int Pn_order = mesh.Pn; + + // save the elem state to an array for exporting to graphics files + DCArrayKokkos elem_scalar_fields(num_elem_scalar_vars, num_elems, "elem_scalars"); + DCArrayKokkos elem_tensor_fields(num_elem_tensor_vars, num_elems, 3, 3, "elem_tensors"); + elem_scalar_fields.set_values(0.0); + elem_tensor_fields.set_values(0.0); + + + // ----------------------------------------------------------------------- + // save the output fields to a single element average array for all state + // ----------------------------------------------------------------------- + for (int mat_id = 0; mat_id < num_mats; mat_id++) { + + // material point and guass point state are concatenated together + concatenate_elem_fields(State.MaterialPoints, + State.GaussPoints, + elem_scalar_fields, + elem_tensor_fields, + State.MaterialToMeshMaps.elem_in_mat_elem, + SimulationParamaters.OutputOptions.output_elem_state, + SimulationParamaters.OutputOptions.output_gauss_pt_state, + State.MaterialToMeshMaps.num_mat_elems.host(mat_id), + mat_id, + num_elems, + den_id, + pres_id, + sie_id, + sspd_id, + mass_id, + stress_id, + strain_id, + vol_id, + div_id, + shock_detector_id, + level_set_id, + vel_grad_id, + conductivity_id, + specific_heat_id); + } // end for mats + + // make specific fields for the element average + if (sie_id>=0){ + FOR_ALL(elem_gid, 0, num_elems, { + // get sie by dividing by the mass + elem_scalar_fields(sie_id, elem_gid) /= (elem_scalar_fields(mass_id, elem_gid)+1.e-20); + }); + } // end if + + Kokkos::fence(); + elem_scalar_fields.update_host(); + elem_tensor_fields.update_host(); + + + // ************************ + // Build the nodal fields + // ************************ + + // save the nodal fields to an array for exporting to graphics files + DCArrayKokkos node_scalar_fields(num_node_scalar_vars, num_nodes, "node_scalars"); + DCArrayKokkos node_vector_fields(num_node_vector_vars, num_nodes, 3, "node_tenors"); + + concatenate_nodal_fields(State.node, + node_scalar_fields, + node_vector_fields, + SimulationParamaters.OutputOptions.output_node_state, + dt, + num_nodes, + num_dims, + node_mass_id, + node_vel_id, + node_accel_id, + node_coord_id, + node_grad_level_set_id, + node_temp_id, + node_disp_id); + + + Kokkos::fence(); + node_scalar_fields.update_host(); + node_vector_fields.update_host(); + + + // ******************************** + // Write the nodal and elem fields + // ******************************** + + if (SimulationParamaters.OutputOptions.format == output_options::viz || + SimulationParamaters.OutputOptions.format == output_options::viz_and_state) { + + int mpi_rank = 0; + int mpi_size = 1; + mesh_io_mpi_detail::query_world_rank_size(mpi_rank, mpi_size); + + MPI_Barrier(MPI_COMM_WORLD); + + if (mpi_rank == 0 && solver_id == 0 && graphics_id == 0) { + struct stat st_rm; + if (stat("vtk", &st_rm) == 0) { + (void)system("rm -f vtk/Fierro*"); + } + if (stat("vtk/data", &st_rm) == 0) { + (void)system("rm -f vtk/data/Fierro*"); + } + } + + MPI_Barrier(MPI_COMM_WORLD); + + struct stat st; + if (stat("vtk", &st) != 0) { + if (system("mkdir vtk") == 1) { + std::cout << "Unable to make vtk directory" << std::endl; + } + } + if (stat("vtk/data", &st) != 0) { + if (system("mkdir vtk/data") == 1) { + std::cout << "Unable to make vtk/data directory" << std::endl; + } + } + + // Per-rank VTU: write owned elements only (rows 0..num_owned_elems-1 of nodes_in_elem), but + // Points must list all local nodes (0..num_nodes-1). Connectivity is full local indexing; + // boundary elements reference ghost nodes — truncating coords to num_owned_nodes (the prior + // bug) misaligned indices and produced inverted/wrong cells in ParaView. + const size_t n_owned_elems = mesh.num_owned_elems; // num_elems; + + + const std::string elem_fields_name = "fields"; + + ViewCArray node_coords_host(&State.node.coords.host(0, 0), num_nodes, num_dims); + ViewCArray nodes_in_elem_host(&mesh.nodes_in_elem.host(0, 0), n_owned_elems, num_nodes_in_elem); + + // VTK diagnostics (optional CellData / PointData in write_vtu): host-side scratch arrays, + // then const pointers so write_vtu can emit one value per owned cell or per owned node. + // diag_mpi_rank_elem — CellData "mpi_rank" (which MPI rank owns each exported element). + // diag_global_elem — CellData "global_elem_id" from mesh.local_to_global_elem_mapping (or local e). + // diag_global_node — PointData "global_node_id" from mesh.local_to_global_node_mapping (or local n). + // p_* — nullptr until filled; passed to write_vtu (nullptr disables that array). + std::vector diag_mpi_rank_elem; + std::vector diag_global_elem; + std::vector diag_global_node; + const double* p_rank_elem = nullptr; + const double* p_glob_elem = nullptr; + const double* p_glob_node = nullptr; + + if (n_owned_elems > 0 && num_nodes > 0) { + diag_mpi_rank_elem.assign(n_owned_elems, static_cast(mpi_rank)); + p_rank_elem = diag_mpi_rank_elem.data(); + + diag_global_elem.resize(n_owned_elems); + if (mpi_size > 1 || mesh.num_elems > mesh.num_owned_elems) { + mesh.local_to_global_elem_mapping.update_host(); + for (size_t e = 0; e < n_owned_elems; e++) { + diag_global_elem[e] = + static_cast(mesh.local_to_global_elem_mapping.host(e)); + } + } + else { + for (size_t e = 0; e < n_owned_elems; e++) { + diag_global_elem[e] = static_cast(e); + } + } + p_glob_elem = diag_global_elem.data(); + + diag_global_node.resize(num_nodes); + if (mpi_size > 1 || mesh.num_nodes > mesh.num_owned_nodes) { + mesh.local_to_global_node_mapping.update_host(); + for (size_t n = 0; n < num_nodes; n++) { + diag_global_node[n] = + static_cast(mesh.local_to_global_node_mapping.host(n)); + } + } + else { + for (size_t n = 0; n < num_nodes; n++) { + diag_global_node[n] = static_cast(n); + } + } + p_glob_node = diag_global_node.data(); + + write_vtu(node_coords_host, + nodes_in_elem_host, + elem_scalar_fields, + elem_tensor_fields, + node_scalar_fields, + node_vector_fields, + elem_scalar_var_names, + elem_tensor_var_names, + node_scalar_var_names, + node_vector_var_names, + elem_fields_name, + graphics_id, + num_nodes, + n_owned_elems, + num_nodes_in_elem, + Pn_order, + num_dims, + solver_id, + mpi_rank, + mpi_size, + p_rank_elem, + p_glob_elem, + p_glob_node); + } + + + // ******************************** + // Build and write the mat fields + // ******************************** + + // note: the file path and folder was created in the elem and node outputs + size_t num_mat_files_written = 0; + if(num_mat_pt_scalar_vars > 0 || num_mat_pt_tensor_vars >0){ + + for (int mat_id = 0; mat_id < num_mats; mat_id++) { + + const size_t num_mat_elems = State.MaterialToMeshMaps.num_mat_elems.host(mat_id); + + // only save material data if the mat lives on the mesh, ie. has state allocated + if (num_mat_elems>0){ + + // set the nodal vars to zero size, we don't write these fields again + node_scalar_var_names.clear(); + node_vector_var_names.clear(); + + // the arrays storing all the material field data + DCArrayKokkos mat_elem_scalar_fields(num_mat_pt_scalar_vars, num_mat_elems, "mat_pt_scalars"); + DCArrayKokkos mat_elem_tensor_fields(num_mat_pt_tensor_vars, num_mat_elems, 3, 3, "mat_pt_tensors"); + + + // concatenate material fields into a single array + concatenate_mat_fields(State.MaterialPoints, + mat_elem_scalar_fields, + mat_elem_tensor_fields, + State.MaterialToMeshMaps.elem_in_mat_elem, + SimulationParamaters.OutputOptions.output_mat_pt_state, + num_mat_elems, + mat_id, + mat_den_id, + mat_pres_id, + mat_sie_id, + mat_sspd_id, + mat_mass_id, + mat_mat_volfrac_id, + mat_geo_volfrac_id, + mat_eroded_id, + mat_stress_id, + mat_strain_id, + mat_conductivity_id, + mat_specific_heat_id, + mat_heat_flux_id); + Kokkos::fence(); + mat_elem_scalar_fields.update_host(); + mat_elem_tensor_fields.update_host(); + + + std::string str_mat_val = std::to_string(mat_id); + std::string mat_fields_name = "mat"; + mat_fields_name += str_mat_val; // add the mat number + + // save the nodes belonging to this part (i.e., the material) + DCArrayKokkos mat_node_coords(num_nodes,num_dims, "mat_node_coords"); + DCArrayKokkos mat_nodes_in_mat_elem(num_mat_elems, num_nodes_in_elem, "mat_nodes_in_mat_elem"); + + // the number of actual nodes belonging to the part (i.e., the material) + size_t num_mat_nodes = 0; + + // build a unique mesh (element and nodes) for the material (i.e., the part) + build_material_elem_node_lists(mesh, + State.node.coords, + mat_node_coords, + mat_nodes_in_mat_elem, + State.MaterialToMeshMaps.elem_in_mat_elem, + mat_id, + num_mat_nodes, + num_mat_elems, + num_nodes_in_elem, + num_dims); + + ViewCArray mat_node_coords_host(&mat_node_coords.host(0,0), num_mat_nodes, num_dims); + ViewCArray mat_nodes_in_elem_host(&mat_nodes_in_mat_elem.host(0,0), num_mat_elems, num_nodes_in_elem); + + // write out a vtu file this + write_vtu(mat_node_coords_host, + mat_nodes_in_elem_host, + mat_elem_scalar_fields, + mat_elem_tensor_fields, + node_scalar_fields, + node_vector_fields, + mat_elem_scalar_var_names, + mat_elem_tensor_var_names, + node_scalar_var_names, + node_vector_var_names, + mat_fields_name, + graphics_id, + num_mat_nodes, + num_mat_elems, + num_nodes_in_elem, + Pn_order, + num_dims, + solver_id, + mpi_rank, + mpi_size, + nullptr, + nullptr, + nullptr); + + // --- Dataset B: Gauss-point cloud for ParaView interpolation --- + if (num_mat_pt_scalar_vars > 0 || num_mat_pt_tensor_vars > 0) { + + std::string pn_name = "mat" + str_mat_val + "_Pn"; + + write_vtu_Pn(mesh, + State, + ref_elem, + mat_elem_scalar_var_names, + mat_elem_tensor_var_names, + pn_name, + graphics_id, + mat_id, + num_mat_elems, + num_nodes_in_elem, + num_dims, + solver_id, + mat_den_id, + mat_pres_id, + mat_sie_id, + mat_sspd_id, + mat_mass_id, + mat_mat_volfrac_id, + mat_geo_volfrac_id, + mat_eroded_id, + mat_stress_id, + mat_strain_id, + mat_conductivity_id, + mat_specific_heat_id); + + } + + + num_mat_files_written++; + + } // end for mat_id + + } // end if material is on the mesh + + } // end if mat variables are to be written + + + // ************************************************* + // write Paraview files to open the graphics files + // ************************************************* + + // save the graphics time + graphics_times(graphics_id) = time_value; + + // check to see if an mesh state was written + bool write_mesh_state = false; + if( num_elem_scalar_vars > 0 || + num_elem_tensor_vars > 0 || + num_node_scalar_vars > 0 || + num_node_vector_vars > 0) + { + write_mesh_state = true; + } + + // check to see if a mat state was written + bool write_mat_pt_state = false; + if( num_mat_pt_scalar_vars > 0 || + num_mat_pt_tensor_vars > 0) + { + write_mat_pt_state = true; + } + + // call the vtm file writer + std::string mat_fields_name = "mat"; + write_vtm_Pn(graphics_times, + elem_fields_name, + mat_fields_name, + time_value, + graphics_id, + num_mats, + write_mesh_state, + write_mat_pt_state, + solver_id); + + // call the pvd file writer + write_pvd(graphics_times, + time_value, + graphics_id, + solver_id, + mpi_rank); + + + // increment graphics id counter + graphics_id++; // this is private variable in the class + + } // end if viz paraview output is to be written + + + // STATE + if (SimulationParamaters.OutputOptions.format == output_options::state || + SimulationParamaters.OutputOptions.format == output_options::viz_and_state) { + + /* write_material_point_state(mesh, + State, + SimulationParamaters, + time_value, + graphics_times, + node_states, + gauss_pt_states, + material_pt_states); */ + write_text_state_Pn(mesh, + State, + ref_elem, + mat_elem_scalar_var_names, + mat_elem_tensor_var_names, + node_scalar_var_names, + node_vector_var_names, + mesh.num_nodes_in_elem, + mesh.num_dims, + time_value, + mat_den_id, + mat_pres_id, + mat_sie_id, + mat_sspd_id, + mat_mass_id, + mat_mat_volfrac_id, + mat_geo_volfrac_id, + mat_eroded_id, + mat_stress_id, + mat_strain_id, + mat_conductivity_id, + mat_specific_heat_id, + node_mass_id, + node_vel_id, + node_coord_id, + node_temp_id, + node_grad_level_set_id, + node_disp_id); + + } // end if state is to be written + + + // will drop ensight outputs in the near future + if (SimulationParamaters.OutputOptions.format == output_options::ensight){ + write_ensight(mesh, + State, + SimulationParamaters, + dt, + time_value, + graphics_times, + node_states, + gauss_pt_states, + material_pt_states); + } + + return; + + } // end write_mesh_Pn + + ///////////////////////////////////////////////////////////////////////////// + /// + /// \fn write_vtu_Pn + /// + /// \brief Writes a vtu ASCII output file for gauss points + /// + /// \param Simulation mesh + /// \param State data + /// \param Simulation parameters + /// \param current time value + /// \param Vector of all graphics output times + /// + ///////////////////////////////////////////////////////////////////////////// + + void write_vtu_Pn( + const swage::Mesh& mesh, + const State_t& State, + elements::fe_ref_elem_t& ref_elem, + const std::vector& mat_scalar_var_names, // length == num_mat_pt_scalar_vars + const std::vector& mat_tensor_var_names, // length == num_mat_pt_tensor_vars + const std::string& partname, // e.g. "mat0_Pn" + const int graphics_id, + const int mat_id, + const size_t num_mat_elems, + const size_t num_nodes_in_elem, + const size_t num_dims, + const size_t solver_id, + // field slot IDs (-1 means "not requested") + const int mat_den_id, + const int mat_pres_id, + const int mat_sie_id, + const int mat_sspd_id, + const int mat_mass_id, + const int mat_mat_volfrac_id, + const int mat_geo_volfrac_id, + const int mat_eroded_id, + const int mat_stress_id, + const int mat_strain_id, + const int mat_conductivity_id, + const int mat_specific_heat_id + ) + { + // ----------------------------------------------------------------------- + // Derived sizes + // ----------------------------------------------------------------------- + + // Number of Gauss points per element: (2*Pn)^num_dims + // Read directly from the reference element so we stay consistent with + // whatever quadrature order was set up for this run. + const size_t num_gp_per_elem = ref_elem.gauss_point_grad_basis.dims(0); + const size_t num_total_gp = num_mat_elems * num_gp_per_elem; + + const size_t num_scalar_vars = mat_scalar_var_names.size(); + const size_t num_tensor_vars = mat_tensor_var_names.size(); + + // ----------------------------------------------------------------------- + // Open file + // ----------------------------------------------------------------------- + + FILE* fp; + char filename[128]; + int max_len = sizeof filename; + + // File lives next to the standard per-material .vtu files. + // The "_Pn" suffix distinguishes this point-cloud dataset from the + // element-averaged dataset written by write_vtu. + snprintf(filename, max_len, + "vtk/data/Fierro.solver%zu.%s.%05d.vtu", + solver_id, partname.c_str(), graphics_id); + + fp = fopen(filename, "w"); + if (!fp) { + std::cerr << "write_vtu_Pn: could not open " << filename << std::endl; + return; + } + + // ----------------------------------------------------------------------- + // VTK file header + // ----------------------------------------------------------------------- + + fprintf(fp, "\n"); + fprintf(fp, "\n"); + fprintf(fp, " \n"); + fprintf(fp, " \n", + num_total_gp, num_total_gp); + + // ----------------------------------------------------------------------- + // Points – physical coordinates of every Gauss point + // + // Mapping: x_phys[dim] = Σ_a N_a(ξ_gp) * x_node_a[dim] + // ----------------------------------------------------------------------- + + fprintf(fp, "\n"); + fprintf(fp, " \n"); + fprintf(fp, " \n"); + fprintf(fp, " \n"); + + DCArrayKokkos x_phys(num_mat_elems, num_gp_per_elem, 3); + x_phys.set_values(0); + FOR_ALL(elem, 0, num_mat_elems, { + const size_t elem_id = + State.MaterialToMeshMaps.elem_in_mat_elem(mat_id, elem); + + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + + for (size_t node_lid = 0; node_lid < num_nodes_in_elem; node_lid++) { + const size_t node_gid = mesh.nodes_in_elem(elem_id, node_lid); + const double N = ref_elem.gauss_point_basis(gp, node_lid); + + for (size_t dim = 0; dim < num_dims; dim++) { + x_phys(elem, gp, dim) += N * State.node.coords(node_gid, dim); + } + } + + } // end gp + }); + x_phys.update_host(); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + + // VTK always expects 3 components; pad Z to 0 for 2-D runs. + fprintf(fp, " %.15e %.15e %.15e\n", + x_phys.host(elem, gp, 0), + x_phys.host(elem, gp, 1), + (num_dims == 3) ? x_phys.host(elem, gp, 2) : 0.0); + + } // end gp + } // end elem + + fprintf(fp, " \n"); + fprintf(fp, " \n"); + + // ----------------------------------------------------------------------- + // Cells – one VTK_VERTEX (type = 1) per Gauss point + // + // connectivity[i] = i (each vertex cell references exactly one point) + // offsets[i] = i+1 + // types[i] = 1 (VTK_VERTEX) + // ----------------------------------------------------------------------- + + fprintf(fp, "\n"); + fprintf(fp, " \n"); + fprintf(fp, " \n"); + + fprintf(fp, " \n"); + for (size_t pt = 0; pt < num_total_gp; pt++) { + fprintf(fp, " %zu\n", pt); + } + fprintf(fp, " \n"); + + fprintf(fp, " \n"); + for (size_t pt = 0; pt < num_total_gp; pt++) { + fprintf(fp, " %zu\n", pt + 1); + } + fprintf(fp, " \n"); + + fprintf(fp, " \n"); + for (size_t pt = 0; pt < num_total_gp; pt++) { + fprintf(fp, " 1\n"); // VTK_VERTEX + } + fprintf(fp, " \n"); + + fprintf(fp, " \n"); + + // ----------------------------------------------------------------------- + // PointData – material-point fields at every Gauss point + // + // The flat material-point index is obtained the same way the device + // kernels do it: pt_id = State.points_in_mat_elem.host(elem, gp) + // ----------------------------------------------------------------------- + + if (num_scalar_vars > 0 || num_tensor_vars > 0) { + + fprintf(fp, "\n"); + fprintf(fp, " \n"); + fprintf(fp, " \n"); + + // ------------------------------------------------------------------- + // Scalar fields + // ------------------------------------------------------------------- + + // Helper lambda-style macro to avoid repeating the loop boilerplate. + // Written as an ordinary block; C++11 lambdas would also work if the + // codebase already uses them. + + // density + if (mat_den_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_den_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.den.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // pressure + if (mat_pres_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_pres_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.pres.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // specific internal energy + if (mat_sie_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_sie_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.sie.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // sound speed + if (mat_sspd_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_sspd_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.sspd.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // mass + if (mat_mass_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_mass_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.mass.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // material volume fraction + if (mat_mat_volfrac_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_mat_volfrac_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.mat_volfrac.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // geometric volume fraction + if (mat_geo_volfrac_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_geo_volfrac_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.geo_volfrac.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // eroded flag (cast to double so VTK Float64 array stays consistent) + if (mat_eroded_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_eroded_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + static_cast( + State.MaterialPoints.eroded.host(mat_id, pt_id))); + } + } + fprintf(fp, " \n"); + } + + // thermal conductivity + if (mat_conductivity_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_conductivity_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.conductivity.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // specific heat + if (mat_specific_heat_id >= 0) { + fprintf(fp, " \n", + mat_scalar_var_names[mat_specific_heat_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " %.15e\n", + State.MaterialPoints.specific_heat.host(mat_id, pt_id)); + } + } + fprintf(fp, " \n"); + } + + // ------------------------------------------------------------------- + // Tensor fields (stored as 9-component Float64 arrays, row-major) + // + // ParaView convention: Txx Txy Txz Tyx Tyy Tyz Tzx Tzy Tzz + // which matches CArray row-major storage (i outer, j inner). + // ------------------------------------------------------------------- + + // stress + if (mat_stress_id >= 0) { + fprintf(fp, " \n", + mat_tensor_var_names[mat_stress_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " "); + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < 3; j++) { + fprintf(fp, " %.15e", + State.MaterialPoints.stress.host( + mat_id, pt_id, i, j)); + } + } + fprintf(fp, "\n"); + } + } + fprintf(fp, " \n"); + } + + // strain + if (mat_strain_id >= 0) { + fprintf(fp, " \n", + mat_tensor_var_names[mat_strain_id].c_str()); + for (size_t elem = 0; elem < num_mat_elems; elem++) { + for (size_t gp = 0; gp < num_gp_per_elem; gp++) { + const size_t pt_id = + State.points_in_mat_elem.host(elem, gp); + fprintf(fp, " "); + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < 3; j++) { + fprintf(fp, " %.15e", + State.MaterialPoints.strain.host( + mat_id, pt_id, i, j)); + } + } + fprintf(fp, "\n"); + } + } + fprintf(fp, " \n"); + } + + fprintf(fp, " \n"); + + } // end if any fields to write + + // ----------------------------------------------------------------------- + // Footer + // ----------------------------------------------------------------------- + + fprintf(fp, " \n"); + fprintf(fp, " \n"); + fprintf(fp, "\n"); + + fclose(fp); + + } // end write_vtu_Pn + + void write_vtm_Pn( + CArray& graphics_times, + const std::string& elem_part_name, + const std::string& mat_part_name, + double time_value, + int graphics_id, + int num_mats, + bool write_mesh_state, + bool write_mat_pt_state, + const size_t solver_id) + { + // The _Pn files are named mat_part_name + mat_id + "_Pn" + // e.g. "mat0_Pn" → Fierro.solver0.mat0_Pn.00001.vtu + // This matches the partname constructed in write_mesh: + // std::string pn_name = "mat" + str_mat_val + "_Pn"; + + for (int file_id = 0; file_id <= graphics_id; file_id++) { + + FILE* out[1]; + char filename[100]; + int max_len = sizeof filename; + int str_output_len; + + str_output_len = snprintf(filename, max_len, + "vtk/data/Fierro.solver%zu.%05d.vtm", solver_id, file_id); + if (str_output_len >= max_len) { + fputs("Filename length exceeded; string truncated", stderr); + } + + out[0] = fopen(filename, "w"); + + fprintf(out[0], "\n"); + fprintf(out[0], "\n"); + fprintf(out[0], " \n"); + + size_t block_id = 0; + + // ------------------------------------------------------------------- + // Block 0: element-averaged mesh (node + element state) + // ------------------------------------------------------------------- + if (write_mesh_state) { + fprintf(out[0], " \n", + block_id); + block_id++; + + fprintf(out[0], " \n"); + fprintf(out[0], " \n", + file_id, solver_id, elem_part_name.c_str(), file_id, + graphics_times(file_id)); + fprintf(out[0], " \n"); + + fprintf(out[0], " \n"); + } + + // ------------------------------------------------------------------- + // Block 1: per-material data + // + // Each material gets its own named sub-block so ParaView shows a + // clean tree. Inside each sub-block there are up to two Pieces: + // • "Elements" — element-averaged material fields (matN.vtu) + // • "GaussPoints" — Gauss-point cloud (matN_Pn.vtu) + // ------------------------------------------------------------------- + if (write_mat_pt_state) { + + fprintf(out[0], " \n", + block_id); + block_id++; + + for (size_t mat_id = 0; mat_id < (size_t)num_mats; mat_id++) { + + // Open a sub-block for this material + fprintf(out[0], + " \n", + mat_id, mat_id); + + size_t piece_id = 0; + + // Piece 0: element-averaged material fields + if (write_mat_pt_state) { + fprintf(out[0], + " \n", + piece_id); + fprintf(out[0], + " \n", + file_id, solver_id, mat_part_name.c_str(), + mat_id, file_id, graphics_times(file_id)); + fprintf(out[0], " \n"); + piece_id++; + } + + // Piece 1: Gauss-point cloud + if (write_mat_pt_state) { + fprintf(out[0], + " \n", + piece_id); + // filename pattern: mat_part_name + mat_id + "_Pn" + // e.g. Fierro.solver0.mat0_Pn.00001.vtu + fprintf(out[0], + " \n", + file_id, solver_id, mat_part_name.c_str(), + mat_id, file_id, graphics_times(file_id)); + fprintf(out[0], " \n"); + piece_id++; + } + + fprintf(out[0], " \n"); // close MatN sub-block + + } // end for mat_id + + fprintf(out[0], " \n"); // close Mat block + + } // end if any material state + + fprintf(out[0], " \n"); + fprintf(out[0], ""); + + fclose(out[0]); + + } // end for file_id + + } // end write_vtm_Pn + + + ///////////////////////////////////////////////////////////////////////////// + /// + /// \fn write_pvd + /// + /// \brief Writes a pvd ASCII output file for the element and nodal fields + /// + /// \param Vector of all graphics output times + /// \param element average field names + /// \param current time value + /// \param graphics index + /// + ///////////////////////////////////////////////////////////////////////////// + void write_pvd(CArray& graphics_times, + double time_value, + int graphics_id, + const size_t solver_id, + int mpi_rank) + { + if (mpi_rank != 0) { + return; + } + FILE* out[20]; char filename[512]; int max_len = static_cast(sizeof filename); @@ -4885,6 +6987,7 @@ class MeshWriter State.MaterialPoints.den.update_host(); State.MaterialPoints.pres.update_host(); State.MaterialPoints.stress.update_host(); + State.MaterialPoints.strain.update_host(); State.MaterialPoints.sspd.update_host(); State.MaterialPoints.sie.update_host(); State.MaterialPoints.mass.update_host(); @@ -5184,6 +7287,823 @@ class MeshWriter return; } // end of state output + + void write_text_state( + const swage::Mesh& mesh, + const State_t& State, + const std::vector& mat_scalar_var_names, + const std::vector& mat_tensor_var_names, + const std::vector& node_scalar_var_names, + const std::vector& node_vector_var_names, + const size_t num_nodes_in_elem, + const size_t num_dims, + const double time_value, + // material-point field slot IDs (-1 means "not requested") + const int mat_den_id, + const int mat_pres_id, + const int mat_sie_id, + const int mat_sspd_id, + const int mat_mass_id, + const int mat_mat_volfrac_id, + const int mat_geo_volfrac_id, + const int mat_eroded_id, + const int mat_stress_id, + const int mat_strain_id, + const int mat_conductivity_id, + const int mat_specific_heat_id, + const int mat_heat_flux_id, + // node field slot IDs (-1 means "not requested") + const int node_mass_id, + const int node_vel_id, + const int node_coord_id, + const int node_temp_id, + const int node_grad_level_set_id, + const int node_disp_id + ) + { + // ----------------------------------------------------------------------- + // Ensure the output directory exists before opening any files. + // ----------------------------------------------------------------------- + struct stat st; + MPI_Barrier(MPI_COMM_WORLD); + if (stat("state", &st) != 0) { + system("mkdir state"); + } + MPI_Barrier(MPI_COMM_WORLD); + + // ----------------------------------------------------------------------- + // FILE 1 – Material-point state + // + // One file per material. One row per element (single integration + // point assumed); the element coordinate is the centroid, computed + // as the simple average of its node positions. + // + // Filename: state/mat_pt_state_t_