Skip to content

Commit

Permalink
Merge pull request #237 from lanl/strength
Browse files Browse the repository at this point in the history 
Strength
  • Loading branch information
@nathanielmorgan
nathanielmorgan authored Oct 25, 2024
2 parents 480b078 + 1358816 commit a8c4cd2
Show file tree
Hide file tree
Showing 30 changed files with 1,468 additions and 578 deletions.
27 changes: 15 additions & 12 deletions single-node-refactor/src/Solvers/SGH_solver_3D/include/sgh_solver_3D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ namespace SGH3D_State
static const std::vector<gauss_pt_state> required_gauss_pt_state =
{
gauss_pt_state::volume,
gauss_pt_state::divergence_velocity
gauss_pt_state::gradient_velocity
};

// Material point state to be initialized for the SGH solver
Expand All @@ -79,7 +79,8 @@ namespace SGH3D_State
material_pt_state::mass,
material_pt_state::volume_fraction,
material_pt_state::specific_internal_energy,
material_pt_state::eroded_flag
material_pt_state::eroded_flag,
material_pt_state::shear_modulii
};

// Material corner state to be initialized for the SGH solver
Expand Down Expand Up @@ -228,7 +229,7 @@ class SGH3D : public Solver
const Material_t& Materials,
const Mesh_t& mesh,
const DCArrayKokkos<double>& GaussPoints_vol,
const DCArrayKokkos<double>& GaussPoints_div,
const DCArrayKokkos<double>& GaussPoints_vel_grad,
const DCArrayKokkos<bool>& GaussPoints_eroded,
const DCArrayKokkos<double>& corner_force,
const DCArrayKokkos<double>& node_coords,
Expand Down Expand Up @@ -267,14 +268,12 @@ class SGH3D : public Solver
const DCArrayKokkos<double>& node_mass,
const DCArrayKokkos<double>& corner_force) const;

KOKKOS_FUNCTION
void get_velgrad(
ViewCArrayKokkos<double>& vel_grad,
const ViewCArrayKokkos<size_t>& elem_node_gids,
const DCArrayKokkos<double>& node_vel,
const ViewCArrayKokkos<double>& b_matrix,
const double GaussPoints_vol,
const size_t elem_gid) const;
DCArrayKokkos<double>& vel_grad,
const Mesh_t mesh,
const DCArrayKokkos<double>& node_coords,
const DCArrayKokkos<double>& node_vel,
const DCArrayKokkos<double>& elem_vol) const;

