Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Consolidating and simplifying command-line options #195

Merged
merged 10 commits into from
Aug 10, 2021
Merged

Consolidating and simplifying command-line options #195

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
9cd1391
Pull all options into TDySetFromOptions.
jeff-cohere Aug 9, 2021
58c3d96
Removed default_ prefix from options.
jeff-cohere Aug 9, 2021
ff62b36
Renamed TDyXYZFunctionDefault -> TDyConstantXYZFunction (material props)
jeff-cohere Aug 9, 2021
82b7210
Reorganized options into TDyOptions container and added some new ones.
jeff-cohere Aug 10, 2021
39f8873
Touched up boundary condition args.
jeff-cohere Aug 10, 2021
8a6171d
Adding mesh generation/reading options. Not finished yet.
jeff-cohere Aug 10, 2021
1de8f72
Fixed a failing test.
jeff-cohere Aug 10, 2021
8f4433a
We now use PETSc's mesh args for TDySetFromOptions.
jeff-cohere Aug 10, 2021
25c73ca
Eliminated TDyCreateDM.
jeff-cohere Aug 10, 2021
8c5c533
Updated richards and th regression tests.
jeff-cohere Aug 10, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions include/private/tdyconditionsimpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
#if !defined(TDYCONDITIONS_H)
#define TDYCONDITIONS_H

#include <petsc.h>

PETSC_INTERN PetscErrorCode TDyConstantBoundaryPressureFn(TDy,PetscReal*,PetscReal*,void*);
PETSC_INTERN PetscErrorCode TDyConstantBoundaryVelocityFn(TDy,PetscReal*,PetscReal*,void*);

#endif

