Merge pull request #174 from Adrian-Diaz/Parallel-Explicit-Solver

Parallel Solvers update

src/Parallel-Solvers/Implicit-Lagrange/FEA_Physics_Modules/FEA_Module_Elasticity.cpp

@@ -2772,6 +2772,13 @@ void FEA_Module_Elasticity::local_matrix_multiply(int ielem, CArrayKokkos<real_t
     else{
       Element_Modulus = module_params->material.elastic_modulus/unit_scaling/unit_scaling;
       Poisson_Ratio = module_params->material.poisson_ratio;
+      if(anisotropic_lattice){
+        for(int idim=0 ; idim < num_dim; idim++){
+          Elastic_Moduli[idim] = module_params->material.elastic_moduli[idim];
+          Shear_Moduli[idim] = module_params->material.shear_moduli[idim];
+          Poisson_Ratios[idim] = module_params->material.poisson_ratios[idim];
+        }
+      }
     }
 
     if(anisotropic_lattice){

src/Parallel-Solvers/Parallel-Explicit/Explicit_Solver.cpp

@@ -90,6 +90,7 @@
 #include "Mass_Constraint.h"
 #include "Moment_of_Inertia_Constraint.h"
 #include "Kinetic_Energy_Minimize.h"
+#include "Area_Normals.h"
 
 #define BUFFER_LINES 20000
 #define MAX_WORD 30
@@ -1402,6 +1403,30 @@ void Explicit_Solver::setup_optimization_problem(){
   problem->finalize(false,true,*fos);
   //problem->check(true,std::cout);
 
+  // //check area normal gradients
+  // ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> rol_debug_coords =
+  // ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(all_node_coords_distributed);
+  // Teuchos::RCP<MV> debug_directions_distributed = Teuchos::rcp(new MV(all_node_map, num_dim));
+  // ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> debug_rol_d =
+  // ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(debug_directions_distributed);
+  // debug_directions_distributed->randomize(-0.00375,0.00375);
+  // if(myrank==0){
+  // for(int inode = 0; inode < 8; inode++){
+  //   for(int idim = 0; idim < num_dim; idim++){
+  //     std::cout<< "AREA NORMALS GRADIENT FOR NODE " << inode << " AND DIM " << idim << std::endl;
+  //     ROL::Ptr<ROL::Objective<real_t>> debug_obj = ROL::makePtr<AreaNormals>(this, inode, idim, nodal_density_flag);
+  //     //debug_directions_distributed->putScalar(-0.1);
+  //     //debug_directions_distributed->norm2(debug_direction_norm);
+  //     //debug_directions_distributed->scale(1/debug_direction_norm(0));
+  //     //set all but first component to 0 for debug
+  //     host_vec_array debug_directions = debug_directions_distributed->getLocalView<HostSpace> (Tpetra::Access::ReadWrite);
+  //     //for(int init = 1; init < nlocal_nodes; init++)
+  //     //directions(4,0) = -0.3;
+  //     debug_obj->checkGradient(*rol_debug_coords, *debug_rol_d);
+  //   }
+  // }
+  // }
+
   //debug checks
   ROL::Ptr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>> rol_x =
    ROL::makePtr<ROL::TpetraMultiVector<real_t,LO,GO,node_type>>(design_node_densities_distributed);

src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/FEA_Module_SGH.cpp

@@ -1419,7 +1419,7 @@ void FEA_Module_SGH::setup(){
     //initialize if topology optimization is used
     if(simparam->topology_optimization_on || simparam->shape_optimization_on){
       init_assembly();
-      assemble_matrix();
+      //assemble_matrix();
     }
 
     // update host copies of arrays modified in this function

src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/FEA_Module_SGH.h

@@ -604,8 +604,8 @@ class FEA_Module_SGH: public FEA_Module{
   //Local FEA data
   DCArrayKokkos<size_t, array_layout, device_type, memory_traits> Global_Gradient_Matrix_Assembly_Map; //Maps element local nodes to columns on ragged right node connectivity graph
   DCArrayKokkos<size_t, array_layout, device_type, memory_traits> Element_Gradient_Matrix_Assembly_Map; //Maps element-node pair to columns on ragged right node to element connectivity
-  RaggedRightArrayKokkos<GO, array_layout, device_type, memory_traits> Graph_Matrix; //stores global indices
-  RaggedRightArrayKokkos<GO, array_layout, device_type, memory_traits> DOF_Graph_Matrix; //stores global indices
+  RaggedRightArrayKokkos<LO, array_layout, device_type, memory_traits> Graph_Matrix; //stores local indices
+  RaggedRightArrayKokkos<LO, array_layout, device_type, memory_traits> DOF_Graph_Matrix; //stores local indices
   RaggedRightArrayKokkos<real_t, Kokkos::LayoutRight, device_type, memory_traits, array_layout> Force_Gradient_Positions;
   RaggedRightArrayKokkos<real_t, Kokkos::LayoutRight, device_type, memory_traits, array_layout> Force_Gradient_Velocities;
   CArrayKokkos<real_t, array_layout, device_type, memory_traits> Force_Gradient_Energies; //transposed such that elem ids correspond to rows

src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/force_gradients_sgh.cpp

@@ -269,7 +269,6 @@ void FEA_Module_SGH::get_force_vgradient_sgh(const DCArrayKokkos <material_t> &m
                 muc(node_lid) = elem_den(elem_gid) *
                                (material(mat_id).q1ex*elem_sspd(elem_gid) + material(mat_id).q2ex*mag_vel);
             } // end if on divergence sign
-
 
             size_t use_shock_dir = 0;
             double mu_term;
@@ -940,7 +939,6 @@ void FEA_Module_SGH::get_force_ugradient_sgh(const DCArrayKokkos <material_t> &m
       }
     }); // end parallel for loop over nodes
     Kokkos::fence();
-
     // --- calculate the forces acting on the nodes from the element ---
     for (size_t elem_gid = 0; elem_gid < rnum_elem; elem_gid++){
     //FOR_ALL_CLASS (elem_gid, 0, rnum_elem, {
@@ -964,7 +962,7 @@ void FEA_Module_SGH::get_force_ugradient_sgh(const DCArrayKokkos <material_t> &m
 
         // the sums in the Riemann solver
         double sum_array[4];
-        double sum_gradient_array[4*max_nodes_per_element*num_dim];
+        double sum_gradient_array[4*max_nodes_per_element*num_dims];
 
         // corner shock impeadance x |corner area normal dot shock_dir|
         double muc_array[8];
@@ -980,7 +978,7 @@ void FEA_Module_SGH::get_force_ugradient_sgh(const DCArrayKokkos <material_t> &m
         // --- Create views of arrays to aid the force calculation ---
 
         ViewCArrayKokkos <double> tau(tau_array, num_dims, num_dims);
-        ViewCArrayKokkos <double> tau_gradient(tau_gradient_array, num_dims, num_dims,max_nodes_per_element, num_dims);
+        ViewCArrayKokkos <double> tau_gradient(tau_gradient_array, num_dims, num_dims, max_nodes_per_element, num_dims);
         ViewCArrayKokkos <double> volume_gradients(volume_gradients_array, num_nodes_in_elem, num_dims);
         ViewCArrayKokkos <double> area_normal(area_normal_array, num_nodes_in_elem, num_dims);
         ViewCArrayKokkos <double> area_normal_gradients(area_normal_gradients_array, num_nodes_in_elem, num_dims, num_nodes_in_elem, num_dims);

src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/geometry.cpp

@@ -1514,7 +1514,7 @@ void FEA_Module_SGH::get_bmatrix_gradients(const ViewCArrayKokkos <double> &B_ma
     gradient_terms(0,7) = y(0) + y(3) - y(5) - y(6);
 
     //y derivative
-    gradient_terms(1,0) = -x(1) + x(3) - x(4) + x(7);
+    gradient_terms(1,0) = -x(1) + x(3) - x(5) + x(7);
     gradient_terms(1,1) = x(0) - x(5);
     gradient_terms(1,2) = 0;
     gradient_terms(1,3) = -x(0) + x(7);