From 2c73caba40ac4238e212b0cd85fc94ef10abe212 Mon Sep 17 00:00:00 2001 From: Weston Ortiz Date: Fri, 19 Apr 2024 09:35:22 -0600 Subject: [PATCH] lagged G term, parallel issues still --- CMakeLists.txt | 2 + include/compute_lagged_variables.h | 29 ++ include/dp_comm.h | 1 + include/linalg/sparse_matrix.h | 1 + include/mm_as.h | 2 +- include/mm_as_structs.h | 2 + include/mm_post_def.h | 9 + src/ad_momentum.cpp | 9 +- src/compute_lagged_variables.cpp | 457 +++++++++++++++++++++++++++++ src/dp_comm.c | 55 ++++ src/load_field_variables.c | 69 ++++- src/mm_as_alloc.c | 5 +- src/mm_fill_momentum.c | 14 +- src/mm_fill_stress.c | 8 +- src/mm_post_proc.c | 85 +++++- src/mm_sol_nonlinear.c | 10 + 16 files changed, 733 insertions(+), 25 deletions(-) create mode 100644 include/compute_lagged_variables.h create mode 100644 src/compute_lagged_variables.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 90d1df835..cc3431fa1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -621,6 +621,7 @@ set(GOMA_INCLUDES include/wr_exo.h include/wr_side_data.h include/wr_soln.h + include/compute_lagged_variables.h include/brkfix/bbb.h include/brkfix/exo_utils.h include/brkfix/fix.h @@ -784,6 +785,7 @@ set(GOMA_SOURCES src/wr_exo.c src/wr_side_data.c src/wr_soln.c + src/compute_lagged_variables.cpp src/brkfix/bbb.c src/brkfix/setup_fix_data.cpp src/brkfix/fix_exo_file.c) diff --git a/include/compute_lagged_variables.h b/include/compute_lagged_variables.h new file mode 100644 index 000000000..5ffb27f0a --- /dev/null +++ b/include/compute_lagged_variables.h @@ -0,0 +1,29 @@ +#ifndef GOMA_COMPUTE_LAGGED_VARIABLES_H +#define GOMA_COMPUTE_LAGGED_VARIABLES_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include "exo_struct.h" +#include "dpi.h" + +struct Lagged_Variables { + double *lagged_variables[MAX_PROB_VAR]; + int local_count; + int global_count; + int local_rows; + int *local_node_to_lagged; + double *exchange_lagged[MAX_PROB_VAR]; + int *mapping_index; + int *index; +}; + +void setup_lagged_variables(Exo_DB *exo, Dpi *dpi, struct Lagged_Variables *lv); +int compute_lagged_variables(Exo_DB * exo, Dpi *dpi, dbl * x, dbl *x_old, dbl *x_older, dbl *xdot, dbl *xdot_old, struct Lagged_Variables *lv); + +#ifdef __cplusplus +} +#endif + +#endif // GOMA_COMPUTE_LAGGED_VARIABLES_H \ No newline at end of file diff --git a/include/dp_comm.h b/include/dp_comm.h index 201dacc94..fa02b798d 100644 --- a/include/dp_comm.h +++ b/include/dp_comm.h @@ -49,4 +49,5 @@ EXTERN void exchange_node(Comm_Ex *cx, /* cx - ptr to communications exchange in Dpi *d, /* dpi - distributed processing info */ double *a); /* x - local processor node-based vector */ +EXTERN void exchange_node_int(Comm_Ex *cx, Dpi *dpi, int *x); #endif /* GOMA_DP_COMM_H */ diff --git a/include/linalg/sparse_matrix.h b/include/linalg/sparse_matrix.h index aa66558d8..807df6a3f 100644 --- a/include/linalg/sparse_matrix.h +++ b/include/linalg/sparse_matrix.h @@ -39,6 +39,7 @@ struct g_GomaSparseMatrix { GomaGlobalOrdinal n_rows; GomaGlobalOrdinal n_cols; GomaGlobalOrdinal nnz; + int *map_to_subsystem; // 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 diff --git a/include/mm_as.h b/include/mm_as.h index b1056dd9a..8c0227bfd 100644 --- a/include/mm_as.h +++ b/include/mm_as.h @@ -52,7 +52,7 @@ extern BASIS_FUNCTIONS_STRUCT **bfd; extern BASIS_FUNCTIONS_STRUCT **bfi; extern BASIS_FUNCTIONS_STRUCT **bfex; extern struct Shell_Block **shell_blocks; -extern struct Field_Variables *fv, *fv_sens; +extern struct Field_Variables *fv, *fv_sens, *fv_lagged; extern struct Diet_Field_Variables *fv_dot_dot, *fv_dot_dot_old; extern struct Diet_Field_Variables *fv_old, *fv_dot, *fv_dot_old; extern struct External_Field_Variables *efv; diff --git a/include/mm_as_structs.h b/include/mm_as_structs.h index 502917ad8..07dd352c0 100644 --- a/include/mm_as_structs.h +++ b/include/mm_as_structs.h @@ -31,6 +31,7 @@ #ifndef GOMA_MM_AS_STRUCTS_H #define GOMA_MM_AS_STRUCTS_H +#include "compute_lagged_variables.h" #include "el_elm.h" #include "mm_elem_block_structs.h" #include "mm_mp_const.h" @@ -912,6 +913,7 @@ struct Uniform_Problem_Description { turbulent_information *turbulent_info; int strong_bc_replace; dbl strong_penalty; + struct Lagged_Variables *lv; }; typedef struct Uniform_Problem_Description UPD_STRUCT; /*____________________________________________________________________________*/ diff --git a/include/mm_post_def.h b/include/mm_post_def.h index 7ed753a65..82ddad77b 100644 --- a/include/mm_post_def.h +++ b/include/mm_post_def.h @@ -416,6 +416,15 @@ typedef struct Post_Processing_Averages { } pp_Average; enum AverageExtraTypes { + AVG_LAGGED_G11, + AVG_LAGGED_G12, + AVG_LAGGED_G13, + AVG_LAGGED_G21, + AVG_LAGGED_G22, + AVG_LAGGED_G23, + AVG_LAGGED_G31, + AVG_LAGGED_G32, + AVG_LAGGED_G33, AVG_DENSITY, AVG_HEAVISIDE, AVG_VISCOSITY, diff --git a/src/ad_momentum.cpp b/src/ad_momentum.cpp index b159e7f34..3f4bfeda8 100644 --- a/src/ad_momentum.cpp +++ b/src/ad_momentum.cpp @@ -808,7 +808,8 @@ void ad_fluid_stress(ADType Pi[DIM][DIM]) { if (evss_f) { for (int a = 0; a < VIM; a++) { for (int b = 0; b < VIM; b++) { - gamma_cont[a][b] = ad_fv->G[a][b] + ad_fv->G[b][a]; + //gamma_cont[a][b] = ad_fv->G[a][b] + ad_fv->G[b][a]; + gamma_cont[a][b] = fv_lagged->G[a][b] + fv_lagged->G[b][a]; } } } else { @@ -1330,7 +1331,7 @@ int ad_calc_pspg(ADType pspg[DIM], } } else { for (p = 0; p < WIM; p++) { - div_G[p] = 0.; + div_G[p] = fv_lagged->div_G[p]; } } @@ -1924,7 +1925,7 @@ int ad_assemble_stress_sqrt_conf(dbl tt, /* parameter to vary time integration f int eqn, var; int peqn, pvar; - int evss_gradv = 0; + int evss_gradv = 1; int i, j, status, mode; ADType v[DIM]; /* Velocity field. */ @@ -2110,7 +2111,7 @@ int ad_assemble_stress_sqrt_conf(dbl tt, /* parameter to vary time integration f g[a][b] = ad_fv->grad_v[a][b]; gt[a][b] = ad_fv->grad_v[b][a]; } else { - g[a][b] = ad_fv->G[a][b]; + g[a][b] = fv_lagged->G[a][b]; gt[b][a] = g[a][b]; } } diff --git a/src/compute_lagged_variables.cpp b/src/compute_lagged_variables.cpp new file mode 100644 index 000000000..51a2f4dab --- /dev/null +++ b/src/compute_lagged_variables.cpp @@ -0,0 +1,457 @@ +#include "compute_lagged_variables.h" + +extern "C" { +#include "load_field_variables.h" +#include "dp_comm.h" +#include "dpi.h" +#include "el_geom.h" +#include "exo_struct.h" +#include "mm_as.h" +#include "mm_fill_ptrs.h" +#include "mm_fill_util.h" +#include "mm_mp.h" +#include "rf_fem.h" +#include "rf_fem_const.h" +#include "rf_node_const.h" +#include "rf_solve.h" +} + +#include +#include +#include +#include + +extern "C" void setup_lagged_variables(Exo_DB *exo, Dpi *dpi, struct Lagged_Variables *lv) { + if (upd->matrix_index[R_STRESS11] >= 0) { + // find + int local_count = 0; + int local_rows = 0; + lv->local_node_to_lagged = (int *)malloc( + (dpi->num_internal_nodes + dpi->num_boundary_nodes + dpi->num_external_nodes) * + sizeof(int)); + for (int inode = 0; + inode < (dpi->num_internal_nodes + dpi->num_boundary_nodes + dpi->num_external_nodes); + inode++) { + NODAL_VARS_STRUCT *nv = Nodes[inode]->Nodal_Vars_Info[pg->imtrx]; + /* + * Fill the vector list which points to the unknowns defined at this + * node... + */ + int inode_varType[MaxVarPerNode], inode_matID[MaxVarPerNode]; + int row_num_unknowns = fill_variable_vector(inode, inode_varType, inode_matID); + /* + * Do a check against the number of unknowns at this + * node stored in the global array + */ + if (row_num_unknowns != nv->Num_Unknowns) { + GOMA_EH(GOMA_ERROR, "Inconsistency counting unknowns."); + } + + /* + * Loop over the unknowns defined at this row node + */ + lv->local_node_to_lagged[inode] = -1; + for (int iunknown = 0; iunknown < row_num_unknowns; iunknown++) { + /* + * Retrieve the var type of the current unknown + */ + int rowVarType = inode_varType[iunknown]; + + if (rowVarType == R_STRESS11) { + lv->local_node_to_lagged[inode] = local_count; + local_count++; + if (inode < (dpi->num_internal_nodes + dpi->num_boundary_nodes)) { + local_rows++; + } + } + } + } + + int global_count; + MPI_Allreduce(&local_rows, &global_count, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + if (global_count > 0) { + lv->global_count = global_count; + lv->local_count = local_count; + lv->local_rows = local_rows; + lv->mapping_index = (int *)malloc(local_count * sizeof(int)); + int RowOffset; + MPI_Scan(&local_rows, &RowOffset, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + RowOffset -= local_rows; + for (int i = 0; i < local_rows; i++) { + lv->mapping_index[i] = RowOffset + i; + } + int *mapping_index_exchange = (int *)malloc(dpi->num_universe_nodes * sizeof(int)); + for (int i = 0; i < dpi->num_universe_nodes; i++) { + if (lv->local_node_to_lagged[i] >= 0) { + mapping_index_exchange[i] = lv->mapping_index[lv->local_node_to_lagged[i]]; + } else { + mapping_index_exchange[i] = -1; + } + } + exchange_node_int(cx[pg->imtrx], dpi, mapping_index_exchange); + for (int i = 0; i < dpi->num_universe_nodes; i++) { + if (lv->local_node_to_lagged[i] >= 0) { + lv->mapping_index[lv->local_node_to_lagged[i]] = mapping_index_exchange[i]; + } + } + free(mapping_index_exchange); + + int v_g[DIM][DIM]; + v_g[0][0] = VELOCITY_GRADIENT11; + v_g[0][1] = VELOCITY_GRADIENT12; + v_g[1][0] = VELOCITY_GRADIENT21; + v_g[1][1] = VELOCITY_GRADIENT22; + v_g[0][2] = VELOCITY_GRADIENT13; + v_g[1][2] = VELOCITY_GRADIENT23; + v_g[2][0] = VELOCITY_GRADIENT31; + v_g[2][1] = VELOCITY_GRADIENT32; + v_g[2][2] = VELOCITY_GRADIENT33; + + lv->index = (int *)malloc(MAX_VARIABLE_TYPES * sizeof(int)); + int offset = 0; + for (int i = V_FIRST; i < V_LAST; i++) { + if (i >= VELOCITY_GRADIENT11 && i <= VELOCITY_GRADIENT33) { + lv->index[i] = offset; + offset++; + } else { + lv->index[i] = -1; + } + } + + for (int i = 0; i < VIM; i++) { + for (int j = 0; j < VIM; j++) { + if (lv->index[v_g[i][j]] >= 0) { + lv->exchange_lagged[lv->index[v_g[i][j]]] = + (double *)calloc(1, dpi->num_universe_nodes * sizeof(double)); + lv->lagged_variables[lv->index[v_g[i][j]]] = + (double *)calloc(1, local_count * sizeof(double)); + } + } + } + } + } +} + +static PetscInt initialize_lagged_matrix(Mat A, + Vec bvec[MAX_PROB_VAR], + Exo_DB *exo, + Dpi *dpi, + dbl *x, + dbl *x_old, + dbl *x_older, + dbl *xdot, + dbl *xdot_old, + struct Lagged_Variables *lv) { + PetscInt *d_nnz = (PetscInt *)calloc(lv->local_rows, sizeof(PetscInt)); + PetscInt *o_nnz = (PetscInt *)calloc(lv->local_rows, sizeof(PetscInt)); + for (int eb_index = 0; eb_index < exo->num_elem_blocks; eb_index++) { + int mn = Matilda[eb_index]; + + pd = pd_glob[mn]; + cr = cr_glob[mn]; + elc = elc_glob[mn]; + elc_rs = elc_rs_glob[mn]; + gn = gn_glob[mn]; + mp = mp_glob[mn]; + vn = vn_glob[mn]; + evpl = evpl_glob[mn]; + + for (int mode = 0; mode < vn->modes; mode++) { + ve[mode] = ve_glob[mn][mode]; + } + + int e_start = exo->eb_ptr[eb_index]; + int e_end = exo->eb_ptr[eb_index + 1]; + + for (int iel = e_start; iel < e_end; iel++) { + int ielem = iel; + + int err = load_elem_dofptr(ielem, exo, x, x_old, xdot, xdot_old, 0); + GOMA_EH(err, "load_elem_dofptr"); + err = bf_mp_init(pd); + GOMA_EH(err, "bf_mp_init"); + + int eqn = R_STRESS11; + for (int i = 0; i < ei[pg->imtrx]->num_local_nodes; i++) { + int gnn_i = lv->local_node_to_lagged[Proc_Elem_Connect[ei[pg->imtrx]->iconnect_ptr + i]]; + if (gnn_i >= 0) { + int ldof_i = ei[upd->matrix_index[eqn]]->ln_to_dof[eqn][i]; + for (int j = 0; j < ei[pg->imtrx]->num_local_nodes; j++) { + int gnn_j = + lv->local_node_to_lagged[Proc_Elem_Connect[ei[pg->imtrx]->iconnect_ptr + j]]; + int ldof_j = ei[upd->matrix_index[eqn]]->ln_to_dof[eqn][j]; + if (gnn_j >= 0) { + if (ldof_i >= 0 && ldof_j >= 0 && ((gnn_i < lv->local_rows) || Num_Proc == 1)) { + if (gnn_i < lv->local_rows) { + if (gnn_j >= lv->local_rows) { + o_nnz[gnn_i] += 1; + } else { + d_nnz[gnn_i] += 1; + } + } + } + } + } + } + } + } /* END for (iel = 0; iel < num_internal_elem; iel++) */ + } /* END for (ieb loop) */ + + if (Num_Proc == 1) { + MatSeqAIJSetPreallocation(A, 0, d_nnz); + } else { + MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz); + } + + for (int eb_index = 0; eb_index < exo->num_elem_blocks; eb_index++) { + int mn = Matilda[eb_index]; + + pd = pd_glob[mn]; + cr = cr_glob[mn]; + elc = elc_glob[mn]; + elc_rs = elc_rs_glob[mn]; + gn = gn_glob[mn]; + mp = mp_glob[mn]; + vn = vn_glob[mn]; + evpl = evpl_glob[mn]; + + for (int mode = 0; mode < vn->modes; mode++) { + ve[mode] = ve_glob[mn][mode]; + } + + int e_start = exo->eb_ptr[eb_index]; + int e_end = exo->eb_ptr[eb_index + 1]; + + for (int iel = e_start; iel < e_end; iel++) { + int ielem = iel; + + int err = load_elem_dofptr(ielem, exo, x, x_old, xdot, xdot_old, 0); + GOMA_EH(err, "load_elem_dofptr"); + err = bf_mp_init(pd); + GOMA_EH(err, "bf_mp_init"); + int ielem_type = ei[pg->imtrx]->ielem_type; + int ip_total = elem_info(NQUAD, ielem_type); /* number of + * quadrature pts */ + + for (int ip = 0; ip < ip_total; ip++) { + dbl xi[3]; + dbl s, t, u; + + find_stu(ip, ielem_type, &s, &t, &u); + xi[0] = s; + xi[1] = t; + xi[2] = u; + + /* + * find quadrature weights for current ip + */ + dbl wt = Gq_weight(ip, ielem_type); + fv->wt = wt; + + /* + * Load up basis function information for ea variable... + * Old usage: fill_shape + */ + err = load_basis_functions(xi, bfd); + GOMA_EH(err, "problem from load_basis_functions"); + + err = load_fv(); + GOMA_EH(err, "load_fv"); + + /* + * This has elemental Jacobian transformation and some + * basic mesh derivatives... + * Old usage: calc_Jac, jelly_belly + */ + err = beer_belly(); + GOMA_EH(err, "beer_belly"); + + err = load_coordinate_scales(pd->CoordinateSystem, fv); + GOMA_EH(err, "load_coordinate_scales(fv)"); + + err = load_bf_grad(); + GOMA_EH(err, "load_bf_grad"); + + err = load_fv_vector(); + GOMA_EH(err, "load_fv_vector"); + + err = load_fv_grads(); + GOMA_EH(err, "load_fv_grads"); + + + int v_g[DIM][DIM]; + v_g[0][0] = VELOCITY_GRADIENT11; + v_g[0][1] = VELOCITY_GRADIENT12; + v_g[1][0] = VELOCITY_GRADIENT21; + v_g[1][1] = VELOCITY_GRADIENT22; + v_g[0][2] = VELOCITY_GRADIENT13; + v_g[1][2] = VELOCITY_GRADIENT23; + v_g[2][0] = VELOCITY_GRADIENT31; + v_g[2][1] = VELOCITY_GRADIENT32; + v_g[2][2] = VELOCITY_GRADIENT33; + int eqn = R_STRESS11; + for (int i = 0; i < ei[pg->imtrx]->num_local_nodes; i++) { + int gnn_i = lv->local_node_to_lagged[Proc_Elem_Connect[ei[pg->imtrx]->iconnect_ptr + i]]; + int ldof_i = ei[upd->matrix_index[eqn]]->ln_to_dof[eqn][i]; + if (gnn_i >= 0 && (gnn_i < lv->local_rows || Num_Proc == 1)) { + for (int a = 0; a < VIM; a++) { + for (int b = 0; b < VIM; b++) { + if (lv->index[v_g[a][b]] >= 0) { + dbl vec_contrib = + fv->grad_v[a][b] * bf[eqn]->phi[ldof_i] * fv->wt * bf[eqn]->detJ; + PetscErrorCode err = + VecSetValue(bvec[lv->index[v_g[a][b]]], lv->mapping_index[gnn_i], + vec_contrib, ADD_VALUES); + CHKERRQ(err); + } + } + } + } + for (int j = 0; j < ei[pg->imtrx]->num_local_nodes; j++) { + int gnn_j = + lv->local_node_to_lagged[Proc_Elem_Connect[ei[pg->imtrx]->iconnect_ptr + j]]; + int ldof_j = ei[upd->matrix_index[eqn]]->ln_to_dof[eqn][j]; + if (ldof_i >= 0 && ldof_j >= 0 && ((gnn_i < lv->local_rows) || Num_Proc == 1)) { + dbl mm_contrib = bf[eqn]->phi[ldof_i] * bf[eqn]->phi[ldof_j] * fv->wt * bf[eqn]->detJ; + PetscInt global_row = lv->mapping_index[gnn_i]; + PetscInt global_col = lv->mapping_index[gnn_j]; + PetscErrorCode err = MatSetValue(A, global_row, global_col, mm_contrib, ADD_VALUES); + CHKERRQ(err); + } + } + } + } /* END for (ip = 0; ip < ip_total; ip++) */ + } /* END for (iel = 0; iel < num_internal_elem; iel++) */ + } /* END for (ieb loop) */ + + free(d_nnz); + free(o_nnz); + return 0; +} + +extern "C" int compute_lagged_variables(Exo_DB *exo, + Dpi *dpi, + dbl *x, + dbl *x_old, + dbl *x_older, + dbl *xdot, + dbl *xdot_old, + struct Lagged_Variables *lv) { + PetscBool petsc_initialized = PETSC_FALSE; + PetscInitialized(&petsc_initialized); + if (!petsc_initialized) { + PetscErrorCode err = PetscInitializeNoArguments(); + CHKERRQ(err); + } + // Setup Petsc solver + PetscErrorCode ierr; + Mat A; + Vec b[MAX_PROB_VAR], sol; + KSP ksp; + PC pc; + ierr = MatCreate(PETSC_COMM_WORLD, &A); + + CHKERRQ(ierr); + ierr = MatSetOptionsPrefix(A, "lagged_proj_A"); + CHKERRQ(ierr); + ierr = MatSetFromOptions(A); + CHKERRQ(ierr); + int v_g[DIM][DIM]; + v_g[0][0] = VELOCITY_GRADIENT11; + v_g[0][1] = VELOCITY_GRADIENT12; + v_g[1][0] = VELOCITY_GRADIENT21; + v_g[1][1] = VELOCITY_GRADIENT22; + v_g[0][2] = VELOCITY_GRADIENT13; + v_g[1][2] = VELOCITY_GRADIENT23; + v_g[2][0] = VELOCITY_GRADIENT31; + v_g[2][1] = VELOCITY_GRADIENT32; + v_g[2][2] = VELOCITY_GRADIENT33; + for (int i = 0; i < VIM; i++) { + for (int j = 0; j < VIM; j++) { + // if (upd->matrix_index[v_g[i][j]] >= 0) { + std::string opt = "lagged_proj_b_" + std::to_string(v_g[i][j]); + ierr = VecCreate(PETSC_COMM_WORLD, &(b[lv->index[v_g[i][j]]])); + CHKERRQ(ierr); + ierr = + VecSetOptionsPrefix(b[lv->index[v_g[i][j]]], opt.c_str()); + CHKERRQ(ierr); + ierr = VecSetFromOptions(b[lv->index[v_g[i][j]]]); + CHKERRQ(ierr); + ierr = VecSetSizes(b[lv->index[v_g[i][j]]], lv->local_rows, + lv->global_count); + CHKERRQ(ierr); + // } + } + } + MatSetSizes(A, lv->local_rows, lv->local_rows, lv->global_count, lv->global_count); + initialize_lagged_matrix(A, b, exo, dpi, x, x_old, x_older, xdot, xdot_old, lv); + + ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + CHKERRQ(ierr); + + for (int i = 0; i < VIM; i++) { + for (int j = 0; j < VIM; j++) { + ierr = VecAssemblyBegin(b[lv->index[v_g[i][j]]]); + CHKERRQ(ierr); + ierr = VecAssemblyEnd(b[lv->index[v_g[i][j]]]); + CHKERRQ(ierr); + } + } + + ierr = VecCreate(PETSC_COMM_WORLD, &sol); + CHKERRQ(ierr); + ierr = VecSetOptionsPrefix(sol, "lagged_proj_sol"); + CHKERRQ(ierr); + ierr = VecSetFromOptions(sol); + CHKERRQ(ierr); + ierr = VecSetSizes(sol, lv->local_rows, lv->global_count); + CHKERRQ(ierr); + + ierr = KSPCreate(MPI_COMM_WORLD, &ksp); + CHKERRQ(ierr); + ierr = KSPSetOptionsPrefix(ksp, "lagged_proj_ksp"); + CHKERRQ(ierr); + ierr = KSPSetFromOptions(ksp); + CHKERRQ(ierr); + + + ierr = KSPGetPC(ksp, &pc); + CHKERRQ(ierr); + ierr = PCSetType(pc, PCBJACOBI); + CHKERRQ(ierr); + ierr = PCSetReusePreconditioner(pc, PETSC_TRUE); + CHKERRQ(ierr); + ierr = KSPSetType(ksp, KSPGMRES); + CHKERRQ(ierr); + ierr = KSPSetTolerances(ksp, 1e-12, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); + CHKERRQ(ierr); + ierr = KSPSetOperators(ksp, A, A); + + for (int i = 0; i < VIM; i++) { + for (int j = 0; j < VIM; j++) { + if (lv->index[v_g[i][j]] >= 0) { + KSPSolve(ksp, b[lv->index[v_g[i][j]]], sol); + VecDestroy(&(b[lv->index[v_g[i][j]]])); + VecGetValues(sol, lv->local_rows, lv->mapping_index, + lv->lagged_variables[lv->index[v_g[i][j]]]); + // Exchange the lagged variables + for (int n = 0; n < dpi->num_universe_nodes; n++) { + if (lv->local_node_to_lagged[n] >= 0) { + lv->exchange_lagged[lv->index[v_g[i][j]]][n] = + lv->lagged_variables[lv->index[v_g[i][j]]] + [lv->local_node_to_lagged[n]]; + } + } + exchange_node(cx[pg->imtrx], dpi, lv->exchange_lagged[lv->index[v_g[i][j]]]); + + } + } + } + MatDestroy(&A); + VecDestroy(&sol); + KSPDestroy(&ksp); + + return 0; +} \ No newline at end of file diff --git a/src/dp_comm.c b/src/dp_comm.c index 5e9199e48..5973c5423 100644 --- a/src/dp_comm.c +++ b/src/dp_comm.c @@ -196,6 +196,61 @@ void exchange_dof_int(Comm_Ex *cx, Dpi *dpi, int *x, int imtrx) /********************************************************************/ /********************************************************************/ /********************************************************************/ +void exchange_node_int(Comm_Ex *cx, Dpi *dpi, int *x) + +/************************************************************ + * + * exchange_dof_int(): + * + * send/recv appropriate pieces of a node-based int 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_node_send[dpi->num_neighbors]; + np_base = alloc_struct_1(COMM_NP_STRUCT, dpi->num_neighbors); + ptrd = (int *)alloc_int_1(total_num_send_unknowns, INT_NOINIT); + ptr_send_list = ptrd; + + /* + * gather up the list of send unknowns + */ + ptr_int = list_node_send; + for (i = total_num_send_unknowns; i > 0; i--) { + *ptrd++ = x[*ptr_int++]; + } + + /* + * store base address for the start of the entries corresponding + * to external nodes in this vector + */ + ptr_recv_list = x + dpi->num_internal_nodes + dpi->num_boundary_nodes; + + 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_node_send[p]); + np_ptr->send_message_length = sizeof(int) * cx[p].num_nodes_send; + np_ptr->recv_message_buf = (void *)ptr_recv_list; + np_ptr->recv_message_length = sizeof(int) * cx[p].num_nodes_recv; + ptr_recv_list += cx[p].num_nodes_recv; + np_ptr++; + } + exchange_neighbor_proc_info(dpi->num_neighbors, np_base); + safer_free((void **)&np_base); + safer_free((void **)&ptr_send_list); +#endif +} void exchange_node(Comm_Ex *cx, Dpi *dpi, double *x) diff --git a/src/load_field_variables.c b/src/load_field_variables.c index 7a941b7d6..0ced6bd66 100644 --- a/src/load_field_variables.c +++ b/src/load_field_variables.c @@ -192,17 +192,15 @@ int load_fv(void) status = stress_eqn_pointer(v_s); GOMA_EH(status, "stress_eqn_pointer(v_s)"); } - if (pdgv[VELOCITY_GRADIENT11]) { - v_g[0][0] = VELOCITY_GRADIENT11; - v_g[0][1] = VELOCITY_GRADIENT12; - v_g[1][0] = VELOCITY_GRADIENT21; - v_g[1][1] = VELOCITY_GRADIENT22; - v_g[0][2] = VELOCITY_GRADIENT13; - v_g[1][2] = VELOCITY_GRADIENT23; - v_g[2][0] = VELOCITY_GRADIENT31; - v_g[2][1] = VELOCITY_GRADIENT32; - v_g[2][2] = VELOCITY_GRADIENT33; - } + v_g[0][0] = VELOCITY_GRADIENT11; + v_g[0][1] = VELOCITY_GRADIENT12; + v_g[1][0] = VELOCITY_GRADIENT21; + v_g[1][1] = VELOCITY_GRADIENT22; + v_g[0][2] = VELOCITY_GRADIENT13; + v_g[1][2] = VELOCITY_GRADIENT23; + v_g[2][0] = VELOCITY_GRADIENT31; + v_g[2][1] = VELOCITY_GRADIENT32; + v_g[2][2] = VELOCITY_GRADIENT33; /* * Since it is possible to have a 1D element in a 2D problem, @@ -1265,6 +1263,20 @@ int load_fv(void) /* * Velocity Gradient (tensor)... */ + // compute lagged variables for velocity gradient + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + int v = R_STRESS11; + + fv_lagged->G[p][q] = 0.; + dofs = ei[upd->matrix_index[v]]->dof[v]; + for (i = 0; i < dofs; i++) { + int lvi = upd->lv->local_node_to_lagged[ei[upd->matrix_index[v]]->gnn_list[v][i]]; + int peqn = upd->lv->index[v_g[p][q]]; + fv_lagged->G[p][q] += upd->lv->lagged_variables[peqn][lvi] * bf[v]->phi[i]; + } + } + } for (p = 0; pdgv[VELOCITY_GRADIENT11] && p < VIM; p++) { if (gn->ConstitutiveEquation == BINGHAM_MIXED) { @@ -1297,8 +1309,12 @@ int load_fv(void) v = v_g[p][q]; if (pdgv[v]) { fv->G[p][q] = fv_old->G[p][q] = fv_dot->G[p][q] = 0.0; + fv_lagged->G[p][q] = 0.; dofs = ei[upd->matrix_index[v]]->dof[v]; for (i = 0; i < dofs; i++) { + int lvi = upd->lv->local_node_to_lagged[ei[upd->matrix_index[v]]->gnn_list[v][i]]; + int peqn = upd->ep[upd->matrix_index[v]][v]; + fv_lagged->G[p][q] += upd->lv->lagged_variables[peqn][lvi] * bf[v]->phi[i]; fv->G[p][q] += *esp->G[p][q][i] * bf[v]->phi[i]; if (pd->TimeIntegration != STEADY) { fv_old->G[p][q] += *esp_old->G[p][q][i] * bf[v]->phi[i]; @@ -2953,6 +2969,37 @@ int load_fv_grads(void) /* * grad(G) */ + int v_g[DIM][DIM]; + v_g[0][0] = VELOCITY_GRADIENT11; + v_g[0][1] = VELOCITY_GRADIENT12; + v_g[1][0] = VELOCITY_GRADIENT21; + v_g[1][1] = VELOCITY_GRADIENT22; + v_g[0][2] = VELOCITY_GRADIENT13; + v_g[1][2] = VELOCITY_GRADIENT23; + v_g[2][0] = VELOCITY_GRADIENT31; + v_g[2][1] = VELOCITY_GRADIENT32; + v_g[2][2] = VELOCITY_GRADIENT33; + for (p = 0; p < VIM; p++) { + for (q = 0; q < VIM; q++) { + for (r = 0; r < VIM; r++) { + int v = R_STRESS11; + fv_lagged->grad_G[r][p][q] = 0.0; + + for (i = 0; i < dofs; i++) { + int lvi = upd->lv->local_node_to_lagged[ei[upd->matrix_index[v]]->gnn_list[v][i]]; + int peqn = upd->lv->index[v_g[p][q]]; + fv_lagged->grad_G[r][p][q] += + upd->lv->lagged_variables[peqn][lvi] * bf[v]->grad_phi[i][r]; + } + } + } + } + for (r = 0; r < dim; r++) { + fv->div_G[r] = 0.0; + for (q = 0; q < dim; q++) { + fv_lagged->div_G[r] += fv_lagged->grad_G[q][q][r]; + } + } if (pd->gv[VELOCITY_GRADIENT11]) { v = VELOCITY_GRADIENT11; diff --git a/src/mm_as_alloc.c b/src/mm_as_alloc.c index 0ac1d6c76..094db77e2 100644 --- a/src/mm_as_alloc.c +++ b/src/mm_as_alloc.c @@ -72,7 +72,7 @@ BASIS_FUNCTIONS_STRUCT **bfd = NULL; BASIS_FUNCTIONS_STRUCT **bfi = NULL; BASIS_FUNCTIONS_STRUCT **bfex = NULL; struct Shell_Block **shell_blocks = NULL; -struct Field_Variables *fv = NULL, *fv_sens = NULL; +struct Field_Variables *fv = NULL, *fv_sens = NULL, *fv_lagged; struct Diet_Field_Variables *fv_dot_dot = NULL, *fv_dot_dot_old = NULL; struct Diet_Field_Variables *fv_old = NULL, *fv_dot = NULL; struct Diet_Field_Variables *fv_dot_old = NULL; @@ -1365,6 +1365,9 @@ int assembly_alloc(Exo_DB *exo) } memset(fv, 0, sizeof(struct Field_Variables)); + fv_lagged = alloc_struct_1(struct Field_Variables, 1); + memset(fv_lagged, 0, sizeof(struct Field_Variables)); + fv_sens = alloc_struct_1(struct Field_Variables, 1); if (Debug_Flag) { DPRINTF(stdout, "%s: Field_Variables Sensitivity @ %p has %lu bytes", yo, (void *)fv_sens, diff --git a/src/mm_fill_momentum.c b/src/mm_fill_momentum.c index 54e3c6533..68cc4942c 100644 --- a/src/mm_fill_momentum.c +++ b/src/mm_fill_momentum.c @@ -2986,9 +2986,17 @@ void fluid_stress(double Pi[DIM][DIM], STRESS_DEPENDENCE_STRUCT *d_Pi) { } if (evss_f) { - for (a = 0; a < VIM; a++) { - for (b = 0; b < VIM; b++) { - gamma_cont[a][b] = fv->G[a][b] + fv->G[b][a]; + if (pd->gv[VELOCITY_GRADIENT11]) { + for (a = 0; a < VIM; a++) { + for (b = 0; b < VIM; b++) { + gamma_cont[a][b] = fv->G[a][b] + fv->G[b][a]; + } + } + } else { + for (a = 0; a < VIM; a++) { + for (b = 0; b < VIM; b++) { + gamma_cont[a][b] = fv_lagged->G[a][b] + fv_lagged->G[b][a]; + } } } } else { diff --git a/src/mm_fill_stress.c b/src/mm_fill_stress.c index 9cb1de717..b7aab71e4 100644 --- a/src/mm_fill_stress.c +++ b/src/mm_fill_stress.c @@ -4720,7 +4720,9 @@ int assemble_gradient(dbl tt, /* parameter to vary time integration from for (a = 0; a < VIM; a++) { for (b = 0; b < VIM; b++) { - grad_v[a][b] = fv->grad_v[a][b]; + // grad_v[a][b] = fv->grad_v[a][b]; + // grad_v[a][b] = fv_lagged->G[a][b]; + grad_v[a][b] = 0; } } @@ -4800,7 +4802,7 @@ int assemble_gradient(dbl tt, /* parameter to vary time integration from /* * J_G_v */ - for (p = 0; p < WIM; p++) { + for (p = WIM; p < WIM; p++) { var = VELOCITY1 + p; if (pd->v[pg->imtrx][var]) { pvar = upd->vp[pg->imtrx][var]; @@ -4838,7 +4840,7 @@ int assemble_gradient(dbl tt, /* parameter to vary time integration from /* * J_G_d */ - for (p = 0; p < dim; p++) { + for (p = dim; p < dim; p++) { var = MESH_DISPLACEMENT1 + p; if (pd->v[pg->imtrx][var]) { pvar = upd->vp[pg->imtrx][var]; diff --git a/src/mm_post_proc.c b/src/mm_post_proc.c index 76254d902..c0b17f61a 100644 --- a/src/mm_post_proc.c +++ b/src/mm_post_proc.c @@ -3903,6 +3903,33 @@ void sum_average_nodal(double **avg_count, double **avg_sum, int global_node, do } } else { switch (pp_average[i]->type) { + case AVG_LAGGED_G11: + avg_sum[i][global_node] += fv_lagged->G[0][0]; + break; + case AVG_LAGGED_G12: + avg_sum[i][global_node] += fv_lagged->G[0][1]; + break; + case AVG_LAGGED_G13: + avg_sum[i][global_node] += fv_lagged->G[0][2]; + break; + case AVG_LAGGED_G21: + avg_sum[i][global_node] += fv_lagged->G[1][0]; + break; + case AVG_LAGGED_G22: + avg_sum[i][global_node] += fv_lagged->G[1][1]; + break; + case AVG_LAGGED_G23: + avg_sum[i][global_node] += fv_lagged->G[1][2]; + break; + case AVG_LAGGED_G31: + avg_sum[i][global_node] += fv_lagged->G[2][0]; + break; + case AVG_LAGGED_G32: + avg_sum[i][global_node] += fv_lagged->G[2][1]; + break; + case AVG_LAGGED_G33: + avg_sum[i][global_node] += fv_lagged->G[2][2]; + break; case AVG_DENSITY: { double rho = density(NULL, time); avg_sum[i][global_node] += rho; @@ -7701,7 +7728,7 @@ void rd_post_process_specs(FILE *ifp, char *input) { char filename[MAX_FNL]; char second_string[MAX_CHAR_IN_INPUT]; char ts[MAX_CHAR_IN_INPUT]; - int i, k, nargs, sz, sens_vec_ct, j; + int i, k, nargs, sz, sens_vec_ct, j=0; long save_position; int pbits[5]; @@ -9081,7 +9108,7 @@ void rd_post_process_specs(FILE *ifp, char *input) { if (!strncasecmp(variable_name, "DENSITY_AVG", strlen(variable_name))) { strcpy(pp_average[i]->type_name, "DENSITY_AVG"); pp_average[i]->non_variable_type = 1; - pp_average[i]->type = AVG_HEAVISIDE; + pp_average[i]->type = AVG_DENSITY; } else if (!strncasecmp(variable_name, "HEAVISIDE", strlen(variable_name))) { strcpy(pp_average[i]->type_name, "HEAVISIDE_AVG"); pp_average[i]->non_variable_type = 1; @@ -9106,6 +9133,60 @@ void rd_post_process_specs(FILE *ifp, char *input) { strcpy(pp_average[i]->type_name, "EM_INC_MAG"); pp_average[i]->non_variable_type = 1; pp_average[i]->type = AVG_EM_INC_MAG; + } else if (!strncasecmp(variable_name, "LAGGED_G2D", strlen(variable_name))) { + int vim = 2; + int new_items = vim*vim-1; + pp_Average **pp_average_tmp = pp_average; + sz = sizeof(pp_Average *); + pp_average = (pp_Average **)array_alloc(1, nn_average + new_items, sz); + for (int k = 0; k < nn_average; k++) { + pp_average[k] = pp_average_tmp[k]; + } + free(pp_average_tmp); + sz = sizeof(pp_Average); + for (int k = nn_average; k < (new_items + nn_average); k++) { + pp_average[k] = (pp_Average *)array_alloc(1, 1, sz); + } + nn_average += new_items; + for (int a = 0; a < vim; a++) { + for (int b = 0; b < vim; b++) { + char type_name[MAX_VAR_NAME_LNGTH]; + snprintf(type_name, MAX_VAR_NAME_LNGTH, "L_G%d%d", a+1, b+1); + strcpy(pp_average[i]->type_name, type_name); + pp_average[i]->non_variable_type = 1; + pp_average[i]->type = AVG_LAGGED_G11 + a * DIM + b; + pp_average[i]->species_index = 0; + i++; + } + } + i-=1; + } else if (!strncasecmp(variable_name, "LAGGED_G3D", strlen(variable_name))) { + int vim = 3; + int new_items = vim*vim-1; + pp_Average **pp_average_tmp = pp_average; + sz = sizeof(pp_Average *); + pp_average = (pp_Average **)array_alloc(1, nn_average + new_items, sz); + for (int k = 0; k < nn_average; k++) { + pp_average[k] = pp_average_tmp[k]; + } + free(pp_average_tmp); + sz = sizeof(pp_Average); + for (int k = nn_average; k < (new_items + nn_average); k++) { + pp_average[k] = (pp_Average *)array_alloc(1, 1, sz); + } + nn_average += new_items; + for (int a = 0; a < vim; a++) { + for (int b = 0; b < vim; b++) { + char type_name[MAX_VAR_NAME_LNGTH]; + snprintf(type_name, MAX_VAR_NAME_LNGTH, "L_G%d%d", a+1, b+1); + strcpy(pp_average[i]->type_name, type_name); + pp_average[i]->non_variable_type = 1; + pp_average[i]->type = AVG_LAGGED_G11 + a * DIM + b; + pp_average[i]->species_index = 0; + i++; + } + } + i-=1; } else if (!strncasecmp(variable_name, "EM", strlen(variable_name))) { int new_items = 6 - 1; pp_Average **pp_average_tmp = pp_average; diff --git a/src/mm_sol_nonlinear.c b/src/mm_sol_nonlinear.c index 94ee1ca50..458d4b6b5 100644 --- a/src/mm_sol_nonlinear.c +++ b/src/mm_sol_nonlinear.c @@ -86,6 +86,7 @@ #include "util/distance_helpers.h" #include "wr_exo.h" #include "wr_side_data.h" +#include "compute_lagged_variables.h" /* * EDW: The prototype for function "mf_sol_lineqn" has been moved @@ -497,6 +498,8 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, struct con_struct *con = con_ptr; struct Level_Set_Data *ls_old; + static struct Lagged_Variables *lv = NULL; + /* Had to add these for newmark-beta updates based on material */ UMI_LIST_STRUCT *curr_mat_list; NODE_INFO_STRUCT *node_ptr; @@ -572,6 +575,11 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, asdv(&res_p, numProcUnknowns); asdv(&res_m, numProcUnknowns); + if (upd->lv == NULL) { + upd->lv = (struct Lagged_Variables *)calloc(1,sizeof(struct Lagged_Variables)); + setup_lagged_variables(exo, dpi, upd->lv); + } + /* * Initialize augmenting condition arrays if needed */ @@ -777,6 +785,8 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams, } get_time(ctod); + compute_lagged_variables(exo,dpi, x, x_old, x_older, xdot, xdot_old, upd->lv); + /* * Mark the CPU time for combined assembly and solve purposes... */