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

Parallel solvers update

src/Parallel-Solvers/CMakeLists.txt

@@ -1,5 +1,3 @@
-cmake_minimum_required(VERSION 3.17)
-project(parallel-src)
 
 find_package(Vector)
 if (CMAKE_VECTOR_NOVEC)

src/Parallel-Solvers/FEA_Module_Inertial.cpp

@@ -663,6 +663,239 @@ void FEA_Module_Inertial::compute_nodal_gradients(const_host_vec_array design_va
 
 void FEA_Module_Inertial::compute_shape_gradients(const_host_vec_array design_coords, host_vec_array design_gradients, bool use_initial_coords)
 {
+    // local number of uniquely assigned elements
+    size_t nonoverlap_nelements = element_map->getLocalNumElements();
+    // local variable for host view in the dual view
+    const_host_vec_array all_node_coords;
+
+    all_node_coords = all_design_node_coords_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    const_host_elem_conn_array nodes_in_elem = global_nodes_in_elem_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    const_host_vec_array       all_node_densities;
+    if (nodal_density_flag)
+    {
+        all_node_densities = all_node_densities_distributed->getLocalView<HostSpace>(Tpetra::Access::ReadOnly);
+    }
+    int    nodes_per_elem = elem->num_basis();
+    int    z_quad, y_quad, x_quad, direct_product_count;
+    size_t local_node_id;
+    LO     ielem;
+    GO     global_element_index;
+
+    real_t Jacobian, weight_multiply, current_density;
+    // CArrayKokkos<real_t> legendre_nodes_1D(num_gauss_points);
+    // CArrayKokkos<real_t> legendre_weights_1D(num_gauss_points);
+    CArray<real_t> legendre_nodes_1D(num_gauss_points);
+    CArray<real_t> legendre_weights_1D(num_gauss_points);
+    real_t pointer_quad_coordinate[num_dim];
+    real_t pointer_quad_coordinate_weight[num_dim];
+    real_t pointer_interpolated_point[num_dim];
+    real_t pointer_JT_row1[num_dim];
+    real_t pointer_JT_row2[num_dim];
+    real_t pointer_JT_row3[num_dim];
+    ViewCArray<real_t> quad_coordinate(pointer_quad_coordinate, num_dim);
+    ViewCArray<real_t> quad_coordinate_weight(pointer_quad_coordinate_weight, num_dim);
+    ViewCArray<real_t> interpolated_point(pointer_interpolated_point, num_dim);
+    ViewCArray<real_t> JT_row1(pointer_JT_row1, num_dim);
+    ViewCArray<real_t> JT_row2(pointer_JT_row2, num_dim);
+    ViewCArray<real_t> JT_row3(pointer_JT_row3, num_dim);
+
+    real_t pointer_basis_values[elem->num_basis()];
+    real_t pointer_basis_derivative_s1[elem->num_basis()];
+    real_t pointer_basis_derivative_s2[elem->num_basis()];
+    real_t pointer_basis_derivative_s3[elem->num_basis()];
+    ViewCArray<real_t> basis_values(pointer_basis_values, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s1(pointer_basis_derivative_s1, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s2(pointer_basis_derivative_s2, elem->num_basis());
+    ViewCArray<real_t> basis_derivative_s3(pointer_basis_derivative_s3, elem->num_basis());
+    CArrayKokkos<real_t, array_layout, device_type, memory_traits> nodal_positions(elem->num_basis(), num_dim);
+    CArrayKokkos<real_t, array_layout, device_type, memory_traits> nodal_density(elem->num_basis());
+    CArrayKokkos<real_t, array_layout, device_type, memory_traits> jacobian_derivatives(elem->num_basis(), num_dim);
+
+    // initialize weights
+    elements::legendre_nodes_1D(legendre_nodes_1D, num_gauss_points);
+    elements::legendre_weights_1D(legendre_weights_1D, num_gauss_points);
+
+    Solver::node_ordering_convention                             active_node_ordering_convention = Solver_Pointer_->active_node_ordering_convention;
+    CArrayKokkos<size_t, array_layout, HostSpace, memory_traits> convert_node_order(max_nodes_per_element);
+    if ((active_node_ordering_convention == Solver::ENSIGHT && num_dim == 3) || (active_node_ordering_convention == Solver::IJK && num_dim == 2))
+    {
+        convert_node_order(0) = 0;
+        convert_node_order(1) = 1;
+        convert_node_order(2) = 3;
+        convert_node_order(3) = 2;
+        if (num_dim == 3)
+        {
+            convert_node_order(4) = 4;
+            convert_node_order(5) = 5;
+            convert_node_order(6) = 7;
+            convert_node_order(7) = 6;
+        }
+    }
+    else
+    {
+        convert_node_order(0) = 0;
+        convert_node_order(1) = 1;
+        convert_node_order(2) = 2;
+        convert_node_order(3) = 3;
+        if (num_dim == 3)
+        {
+            convert_node_order(4) = 4;
+            convert_node_order(5) = 5;
+            convert_node_order(6) = 6;
+            convert_node_order(7) = 7;
+        }
+    }
+
+    // initialize design gradients to 0
+    for (int init = 0; init < nlocal_nodes; init++)
+    {
+        design_gradients(init, 0) = design_gradients(init, 1) = design_gradients(init, 2) = 0;
+    }
+
+    // loop over elements and use quadrature rule to compute volume from Jacobian determinant
+    for (int ielem = 0; ielem < rnum_elem; ielem++)
+    {
+        // acquire set of nodes for this local element
+        for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+        {
+            local_node_id = all_node_map->getLocalElement(nodes_in_elem(ielem, convert_node_order(node_loop)));
+            nodal_positions(node_loop, 0) = all_node_coords(local_node_id, 0);
+            nodal_positions(node_loop, 1) = all_node_coords(local_node_id, 1);
+            nodal_positions(node_loop, 2) = all_node_coords(local_node_id, 2);
+            if (nodal_density_flag)
+            {
+                nodal_density(node_loop) = all_node_densities(local_node_id, 0);
+            }
+
+            //init derivatives
+            jacobian_derivatives(node_loop,0) = 0;
+            jacobian_derivatives(node_loop,1) = 0;
+            if(num_dim==3){
+                jacobian_derivatives(node_loop,2) = 0;
+            }
+        }
+
+        direct_product_count = std::pow(num_gauss_points, num_dim);
+
+        // loop over quadrature points
+        for (int iquad = 0; iquad < direct_product_count; iquad++)
+        {
+            // set current quadrature point
+            if (num_dim == 3)
+            {
+                z_quad = iquad / (num_gauss_points * num_gauss_points);
+            }
+            y_quad = (iquad % (num_gauss_points * num_gauss_points)) / num_gauss_points;
+            x_quad = iquad % num_gauss_points;
+            quad_coordinate(0) = legendre_nodes_1D(x_quad);
+            quad_coordinate(1) = legendre_nodes_1D(y_quad);
+            if (num_dim == 3)
+            {
+                quad_coordinate(2) = legendre_nodes_1D(z_quad);
+            }
+
+            // set current quadrature weight
+            quad_coordinate_weight(0) = legendre_weights_1D(x_quad);
+            quad_coordinate_weight(1) = legendre_weights_1D(y_quad);
+            if (num_dim == 3)
+            {
+                quad_coordinate_weight(2) = legendre_weights_1D(z_quad);
+            }
+            else
+            {
+                quad_coordinate_weight(2) = 1;
+            }
+            weight_multiply = quad_coordinate_weight(0) * quad_coordinate_weight(1) * quad_coordinate_weight(2);
+
+            // compute shape functions at this point for the element type
+            elem->basis(basis_values, quad_coordinate);
+
+            // compute all the necessary coordinates and derivatives at this point
+
+            // compute shape function derivatives
+            elem->partial_xi_basis(basis_derivative_s1, quad_coordinate);
+            elem->partial_eta_basis(basis_derivative_s2, quad_coordinate);
+            elem->partial_mu_basis(basis_derivative_s3, quad_coordinate);
+
+            // compute derivatives of x,y,z w.r.t the s,t,w isoparametric space needed by JT (Transpose of the Jacobian)
+            // derivative of x,y,z w.r.t s
+            JT_row1(0) = 0;
+            JT_row1(1) = 0;
+            JT_row1(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row1(0) += nodal_positions(node_loop, 0) * basis_derivative_s1(node_loop);
+                JT_row1(1) += nodal_positions(node_loop, 1) * basis_derivative_s1(node_loop);
+                JT_row1(2) += nodal_positions(node_loop, 2) * basis_derivative_s1(node_loop);
+            }
+
+            // derivative of x,y,z w.r.t t
+            JT_row2(0) = 0;
+            JT_row2(1) = 0;
+            JT_row2(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row2(0) += nodal_positions(node_loop, 0) * basis_derivative_s2(node_loop);
+                JT_row2(1) += nodal_positions(node_loop, 1) * basis_derivative_s2(node_loop);
+                JT_row2(2) += nodal_positions(node_loop, 2) * basis_derivative_s2(node_loop);
+            }
+
+            // derivative of x,y,z w.r.t w
+            JT_row3(0) = 0;
+            JT_row3(1) = 0;
+            JT_row3(2) = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                JT_row3(0) += nodal_positions(node_loop, 0) * basis_derivative_s3(node_loop);
+                JT_row3(1) += nodal_positions(node_loop, 1) * basis_derivative_s3(node_loop);
+                JT_row3(2) += nodal_positions(node_loop, 2) * basis_derivative_s3(node_loop);
+                // debug print
+                /*if(myrank==1&&nodal_positions(node_loop,2)*basis_derivative_s3(node_loop)<-10000000){
+                std::cout << " ELEMENT VOLUME JACOBIAN DEBUG ON TASK " << myrank << std::endl;
+                std::cout << node_loop+1 << " " << JT_row3(2) << " "<< nodal_positions(node_loop,2) <<" "<< basis_derivative_s3(node_loop) << std::endl;
+                std::fflush(stdout);
+                }*/
+            }
+
+            // compute the determinant of the Jacobian
+            Jacobian = JT_row1(0) * (JT_row2(1) * JT_row3(2) - JT_row3(1) * JT_row2(2)) -
+                       JT_row1(1) * (JT_row2(0) * JT_row3(2) - JT_row3(0) * JT_row2(2)) +
+                       JT_row1(2) * (JT_row2(0) * JT_row3(1) - JT_row3(0) * JT_row2(1));
+            if (Jacobian < 0)
+            {
+                Jacobian = -Jacobian;
+            }
+
+            current_density = 0;
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                current_density += nodal_density(node_loop) * basis_values(node_loop);
+                jacobian_derivatives(node_loop,0) = basis_derivative_s1(node_loop) * (JT_row2(1) * JT_row3(2) - JT_row3(1) * JT_row2(2)) -
+                                                    JT_row1(1) * (basis_derivative_s2(node_loop) * JT_row3(2) - basis_derivative_s3(node_loop) * JT_row2(2)) +
+                                                    JT_row1(2) * (basis_derivative_s2(node_loop) * JT_row3(1) - basis_derivative_s3(node_loop) * JT_row2(1));
+
+                jacobian_derivatives(node_loop,1) = JT_row1(0) * (basis_derivative_s2(node_loop) * JT_row3(2) - basis_derivative_s3(node_loop) * JT_row2(2)) -
+                                                    basis_derivative_s1(node_loop) * (JT_row2(0) * JT_row3(2) - JT_row3(0) * JT_row2(2)) +
+                                                    JT_row1(2) * (JT_row2(0) * basis_derivative_s3(node_loop) - JT_row3(0) * basis_derivative_s2(node_loop));
+
+                jacobian_derivatives(node_loop,2) = JT_row1(0) * (JT_row2(1) * basis_derivative_s3(node_loop) - JT_row3(1) * basis_derivative_s2(node_loop)) -
+                                                    JT_row1(1) * (JT_row2(0) * basis_derivative_s3(node_loop) - JT_row3(0) * basis_derivative_s2(node_loop)) +
+                                                    basis_derivative_s1(node_loop) * (JT_row2(0) * JT_row3(1) - JT_row3(0) * JT_row2(1));
+            }
+
+            // assign contribution to every local node this element has
+            for (int node_loop = 0; node_loop < elem->num_basis(); node_loop++)
+            {
+                if (map->isNodeGlobalElement(nodes_in_elem(ielem, node_loop)))
+                {
+                    local_node_id = map->getLocalElement(nodes_in_elem(ielem, node_loop));
+                    design_gradients(local_node_id, 0) += weight_multiply * current_density * jacobian_derivatives(node_loop,0);
+                    design_gradients(local_node_id, 1) += weight_multiply * current_density * jacobian_derivatives(node_loop,1);
+                    design_gradients(local_node_id, 2) += weight_multiply * current_density * jacobian_derivatives(node_loop,2);
+                }
+            }
+        }
+    }
 
 }
 

src/Parallel-Solvers/Parallel-Explicit/CMakeLists.txt

@@ -1,9 +1,3 @@
-cmake_minimum_required(VERSION 3.1.3)
-
-set(CMAKE_CXX_FLAGS  "${Trilinos_CXX_COMPILER_FLAGS} ${CMAKE_CXX_FLAGS}")
-set(CMAKE_C_FLAGS  "${Trilinos_C_COMPILER_FLAGS} ${CMAKE_C_FLAGS}")
-set(CMAKE_Fortran_FLAGS  "${Trilinos_Fortran_COMPILER_FLAGS} ${CMAKE_Fortran_FLAGS}")
-
 include_directories(SGH_Solver/include)
 include_directories(Dynamic_Elastic_Solver)
 #include_directories(Eulerian_Solver)

src/Parallel-Solvers/Parallel-Explicit/Dynamic_Elastic_Solver/CMakeLists.txt

@@ -1,8 +1,3 @@
-cmake_minimum_required(VERSION 3.1.3)
-
-set(CMAKE_CXX_FLAGS  "${Trilinos_CXX_COMPILER_FLAGS} ${CMAKE_CXX_FLAGS}")
-set(CMAKE_C_FLAGS  "${Trilinos_C_COMPILER_FLAGS} ${CMAKE_C_FLAGS}")
-set(CMAKE_Fortran_FLAGS  "${Trilinos_Fortran_COMPILER_FLAGS} ${CMAKE_Fortran_FLAGS}")
 set(SRC_Files geometry.cpp boundary.cpp time_integration.cpp momentum.cpp
 elastic_optimization.cpp force_gradients_elastic.cpp force_elastic.cpp properties.cpp)
 set(FEA_Module_SRC FEA_Module_Dynamic_Elasticity.cpp )

src/Parallel-Solvers/Parallel-Explicit/SGH_Solver/CMakeLists.txt

@@ -1,21 +1,16 @@
-cmake_minimum_required(VERSION 3.1.3)
-
-set(CMAKE_CXX_FLAGS  "${Trilinos_CXX_COMPILER_FLAGS} ${CMAKE_CXX_FLAGS}")
-set(CMAKE_C_FLAGS  "${Trilinos_C_COMPILER_FLAGS} ${CMAKE_C_FLAGS}")
-set(CMAKE_Fortran_FLAGS  "${Trilinos_Fortran_COMPILER_FLAGS} ${CMAKE_Fortran_FLAGS}")
-
 include_directories(include)
 
-set(SRC_Files 
-src/geometry.cpp 
-src/boundary.cpp 
-src/time_integration.cpp 
-src/momentum.cpp 
+set(SRC_Files
+src/geometry.cpp
+src/boundary.cpp
+src/time_integration.cpp
+src/momentum.cpp
 src/force_sgh.cpp
-src/sgh_optimization.cpp 
-src/force_gradients_sgh.cpp 
-src/power_gradients_sgh.cpp 
-src/energy_sgh.cpp 
+src/sgh_optimization.cpp
+src/sgh_shape_optimization.cpp
+src/force_gradients_sgh.cpp
+src/power_gradients_sgh.cpp
+src/energy_sgh.cpp
 src/properties.cpp
 src/setup_sgh.cpp)
 

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

@@ -424,6 +424,10 @@ class FEA_Module_SGH : public FEA_Module
 
     virtual void update_forward_solve(Teuchos::RCP<const MV> zp, bool print_design=false);
 
+    void update_forward_solve_TO(Teuchos::RCP<const MV> zp);
+
+    void update_forward_solve_SO(Teuchos::RCP<const MV> zp);
+
     void comm_node_masses();
 
     void comm_adjoint_vector(int cycle);
@@ -500,6 +504,8 @@ class FEA_Module_SGH : public FEA_Module
 
     void compute_topology_optimization_adjoint_full(Teuchos::RCP<const MV> design_densities_distributed); // Force depends on node coords, velocity, and sie
 
+    void compute_shape_optimization_adjoint_full(Teuchos::RCP<const MV> design_densities_distributed); // Force depends on node coords, velocity, and sie
+
     void compute_topology_optimization_gradient_full(Teuchos::RCP<const MV> design_densities_distributed, Teuchos::RCP<MV> design_gradients_distributed);
 
     void compute_topology_optimization_gradient_tally(Teuchos::RCP<const MV> design_densities_distributed, Teuchos::RCP<MV> design_gradients_distributed,
@@ -508,15 +514,13 @@ class FEA_Module_SGH : public FEA_Module
     void compute_topology_optimization_gradient_IVP(Teuchos::RCP<const MV> design_densities_distributed, Teuchos::RCP<MV> design_gradients_distributed,
                                                       unsigned long cycle, real_t global_dt);
 
-    void compute_shape_optimization_adjoint_full(Teuchos::RCP<const MV> design_coordinates_distributed){} // Force depends on node coords, velocity, and sie
-
-    void compute_shape_optimization_gradient_full(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed){}
+    void compute_shape_optimization_gradient_full(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed);
 
     void compute_shape_optimization_gradient_tally(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed,
-                                                      unsigned long cycle, real_t global_dt){}
+                                                      unsigned long cycle, real_t global_dt);
 
     void compute_shape_optimization_gradient_IVP(Teuchos::RCP<const MV> design_coordinates_distributed, Teuchos::RCP<MV> design_gradients_distributed,
-                                                      unsigned long cycle, real_t global_dt){}
+                                                      unsigned long cycle, real_t global_dt);
 
     void boundary_adjoint(const mesh_t& mesh,
                           const DCArrayKokkos<boundary_t>& boundary,