21 changes: 1 addition & 20 deletions include/private/tdycoreimpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,8 @@ struct _p_TDy {
TDyTimeIntegrator ti;
TDySetupFlags setupflags;
TDyIO io;
PetscBool init_with_random_field;
PetscBool init_from_file;
char init_file[PETSC_MAX_PATH_LEN];

// options that determine the behavior(s) of the dycore
TDyOptions options;

/* arrays of the size of the Hasse diagram */
Expand All @@ -56,9 +54,6 @@ struct _p_TDy {
PetscInt ncv,nfv; /* number of {cell|face} vertices */

/* non-linear function of liquid pressure */
PetscInt rho_type;
PetscInt mu_type;
PetscInt enthalpy_type;
PetscReal *rho, *drho_dP, *d2rho_dP2; /* density of water [kg m-3]*/
PetscReal *vis, *dvis_dP, *d2vis_dP2; /* viscosity of water [Pa s] */
PetscReal *h, *dh_dP, *dh_dT; /* enthalpy of water */
Expand Down Expand Up @@ -101,12 +96,6 @@ struct _p_TDy {
void *soildensityctx;
void *soilspecificheatctx;

/* method-specific information*/
TDyMethod method;
TDyMode mode;
TDyQuadratureType qtype;
PetscBool allow_unsuitable_mesh;

/* Wheeler-Yotov */
PetscInt *vmap; /* [cell,local_vertex] --> global_vertex */
PetscInt *emap; /* [cell,local_vertex,direction] --> global_face */
Expand All @@ -122,9 +111,6 @@ struct _p_TDy {
PetscInt *faces;

/* MPFA-O */
PetscInt mpfao_gmatrix_method;
PetscInt mpfao_bc_type;

TDyMesh *mesh;
PetscReal ****subc_Gmatrix; /* Gmatrix for subcells */
PetscReal ***Trans;
Expand All @@ -149,12 +135,7 @@ struct _p_TDy {

PetscInt *closureSize, **closure, maxClosureSize;

PetscBool output_mesh;
PetscBool regression_testing;
TDyRegression *regression;

/* Timers enabled? */
PetscBool enable_timers;
};


Expand Down
10 changes: 5 additions & 5 deletions include/private/tdymaterialpropertiesimpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ PETSC_INTERN PetscBool TDyIsThermalConductivytSet(TDy);
PETSC_INTERN PetscBool TDyIsSoilSpecificHeatSet(TDy);
PETSC_INTERN PetscBool TDyIsSoilDensitySet(TDy);

PETSC_INTERN PetscErrorCode TDySoilDensityFunctionDefault(TDy,PetscReal*,PetscReal*,void*);
PETSC_INTERN PetscErrorCode TDySoilSpecificHeatFunctionDefault(TDy,PetscReal*,PetscReal*,void*);
PETSC_INTERN PetscErrorCode TDyPermeabilityFunctionDefault(TDy,PetscReal *,PetscReal *,void*);
PETSC_INTERN PetscErrorCode TDyThermalConductivityFunctionDefault(TDy,PetscReal *,PetscReal *,void*);
PETSC_INTERN PetscErrorCode TDyPorosityFunctionDefault(TDy,PetscReal *,PetscReal *,void*);
PETSC_INTERN PetscErrorCode TDyConstantSoilDensityFunction(TDy,PetscReal*,PetscReal*,void*);
PETSC_INTERN PetscErrorCode TDyConstantSoilSpecificHeatFunction(TDy,PetscReal*,PetscReal*,void*);
PETSC_INTERN PetscErrorCode TDyConstantPermeabilityFunction(TDy,PetscReal *,PetscReal *,void*);
PETSC_INTERN PetscErrorCode TDyConstantThermalConductivityFunction(TDy,PetscReal *,PetscReal *,void*);
PETSC_INTERN PetscErrorCode TDyConstantPorosityFunction(TDy,PetscReal *,PetscReal *,void*);
Comment on lines +29 to +33
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not pertinent to the name change, but we may want to make these vectorizable, taking arguments (TDy, PetscInt, const PetscReal *, PetscReal *, void *).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good idea and one we should adopt.


#endif

66 changes: 50 additions & 16 deletions include/private/tdyoptions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,22 +3,56 @@

#include <petsc.h>

typedef struct{
PetscReal gravity_constant;

// Default values for material properties
PetscReal default_porosity;
PetscReal default_permeability;
PetscReal default_soil_density;
PetscReal default_soil_specific_heat;
PetscReal default_thermal_conductivity;

// Default values for characteristic curves
PetscReal default_residual_saturation;
PetscReal default_gardner_n;
PetscReal default_vangenuchten_m;
PetscReal default_vangenutchen_alpha;
typedef struct {

// Model settings
PetscReal gravity_constant;
TDyMode mode;
TDyWaterDensityType rho_type;
PetscInt mu_type;
PetscInt enthalpy_type;

// Numerics settings
TDyMethod method;
TDyQuadratureType qtype;
PetscInt mpfao_gmatrix_method;
PetscInt mpfao_bc_type;
PetscBool tpf_allow_all_meshes;

// Constant material properties
PetscReal porosity;
PetscReal permeability;
PetscReal soil_density;
PetscReal soil_specific_heat;
PetscReal thermal_conductivity;

// Characteristic curve parameters
PetscReal residual_saturation;
PetscReal gardner_n;
PetscReal vangenuchten_m;
PetscReal vangenuchten_alpha;

// Constant boundary values
PetscReal boundary_pressure;
PetscReal boundary_temperature;
PetscReal boundary_velocity; // (normal component)

// Initial conditions
PetscBool init_with_random_field;
PetscBool init_from_file;
char init_file[PETSC_MAX_PATH_LEN];

// Mesh-related options
PetscBool generate_mesh;
PetscBool read_mesh;
char mesh_file[PETSC_MAX_PATH_LEN];

// I/O settings
PetscBool output_mesh;
PetscBool regression_testing;

// Timers enabled?
PetscBool enable_timers;
} TDyOptions;

#endif
#endif
36 changes: 19 additions & 17 deletions src/mpfao/3D/tdympfao3D_core.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ PetscErrorCode TDyComputeGMatrixMPFAOFor3DMesh(TDy tdy) {
// extract thermal conductivity tensor
PetscReal Kappa[3][3];

if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
for (PetscInt ii=0; ii<dim; ii++) {
for (PetscInt jj=0; jj<dim; jj++) {
Kappa[ii][jj] = matprop->Kappa0[icell*dim*dim + ii*dim + jj];
Expand Down Expand Up @@ -102,7 +102,7 @@ PetscErrorCode TDyComputeGMatrixMPFAOFor3DMesh(TDy tdy) {
ierr = TDyComputeEntryOfGMatrix3D(area, normal, K, nu, subcells->T[subcell_id], dim,
&(tdy->subc_Gmatrix[icell][isubcell][ii][jj])); CHKERRQ(ierr);

if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
ierr = TDySubCell_GetIthNuVector(subcells, subcell_id, jj, dim, &nu[0]); CHKERRQ(ierr);

ierr = TDyComputeEntryOfGMatrix3D(area, normal, Kappa,
Expand Down Expand Up @@ -151,7 +151,7 @@ PetscErrorCode TDyComputeGMatrixTPFFor3DMesh(TDy tdy) {
// extract thermal conductivity tensor
PetscReal Kappa[3][3];

if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
for (ii=0; ii<dim; ii++) {
for (jj=0; jj<dim; jj++) {
Kappa[ii][jj] = matprop->Kappa0[icell*dim*dim + ii*dim + jj];
Expand Down Expand Up @@ -243,7 +243,7 @@ PetscErrorCode TDyComputeGMatrixTPFFor3DMesh(TDy tdy) {
tdy->subc_Gmatrix[icell][isubcell][ii][jj] = area * (dot_prod) * K_aveg/(dist);
}

if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
if (ii == jj) {
ierr = TDySubCell_GetIthNuStarVector(subcells, subcell_id, jj, dim, &nu[0]); CHKERRQ(ierr);

Expand Down Expand Up @@ -563,7 +563,8 @@ PetscErrorCode ComputeTransmissibilityMatrix_ForNonCornerVertex(TDy tdy,

PetscInt npitf_dir_bc_all, npitf_neu_bc_all;

if (tdy->mpfao_bc_type == MPFAO_DIRICHLET_BC || tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC) {
if (tdy->options.mpfao_bc_type == MPFAO_DIRICHLET_BC ||
tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC) {
nflux_dir_bc_up = nflux_all_bc_up;
nflux_dir_bc_dn = nflux_all_bc_dn;
npitf_dir_bc_all= npitf_bc_all;
Expand Down Expand Up @@ -1020,7 +1021,7 @@ PetscErrorCode ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInte
// Need to swap entries in (*Trans)[vertices->id][:][:] such that
// Step-1. Rows correspond to faces saved in the order of vertices->face_ids[:]
// Step-2. Columns are such that boundary cells are followed by the internal cell.
// The boundary cells are in the order of vertices->cell_ids[1:3]
// The boundary cells are in the order of vertices->cell_ids[1:3]
// Note: vertices->cell_ids[0] correspond to internal cell
// Step-3. Finally move the last column as the first column. So, columns
// are in the order of vertices->cell_ids[0:3]
Expand Down Expand Up @@ -1137,7 +1138,7 @@ PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy tdy) {
PetscErrorCode ierr;

PetscFunctionBegin;
TDY_START_FUNCTION_TIMER()
TDY_START_FUNCTION_TIMER()

// If a face is shared by two cells that belong to different
// regions, zero the rows in the transmissiblity matrix
Expand All @@ -1155,7 +1156,7 @@ PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy tdy) {
PetscInt row[1];
row[0] = iface*num_subfaces + isubface;
ierr = MatZeroRows(tdy->Trans_mat,1,row,0.0,0,0); CHKERRQ(ierr);
if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
ierr = MatZeroRows(tdy->Temp_Trans_mat,1,row,0.0,0,0); CHKERRQ(ierr);
}
}
Expand All @@ -1180,22 +1181,23 @@ PetscErrorCode TDyComputeTransmissibilityMatrix3DMesh(TDy tdy) {
PetscFunctionBegin;
TDY_START_FUNCTION_TIMER()


for (ivertex=0; ivertex<mesh->num_vertices; ivertex++) {

if (!vertices->is_local[ivertex]) continue;

if (vertices->num_boundary_faces[ivertex] == 0 || vertices->num_internal_cells[ivertex] > 1) {
ierr = ComputeTransmissibilityMatrix_ForNonCornerVertex(tdy, ivertex, cells, 0); CHKERRQ(ierr);
if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
ierr = ComputeTransmissibilityMatrix_ForNonCornerVertex(tdy, ivertex, cells, 1); CHKERRQ(ierr);
}
} else {
// It is assumed that neumann boundary condition is a zero-flux boundary condition.
// Thus, compute transmissiblity entries only for dirichlet boundary condition.
if (tdy->mpfao_bc_type == MPFAO_DIRICHLET_BC || tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC) {
if (tdy->options.mpfao_bc_type == MPFAO_DIRICHLET_BC ||
tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC) {
ierr = ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInternalVertices(tdy, ivertex, cells, 0); CHKERRQ(ierr);
if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
ierr = ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInternalVertices(tdy, ivertex, cells, 1); CHKERRQ(ierr);
}
}
Expand All @@ -1205,14 +1207,14 @@ PetscErrorCode TDyComputeTransmissibilityMatrix3DMesh(TDy tdy) {
ierr = MatAssemblyBegin(tdy->Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(tdy->Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

if (tdy->mode == TH) {
if (tdy->options.mode == TH) {
ierr = MatAssemblyBegin(tdy->Temp_Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(tdy->Temp_Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}

TDyRegion *region = &mesh->region_connected;
if (region->num_cells > 0){
if (tdy->mpfao_gmatrix_method == MPFAO_GMATRIX_TPF ) {
if (tdy->options.mpfao_gmatrix_method == MPFAO_GMATRIX_TPF ) {
ierr = TDyUpdateTransmissibilityMatrix(tdy); CHKERRQ(ierr);
} else {
PetscPrintf(PETSC_COMM_WORLD,"WARNING -- Connected region option is only supported with MPFA-O TPF\n");
Expand Down Expand Up @@ -1437,11 +1439,11 @@ PetscErrorCode ComputeFacePermeabilityTensor(TDy tdy, PetscInt face_id, PetscRea
/// g = gravity vector
/// GravDis = gravity discretization term that is not dependent on the unknown variable(s)
/// such as pressure, temperautre. Thus, this term is precomputed.
///
///
/// Starnoni, M., Berre, I., Keilegavlen, E., & Nordbotten, J. M. (2019).
/// Consistent mpfa discretization for flow in the presence of gravity. Water
/// Resources Research, 55, 10105– 10118. https://doi.org/10.1029/2019WR025384
///
///
/// @param [in] tdy A TDy struct
/// @returns 0 on success, or a non-zero error code on failure
PetscErrorCode TDyComputeGravityDiscretizationFor3DMesh(TDy tdy) {
Expand Down Expand Up @@ -1492,7 +1494,7 @@ PetscErrorCode TDyComputeGravityDiscretizationFor3DMesh(TDy tdy) {

// Currently, only zero-flux neumann boundary condition is implemented.
// If the boundary condition is neumann, then gravity discretization term is zero
if (tdy->mpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up < 0 || cell_id_dn < 0)) continue;
if (tdy->options.mpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up < 0 || cell_id_dn < 0)) continue;

PetscInt cell_id;
if (cell_id_up < 0) {
Expand Down
Loading