void get_divergence(
DCArrayKokkos<double>& GaussPoints_div,
Expand All @@ -295,6 +294,7 @@ class SGH3D : public Solver
const Mesh_t& mesh,
const DCArrayKokkos<double>& node_coords,
const DCArrayKokkos<double>& node_vel,
const DCArrayKokkos<double>& GaussPoints_vel_grad,
const DCArrayKokkos<double>& MaterialPoints_den,
const DCArrayKokkos<double>& MaterialPoints_pres,
const DCArrayKokkos<double>& MaterialPoints_stress,
Expand All @@ -304,7 +304,8 @@ class SGH3D : public Solver
const DCArrayKokkos<double>& MaterialPoints_mass,
const DCArrayKokkos<double>& MaterialPoints_eos_state_vars,
const DCArrayKokkos<double>& MaterialPoints_strength_state_vars,
const DCArrayKokkos<bool>& GaussPoints_eroded,
const DCArrayKokkos<bool>& MaterialPoints_eroded,
const DCArrayKokkos<double>& MaterialPoints_shear_modulii,
const DCArrayKokkos<size_t>& MaterialToMeshMaps_elem,
const double time_value,
const double dt,
Expand All @@ -319,13 +320,15 @@ class SGH3D : public Solver
const DCArrayKokkos<double>& GaussPoints_vol,
const DCArrayKokkos<double>& node_coords,
const DCArrayKokkos<double>& node_vel,
const DCArrayKokkos<double>& GaussPoints_vel_grad,
const DCArrayKokkos<double>& MaterialPoints_den,
const DCArrayKokkos<double>& MaterialPoints_sie,
const DCArrayKokkos<double>& MaterialPoints_pres,
const DCArrayKokkos<double>& MaterialPoints_stress,
const DCArrayKokkos<double>& MaterialPoints_sspd,
const DCArrayKokkos<double>& MaterialPoints_eos_state_vars,
const DCArrayKokkos<double>& MaterialPoints_strength_state_vars,
const DCArrayKokkos<double>& MaterialPoints_shear_modulii,
const DCArrayKokkos<size_t>& MaterialToMeshMaps_elem,
const size_t num_mat_elems,
const size_t mat_id,
Expand Down Expand Up @@ -389,7 +392,7 @@ class SGH3D : public Solver
const DCArrayKokkos<double>& MaterialPoints_sspd,
const double den,
const double sie,
const ViewCArrayKokkos<double>& vel_grad,
const DCArrayKokkos<double>& GaussPoints_vel_grad,
const ViewCArrayKokkos<size_t>& elem_node_gids,
const DCArrayKokkos<double>& node_coords,
const DCArrayKokkos<double>& node_vel,
Expand Down
28 changes: 12 additions & 16 deletions single-node-refactor/src/Solvers/SGH_solver_3D/src/force_sgh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
void SGH3D::get_force(const Material_t& Materials,
const Mesh_t& mesh,
const DCArrayKokkos<double>& GaussPoints_vol,
const DCArrayKokkos<double>& GaussPoints_div,
const DCArrayKokkos<double>& GaussPoints_vel_grad,
const DCArrayKokkos<bool>& MaterialPoints_eroded,
const DCArrayKokkos<double>& corner_force,
const DCArrayKokkos<double>& node_coords,
Expand Down Expand Up @@ -120,8 +120,7 @@ void SGH3D::get_force(const Material_t& Materials,
// Riemann velocity
double vel_star_array[3];

// velocity gradient
double vel_grad_array[9];


// --- Create views of arrays to aid the force calculation ---

Expand All @@ -131,7 +130,7 @@ void SGH3D::get_force(const Material_t& Materials,
ViewCArrayKokkos<double> sum(sum_array, 4);
ViewCArrayKokkos<double> muc(muc_array, num_nodes_in_elem);
ViewCArrayKokkos<double> vel_star(vel_star_array, num_dims);
ViewCArrayKokkos<double> vel_grad(vel_grad_array, num_dims, num_dims);


// element volume
double vol = GaussPoints_vol(elem_gid);
Expand All @@ -149,13 +148,6 @@ void SGH3D::get_force(const Material_t& Materials,
node_coords,
elem_node_gids);

// --- Calculate the velocity gradient ---
SGH3D::get_velgrad(vel_grad,
elem_node_gids,
node_vel,
area_normal,
vol,
elem_gid);

// the -1 is for the inward surface area normal,
for (size_t node_lid = 0; node_lid < num_nodes_in_elem; node_lid++) {
Expand All @@ -164,16 +156,18 @@ void SGH3D::get_force(const Material_t& Materials,
} // end for
} // end for

double div = GaussPoints_div(elem_gid);
double div = GaussPoints_vel_grad(elem_gid, 0, 0) +
GaussPoints_vel_grad(elem_gid, 1, 1) +
GaussPoints_vel_grad(elem_gid, 2, 2);

// vel = [u,v,w]
// [du/dx, du/dy, du/dz]
// vel_grad = [dv/dx, dv/dy, dv/dz]
// [dw/dx, dw/dy, dw/dz]
double curl[3];
curl[0] = vel_grad(2, 1) - vel_grad(1, 2); // dw/dy - dv/dz
curl[1] = vel_grad(0, 2) - vel_grad(2, 0); // du/dz - dw/dx
curl[2] = vel_grad(1, 0) - vel_grad(0, 1); // dv/dx - du/dy
curl[0] = GaussPoints_vel_grad(elem_gid, 2, 1) - GaussPoints_vel_grad(elem_gid, 1, 2); // dw/dy - dv/dz
curl[1] = GaussPoints_vel_grad(elem_gid, 0, 2) - GaussPoints_vel_grad(elem_gid, 2, 0); // du/dz - dw/dx
curl[2] = GaussPoints_vel_grad(elem_gid, 1, 0) - GaussPoints_vel_grad(elem_gid, 0, 1); // dv/dx - du/dy

double mag_curl = sqrt(curl[0] * curl[0] + curl[1] * curl[1] + curl[2] * curl[2]);

Expand Down Expand Up @@ -336,7 +330,9 @@ void SGH3D::get_force(const Material_t& Materials,
size_t neighbor_gid = mesh.elems_in_elem(elem_gid, elem_lid);

// calculate the velocity divergence in neighbor
double div_neighbor = GaussPoints_div(neighbor_gid);
double div_neighbor = GaussPoints_vel_grad(neighbor_gid, 0, 0) +
GaussPoints_vel_grad(neighbor_gid, 1, 1) +
GaussPoints_vel_grad(neighbor_gid, 2, 2);

r_face = r_coef * (div_neighbor + small) / (div + small);

Expand Down
155 changes: 83 additions & 72 deletions single-node-refactor/src/Solvers/SGH_solver_3D/src/momentum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,92 +93,103 @@ void SGH3D::update_velocity(double rk_alpha,
/// \brief This function calculates the velocity gradient
///
/// \param Velocity gradient
/// \param Global ids of the nodes in this element
/// \param View of the nodal velocity data
/// \param The finite element B matrix
/// \param mesh object
/// \param node coordinates
/// \param node velocity
/// \param The volume of the particular element
/// \param The global id of this particular element
///
/////////////////////////////////////////////////////////////////////////////
KOKKOS_FUNCTION
void SGH3D::get_velgrad(ViewCArrayKokkos<double>& vel_grad,
const ViewCArrayKokkos<size_t>& elem_node_gids,
const DCArrayKokkos<double>& node_vel,
const ViewCArrayKokkos<double>& b_matrix,
const double elem_vol,
const size_t elem_gid) const
void SGH3D::get_velgrad(DCArrayKokkos<double>& vel_grad,
const Mesh_t mesh,
const DCArrayKokkos<double>& node_coords,
const DCArrayKokkos<double>& node_vel,
const DCArrayKokkos<double>& elem_vol) const
{
const size_t num_nodes_in_elem = 8;
const size_t num_dims = 3;

double u_array[num_nodes_in_elem];
double v_array[num_nodes_in_elem];
double w_array[num_nodes_in_elem];
ViewCArrayKokkos<double> u(u_array, num_nodes_in_elem); // x-direction velocity component
ViewCArrayKokkos<double> v(v_array, num_nodes_in_elem); // y-direction velocity component
ViewCArrayKokkos<double> w(w_array, num_nodes_in_elem); // z-direction velocity component

// get the vertex velocities for the cell
for (size_t node_lid = 0; node_lid < num_nodes_in_elem; node_lid++) {
// Get node gid
size_t node_gid = elem_node_gids(node_lid);

u(node_lid) = node_vel(1, node_gid, 0);
v(node_lid) = node_vel(1, node_gid, 1);
w(node_lid) = node_vel(1, node_gid, 2);
} // end for
// --- calculate the forces acting on the nodes from the element ---
FOR_ALL(elem_gid, 0, mesh.num_elems, {

// --- calculate the velocity gradient terms ---
double inverse_vol = 1.0 / elem_vol;

// x-dir
vel_grad(0, 0) = (u(0) * b_matrix(0, 0) + u(1) * b_matrix(1, 0)
+ u(2) * b_matrix(2, 0) + u(3) * b_matrix(3, 0)
+ u(4) * b_matrix(4, 0) + u(5) * b_matrix(5, 0)
+ u(6) * b_matrix(6, 0) + u(7) * b_matrix(7, 0)) * inverse_vol;
double u_array[num_nodes_in_elem];
double v_array[num_nodes_in_elem];
double w_array[num_nodes_in_elem];
ViewCArrayKokkos<double> u(u_array, num_nodes_in_elem); // x-dir vel component
ViewCArrayKokkos<double> v(v_array, num_nodes_in_elem); // y-dir vel component
ViewCArrayKokkos<double> w(w_array, num_nodes_in_elem); // z-dir vel component

// cut out the node_gids for this element
ViewCArrayKokkos<size_t> elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem);

// The b_matrix are the outward corner area normals
double b_matrix_array[24];
ViewCArrayKokkos<double> b_matrix(b_matrix_array, num_nodes_in_elem, num_dims);
geometry::get_bmatrix(b_matrix, elem_gid, node_coords, elem_node_gids);

vel_grad(0, 1) = (u(0) * b_matrix(0, 1) + u(1) * b_matrix(1, 1)
+ u(2) * b_matrix(2, 1) + u(3) * b_matrix(3, 1)
+ u(4) * b_matrix(4, 1) + u(5) * b_matrix(5, 1)
+ u(6) * b_matrix(6, 1) + u(7) * b_matrix(7, 1)) * inverse_vol;
// get the vertex velocities for the elem
for (size_t node_lid = 0; node_lid < num_nodes_in_elem; node_lid++) {
// Get node gid
size_t node_gid = elem_node_gids(node_lid);

vel_grad(0, 2) = (u(0) * b_matrix(0, 2) + u(1) * b_matrix(1, 2)
+ u(2) * b_matrix(2, 2) + u(3) * b_matrix(3, 2)
+ u(4) * b_matrix(4, 2) + u(5) * b_matrix(5, 2)
+ u(6) * b_matrix(6, 2) + u(7) * b_matrix(7, 2)) * inverse_vol;
u(node_lid) = node_vel(1, node_gid, 0);
v(node_lid) = node_vel(1, node_gid, 1);
w(node_lid) = node_vel(1, node_gid, 2);
} // end for

// y-dir
vel_grad(1, 0) = (v(0) * b_matrix(0, 0) + v(1) * b_matrix(1, 0)
+ v(2) * b_matrix(2, 0) + v(3) * b_matrix(3, 0)
+ v(4) * b_matrix(4, 0) + v(5) * b_matrix(5, 0)
+ v(6) * b_matrix(6, 0) + v(7) * b_matrix(7, 0)) * inverse_vol;
// --- calculate the velocity gradient terms ---
double inverse_vol = 1.0 / elem_vol(elem_gid);

vel_grad(1, 1) = (v(0) * b_matrix(0, 1) + v(1) * b_matrix(1, 1)
+ v(2) * b_matrix(2, 1) + v(3) * b_matrix(3, 1)
+ v(4) * b_matrix(4, 1) + v(5) * b_matrix(5, 1)
+ v(6) * b_matrix(6, 1) + v(7) * b_matrix(7, 1)) * inverse_vol;
vel_grad(1, 2) = (v(0) * b_matrix(0, 2) + v(1) * b_matrix(1, 2)
+ v(2) * b_matrix(2, 2) + v(3) * b_matrix(3, 2)
+ v(4) * b_matrix(4, 2) + v(5) * b_matrix(5, 2)
+ v(6) * b_matrix(6, 2) + v(7) * b_matrix(7, 2)) * inverse_vol;

// z-dir
vel_grad(2, 0) = (w(0) * b_matrix(0, 0) + w(1) * b_matrix(1, 0)
+ w(2) * b_matrix(2, 0) + w(3) * b_matrix(3, 0)
+ w(4) * b_matrix(4, 0) + w(5) * b_matrix(5, 0)
+ w(6) * b_matrix(6, 0) + w(7) * b_matrix(7, 0)) * inverse_vol;

vel_grad(2, 1) = (w(0) * b_matrix(0, 1) + w(1) * b_matrix(1, 1)
+ w(2) * b_matrix(2, 1) + w(3) * b_matrix(3, 1)
+ w(4) * b_matrix(4, 1) + w(5) * b_matrix(5, 1)
+ w(6) * b_matrix(6, 1) + w(7) * b_matrix(7, 1)) * inverse_vol;

vel_grad(2, 2) = (w(0) * b_matrix(0, 2) + w(1) * b_matrix(1, 2)
+ w(2) * b_matrix(2, 2) + w(3) * b_matrix(3, 2)
+ w(4) * b_matrix(4, 2) + w(5) * b_matrix(5, 2)
+ w(6) * b_matrix(6, 2) + w(7) * b_matrix(7, 2)) * inverse_vol;
// x-dir
vel_grad(elem_gid, 0, 0) = (u(0) * b_matrix(0, 0) + u(1) * b_matrix(1, 0)
+ u(2) * b_matrix(2, 0) + u(3) * b_matrix(3, 0)
+ u(4) * b_matrix(4, 0) + u(5) * b_matrix(5, 0)
+ u(6) * b_matrix(6, 0) + u(7) * b_matrix(7, 0)) * inverse_vol;

vel_grad(elem_gid, 0, 1) = (u(0) * b_matrix(0, 1) + u(1) * b_matrix(1, 1)
+ u(2) * b_matrix(2, 1) + u(3) * b_matrix(3, 1)
+ u(4) * b_matrix(4, 1) + u(5) * b_matrix(5, 1)
+ u(6) * b_matrix(6, 1) + u(7) * b_matrix(7, 1)) * inverse_vol;

vel_grad(elem_gid, 0, 2) = (u(0) * b_matrix(0, 2) + u(1) * b_matrix(1, 2)
+ u(2) * b_matrix(2, 2) + u(3) * b_matrix(3, 2)
+ u(4) * b_matrix(4, 2) + u(5) * b_matrix(5, 2)
+ u(6) * b_matrix(6, 2) + u(7) * b_matrix(7, 2)) * inverse_vol;

// y-dir
vel_grad(elem_gid, 1, 0) = (v(0) * b_matrix(0, 0) + v(1) * b_matrix(1, 0)
+ v(2) * b_matrix(2, 0) + v(3) * b_matrix(3, 0)
+ v(4) * b_matrix(4, 0) + v(5) * b_matrix(5, 0)
+ v(6) * b_matrix(6, 0) + v(7) * b_matrix(7, 0)) * inverse_vol;

vel_grad(elem_gid, 1, 1) = (v(0) * b_matrix(0, 1) + v(1) * b_matrix(1, 1)
+ v(2) * b_matrix(2, 1) + v(3) * b_matrix(3, 1)
+ v(4) * b_matrix(4, 1) + v(5) * b_matrix(5, 1)
+ v(6) * b_matrix(6, 1) + v(7) * b_matrix(7, 1)) * inverse_vol;
vel_grad(elem_gid, 1, 2) = (v(0) * b_matrix(0, 2) + v(1) * b_matrix(1, 2)
+ v(2) * b_matrix(2, 2) + v(3) * b_matrix(3, 2)
+ v(4) * b_matrix(4, 2) + v(5) * b_matrix(5, 2)
+ v(6) * b_matrix(6, 2) + v(7) * b_matrix(7, 2)) * inverse_vol;

// z-dir
vel_grad(elem_gid, 2, 0) = (w(0) * b_matrix(0, 0) + w(1) * b_matrix(1, 0)
+ w(2) * b_matrix(2, 0) + w(3) * b_matrix(3, 0)
+ w(4) * b_matrix(4, 0) + w(5) * b_matrix(5, 0)
+ w(6) * b_matrix(6, 0) + w(7) * b_matrix(7, 0)) * inverse_vol;

vel_grad(elem_gid, 2, 1) = (w(0) * b_matrix(0, 1) + w(1) * b_matrix(1, 1)
+ w(2) * b_matrix(2, 1) + w(3) * b_matrix(3, 1)
+ w(4) * b_matrix(4, 1) + w(5) * b_matrix(5, 1)
+ w(6) * b_matrix(6, 1) + w(7) * b_matrix(7, 1)) * inverse_vol;

vel_grad(elem_gid, 2, 2) = (w(0) * b_matrix(0, 2) + w(1) * b_matrix(1, 2)
+ w(2) * b_matrix(2, 2) + w(3) * b_matrix(3, 2)
+ w(4) * b_matrix(4, 2) + w(5) * b_matrix(5, 2)
+ w(6) * b_matrix(6, 2) + w(7) * b_matrix(7, 2)) * inverse_vol;
}); // end parallel for over elem_gid

return;
} // end function
} // end subroutine


/////////////////////////////////////////////////////////////////////////////
Expand Down
Loading

0 comments on commit a8c4cd2

Please sign in to comment.