tpetra

CMakeLists.txt

@@ -439,6 +439,8 @@ set(GOMA_INCLUDES
     include/el_quality.h
     include/exo_conn.h
     include/exo_struct.h
+    include/linalg/sparse_matrix.h
+    include/linalg/sparse_matrix_tpetra.h
     include/load_field_variables.h
     include/loca_const.h
     include/loca_util_const.h
@@ -591,6 +593,8 @@ set(GOMA_SOURCES
     src/el_quality.c
     src/exo_conn.c
     src/globals.c
+    src/linalg/sparse_matrix.cpp
+    src/linalg/sparse_matrix_tpetra.cpp
     src/load_field_variables.c
     src/loca_bord.c
     src/loca_eigen_c2f.F

include/dp_comm.h

@@ -35,6 +35,16 @@ EXTERN void exchange_dof(Comm_Ex *, /* cx - ptr to communications exchange info
                          double *,  /* x - local processor dof-based vector */
                          int);
 
+EXTERN void exchange_dof_int(Comm_Ex *, /* cx - ptr to communications exchange info */
+                         Dpi *,     /* dpi - distributed processing info */
+                         int *,  /* x - local processor dof-based vector */
+                         int);
+
+EXTERN void exchange_dof_long_long(Comm_Ex *, /* cx - ptr to communications exchange info */
+                         Dpi *,     /* dpi - distributed processing info */
+                         long long *,  /* x - local processor dof-based vector */
+                         int);
+
 EXTERN void exchange_node(Comm_Ex *cx, /* cx - ptr to communications exchange info */
                           Dpi *d,      /* dpi - distributed processing info */
                           double *a);  /* x - local processor node-based vector */

include/linalg/sparse_matrix.h

@@ -0,0 +1,108 @@
+#ifndef GOMA_SPARSE_MATRIX
+#define GOMA_SPARSE_MATRIX
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#ifdef GOMA_MATRIX_GO_LONG_LONG
+typedef long long GomaGlobalOrdinal;
+#define MPI_GOMA_ORDINAL MPI_LONG_LONG
+#else
+typedef int GomaGlobalOrdinal;
+#define MPI_GOMA_ORDINAL MPI_INT
+#endif
+
+#ifdef DISABLE_CPP
+#define CPP_ALREADY_DISABLED
+#endif
+#define DISABLE_CPP
+#include "dp_types.h"
+#include "exo_struct.h"
+#include "mm_eh.h"
+#include "rf_node_const.h"
+#include "sl_util_structs.h"
+#ifndef CPP_ALREADY_DISABLED
+#undef DISABLE_CPP
+#endif
+
+
+enum GomaSparseMatrixType {
+  GOMA_SPARSE_MATRIX_TYPE_EPETRA,
+  GOMA_SPARSE_MATRIX_TYPE_AZTEC_MRS,
+  GOMA_SPARSE_MATRIX_TYPE_TPETRA,
+  GOMA_SPARSE_MATRIX_TYPE_PETSC
+};
+
+struct g_GomaSparseMatrix {
+  enum GomaSparseMatrixType type;
+  void *data;
+  GomaGlobalOrdinal *global_ids;
+  GomaGlobalOrdinal n_rows;
+  GomaGlobalOrdinal n_cols;
+  // Create matrix with given rows and columns
+  // coo_rows and coo_cols are the COO representation of the matrix ordered by row
+  // local_nnz is the number of non-zero entries in the local partition
+  // max_nz_per_row is the maximum number of non-zero entries in any row
+  goma_error (*create_graph)(struct g_GomaSparseMatrix *matrix,
+                             GomaGlobalOrdinal n_rows,
+                             GomaGlobalOrdinal *row_list,
+                                      GomaGlobalOrdinal n_cols,
+                                      GomaGlobalOrdinal *col_list,
+                             GomaGlobalOrdinal local_nnz,
+                             GomaGlobalOrdinal max_nz_per_row,
+                             GomaGlobalOrdinal *coo_rows,
+                             GomaGlobalOrdinal *coo_cols);
+  // optional, run after create graph to finalize structure
+  goma_error (*complete_graph)(struct g_GomaSparseMatrix *matrix);
+  // Insert values into matrix row replacing existing values
+  goma_error (*insert_row_values)(struct g_GomaSparseMatrix *matrix,
+                                  GomaGlobalOrdinal global_row,
+                                  GomaGlobalOrdinal num_entries,
+                                  double *values,
+                                  GomaGlobalOrdinal *indices);
+  // Add values into a matrix row, sums existing values
+  goma_error (*sum_into_row_values)(struct g_GomaSparseMatrix *matrix,
+                                    GomaGlobalOrdinal global_row,
+                                    GomaGlobalOrdinal num_entries,
+                                    double *values,
+                                    GomaGlobalOrdinal *indices);
+  // set matrix non-zeros to specified scalar value (commonly re-zero for next assembly)
+  goma_error (*put_scalar)(struct g_GomaSparseMatrix *matrix, double scalar);
+  // row sum scaling, compute row sum scale, scale matrix and b, and return scaling vector
+  goma_error (*row_sum_scaling)(struct g_GomaSparseMatrix *matrix, double *b, double *scale);
+  // delete the allocated matrix;
+  goma_error (*destroy)(struct g_GomaSparseMatrix *matrix);
+};
+typedef struct g_GomaSparseMatrix *GomaSparseMatrix;
+
+goma_error GomaSparseMatrix_Create(GomaSparseMatrix *matrix, enum GomaSparseMatrixType type);
+// populates:
+//   - global_ids
+//   - n_rows
+//   - n_cols
+goma_error GomaSparseMatrix_SetProblemGraph(
+    GomaSparseMatrix matrix,
+    int num_internal_dofs,
+    int num_boundary_dofs,
+    int num_external_dofs,
+    int local_nodes,
+    NODE_INFO_STRUCT **Nodes,
+    int MaxVarPerNode,
+    int *Matilda,
+    int Inter_Mask[MAX_NUM_MATRICES][MAX_VARIABLE_TYPES][MAX_VARIABLE_TYPES],
+    Exo_DB *exo,
+    Dpi *dpi,
+    Comm_Ex *cx,
+    int imtrx,
+    int Debug_Flag,
+    struct GomaLinearSolverData *ams);
+
+goma_error GomaSparseMatrix_LoadLec(GomaSparseMatrix matrix, int ielem, struct Local_Element_Contributions *lec, double resid_vector[]);
+
+goma_error GomaSparseMatrix_Destroy(GomaSparseMatrix *matrix);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // GOMA_SPARSE_MATRIX

include/linalg/sparse_matrix_tpetra.h

@@ -0,0 +1,60 @@
+#ifndef GOMA_SPARSE_MATRIX_TPETRA
+#define GOMA_SPARSE_MATRIX_TPETRA
+#ifdef __cplusplus
+#include "Teuchos_RCP.hpp"
+#include <Tpetra_FECrsMatrix.hpp>
+
+#include "linalg/sparse_matrix.h"
+
+using LO = int;
+using GO = GomaGlobalOrdinal;
+
+struct TpetraSparseMatrix {
+  Teuchos::RCP<Tpetra::FECrsMatrix<double, LO, GO>> matrix;
+  Teuchos::RCP<Tpetra::Map<LO, GO>> row_map;
+  Teuchos::RCP<Tpetra::Map<LO, GO>> col_map;
+  Teuchos::RCP<Tpetra::FECrsGraph<LO, GO>> crs_graph;
+  TpetraSparseMatrix() = default;
+};
+
+extern "C" {
+#endif
+
+
+goma_error GomaSparseMatrix_Tpetra_Create(GomaSparseMatrix *matrix);
+
+goma_error g_tpetra_create_graph(GomaSparseMatrix matrix,
+                           GomaGlobalOrdinal n_rows,
+                           GomaGlobalOrdinal *row_list,
+                                      GomaGlobalOrdinal n_cols,
+                                      GomaGlobalOrdinal *col_list,
+                           GomaGlobalOrdinal local_nnz,
+                           GomaGlobalOrdinal max_per_row,
+                           GomaGlobalOrdinal *coo_rows,
+                           GomaGlobalOrdinal *coo_cols);
+
+goma_error complete_graph(GomaSparseMatrix matrix);
+
+goma_error g_tpetra_insert_row_values(GomaSparseMatrix matrix,
+                                           GomaGlobalOrdinal global_row,
+                                           GomaGlobalOrdinal num_entries,
+                                           double *values,
+                                           GomaGlobalOrdinal *indices);
+
+goma_error g_tpetra_sum_into_row_values(GomaSparseMatrix matrix,
+                                  GomaGlobalOrdinal global_row,
+                                  GomaGlobalOrdinal num_entries,
+                                  double *values,
+                                  GomaGlobalOrdinal *indices);
+
+goma_error g_tpetra_put_scalar(GomaSparseMatrix matrix, double scalar);
+
+goma_error g_tpetra_row_sum_scaling(GomaSparseMatrix matrix, double *b, double *scale);
+
+goma_error g_tpetra_destroy(GomaSparseMatrix matrix);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // GOMA_SPARSE_MATRIX_TPETRA

include/mm_mp_const.h

@@ -563,7 +563,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri
 #define SOLID_DIFFUSION_ELECTRONEUTRALITY        911
 #define SOLID_DIFFUSION                          912
 #define GAS_DIFFUSION                            913
-#define FULL                                     914
+#define METAL_CORROSION_FULL                     914
 #define ANNIHILATION_ELECTRONEUTRALITY           915
 #define ANNIHILATION                             916
 #define NET_CHARGE                               917 /* refer to F multiplied by sum of ci zi   */

include/sl_stratimikos_interface.h

@@ -9,6 +9,12 @@ extern "C" {
 #include "rf_fem_const.h"
 #include "rf_io_const.h"
 
+int stratimikos_solve_tpetra(struct GomaLinearSolverData *ams,
+                      double *x_,
+                      double *b_,
+                      int *iterations,
+                      char stratimikos_file[MAX_NUM_MATRICES][MAX_CHAR_IN_INPUT],
+                      int imtrx);
 int stratimikos_solve(struct GomaLinearSolverData *ams,
                       double *x_,
                       double *b_,

include/sl_util_structs.h

@@ -51,19 +51,6 @@
 #define JAC      0
 #define NUM_ALSS 1
 
-struct Matrix_Data {
-  struct GomaLinearSolverData *ams;
-  double *x;       /* Solution vector */
-  double *x_old;   /* Solution vector , previous last time step */
-  double *x_older; /* Solution vector , previous prev time step */
-  double *x_oldest;
-  double *xdot;     /* xdot of current solution                  */
-  double *xdot_old; /* xdot_old of current solution              */
-  double *xdot_older;
-  double *x_update;     /* last update vector */
-  double *resid_vector; /* Residual vector */
-  double *scale;
-};
 
 struct GomaLinearSolverData {
   int proc_config[AZ_PROC_SIZE];
@@ -128,6 +115,20 @@ struct GomaLinearSolverData {
   int *GlobalIDs;                  /* Pointer to global ids of DOFs (only available with epetra) */
 
   void *PetscMatrixData;
+  void *GomaMatrixData;
 };
 
+struct Matrix_Data {
+  struct GomaLinearSolverData *ams;
+  double *x;       /* Solution vector */
+  double *x_old;   /* Solution vector , previous last time step */
+  double *x_older; /* Solution vector , previous prev time step */
+  double *x_oldest;
+  double *xdot;     /* xdot of current solution                  */
+  double *xdot_old; /* xdot_old of current solution              */
+  double *xdot_older;
+  double *x_update;     /* last update vector */
+  double *resid_vector; /* Residual vector */
+  double *scale;
+};
 #endif

src/dp_comm.c

@@ -100,6 +100,62 @@ void exchange_dof(Comm_Ex *cx, Dpi *dpi, double *x, int imtrx)
   safer_free((void **)&ptr_send_list);
 #endif /* PARALLEL */
 }
+
+void exchange_dof_int(Comm_Ex *cx, Dpi *dpi, int *x, int imtrx)
+
+/************************************************************
+ *
+ *  exchange_dof():
+ *
+ *  send/recv appropriate pieces of a dof-based double array
+ ************************************************************/
+{
+  COMM_NP_STRUCT *np_base, *np_ptr;
+  int *ptr_send_list, *ptr_recv_list;
+  register int *ptrd;
+  register int *ptr_int, i;
+  int p;
+  int num_neighbors = dpi->num_neighbors;
+  int total_num_send_unknowns;
+
+  if (num_neighbors == 0)
+    return;
+
+#ifdef PARALLEL
+  total_num_send_unknowns = ptr_dof_send[imtrx][dpi->num_neighbors];
+  np_base = alloc_struct_1(COMM_NP_STRUCT, dpi->num_neighbors);
+  ptrd = malloc(total_num_send_unknowns * sizeof(int));
+  ptr_send_list = ptrd;
+
+  /*
+   * gather up the list of send unknowns
+   */
+  ptr_int = list_dof_send[imtrx];
+  for (i = total_num_send_unknowns; i > 0; i--) {
+    *ptrd++ = x[*ptr_int++];
+  }
+
+  /*
+   * store base address for the start of the external degrees of freedom
+   * in this vector
+   */
+  ptr_recv_list = x + num_internal_dofs[imtrx] + num_boundary_dofs[imtrx];
+
+  np_ptr = np_base;
+  for (p = 0; p < dpi->num_neighbors; p++) {
+    np_ptr->neighbor_ProcID = cx[p].neighbor_name;
+    np_ptr->send_message_buf = (void *)(ptr_send_list + ptr_dof_send[imtrx][p]);
+    np_ptr->send_message_length = sizeof(int) * cx[p].num_dofs_send;
+    np_ptr->recv_message_buf = (void *)ptr_recv_list;
+    np_ptr->recv_message_length = sizeof(int) * cx[p].num_dofs_recv;
+    ptr_recv_list += cx[p].num_dofs_recv;
+    np_ptr++;
+  }
+  exchange_neighbor_proc_info(dpi->num_neighbors, np_base);
+  safer_free((void **)&np_base);
+  safer_free((void **)&ptr_send_list);
+#endif /* PARALLEL */
+}
 /********************************************************************/
 /********************************************************************/
 /********************************************************************/
@@ -196,6 +252,62 @@ void exchange_node(Comm_Ex *cx, Dpi *dpi, double *x)
   safer_free((void **)&ptr_send_list);
 #endif
 }
+
+void exchange_dof_long_long(Comm_Ex *cx, Dpi *dpi, long long *x, int imtrx)
+
+/************************************************************
+ *
+ *  exchange_dof():
+ *
+ *  send/recv appropriate pieces of a dof-based double array
+ ************************************************************/
+{
+  COMM_NP_STRUCT *np_base, *np_ptr;
+  long long *ptr_send_list, *ptr_recv_list;
+  register long long *ptrd;
+  register int *ptr_int, i;
+  int p;
+  int num_neighbors = dpi->num_neighbors;
+  int total_num_send_unknowns;
+
+  if (num_neighbors == 0)
+    return;
+
+#ifdef PARALLEL
+  total_num_send_unknowns = ptr_dof_send[imtrx][dpi->num_neighbors];
+  np_base = alloc_struct_1(COMM_NP_STRUCT, dpi->num_neighbors);
+  ptrd = malloc(total_num_send_unknowns * sizeof(long long));
+  ptr_send_list = ptrd;
+
+  /*
+   * gather up the list of send unknowns
+   */
+  ptr_int = list_dof_send[imtrx];
+  for (i = total_num_send_unknowns; i > 0; i--) {
+    *ptrd++ = x[*ptr_int++];
+  }
+
+  /*
+   * store base address for the start of the external degrees of freedom
+   * in this vector
+   */
+  ptr_recv_list = x + num_internal_dofs[imtrx] + num_boundary_dofs[imtrx];
+
+  np_ptr = np_base;
+  for (p = 0; p < dpi->num_neighbors; p++) {
+    np_ptr->neighbor_ProcID = cx[p].neighbor_name;
+    np_ptr->send_message_buf = (void *)(ptr_send_list + ptr_dof_send[imtrx][p]);
+    np_ptr->send_message_length = sizeof(long long) * cx[p].num_dofs_send;
+    np_ptr->recv_message_buf = (void *)ptr_recv_list;
+    np_ptr->recv_message_length = sizeof(long long) * cx[p].num_dofs_recv;
+    ptr_recv_list += cx[p].num_dofs_recv;
+    np_ptr++;
+  }
+  exchange_neighbor_proc_info(dpi->num_neighbors, np_base);
+  safer_free((void **)&np_base);
+  safer_free((void **)&ptr_send_list);
+#endif /* PARALLEL */
+}
 /********************************************************************/
 /********************************************************************/
 /********************************************************************/