Skip to content

Commit

Permalink
Merge pull request #199 from TDycores-Project/jeff-cohere/boundary-label
Browse files Browse the repository at this point in the history 
We now construct a DMLabel for boundary faces.
  • Loading branch information
@jeff-cohere
jeff-cohere authored Aug 14, 2021
2 parents 4586201 + 16c8173 commit 974f999
Show file tree
Hide file tree
Showing 5 changed files with 151 additions and 57 deletions.
2 changes: 2 additions & 0 deletions include/tdycore.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@ PETSC_EXTERN PetscErrorCode TDyView(TDy,PetscViewer viewer);

PETSC_EXTERN PetscErrorCode TDyGetDimension(TDy,PetscInt*);
PETSC_EXTERN PetscErrorCode TDyGetDM(TDy,DM*);
PETSC_EXTERN PetscErrorCode TDyGetBoundaryFaces(TDy,PetscInt*, const PetscInt**);
PETSC_EXTERN PetscErrorCode TDyRestoreBoundaryFaces(TDy,PetscInt*, const PetscInt**);
PETSC_EXTERN PetscErrorCode TDyGetCentroidArray(TDy,PetscReal**);

PETSC_EXTERN PetscErrorCode TDySetGravityVector(TDy,PetscReal*);
Expand Down
27 changes: 16 additions & 11 deletions src/mesh/tdycoremesh.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -729,6 +729,11 @@ PetscErrorCode SaveMeshConnectivityInfo(TDy tdy) {
ierr = DMPlexGetSupportSize(dm, f, &support_size); CHKERRQ(ierr);
ierr = DMPlexGetSupport(dm, f, &support); CHKERRQ(ierr);

// TODO: This is where we decide whether a face belongs to the domain
// TODO: boundary. It's logically consistent with the way that DMPlex
// TODO: decides on the domain boundary, so we can leave it like this
// TODO: for now, but we should favor the use of the "boundary" DMLabel
// TODO: in future efforts.
if (support_size == 2) {
faces->is_internal[iface] = PETSC_TRUE;
} else {
Expand All @@ -737,19 +742,19 @@ PetscErrorCode SaveMeshConnectivityInfo(TDy tdy) {

for (PetscInt s=0; s<support_size; s++) {
PetscInt icell = support[s] - c_start;
PetscBool found = PETSC_FALSE;
PetscInt cOffsetFace = cells->face_offset[icell];
for (PetscInt ii=0; ii<cells->num_faces[icell]; ii++) {
if (cells->face_ids[cOffsetFace+ii] == f-f_start) {
found = PETSC_TRUE;
break;
}
}
if (!found) {
cells->face_ids[cOffsetFace + cells->num_faces[icell]] = f-f_start;
cells->num_faces[icell]++;
PetscBool found = PETSC_FALSE;
PetscInt cOffsetFace = cells->face_offset[icell];
for (PetscInt ii=0; ii<cells->num_faces[icell]; ii++) {
if (cells->face_ids[cOffsetFace+ii] == f-f_start) {
found = PETSC_TRUE;
break;
}
}
if (!found) {
cells->face_ids[cOffsetFace + cells->num_faces[icell]] = f-f_start;
cells->num_faces[icell]++;
found = PETSC_TRUE;
}
}

// If it is a boundary face, increment the number of boundary
Expand Down
100 changes: 84 additions & 16 deletions src/tdycore.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,8 @@ PetscErrorCode TDyCreate(TDy *_tdy) {
tdy->Tref = 25;
tdy->gravity[0] = 0; tdy->gravity[1] = 0; tdy->gravity[2] = 0;
tdy->dm = NULL;
tdy->solution = NULL;
tdy->J = NULL;

/* initialize method information to null */
tdy->vmap = NULL; tdy->emap = NULL; tdy->Alocal = NULL; tdy->Flocal = NULL;
Expand Down Expand Up @@ -376,7 +378,15 @@ PetscErrorCode TDyCreateGrid(TDy tdy) {
tdy->dm = dm;
}

/* compute/store plex geometry */
// Mark the grid's boundary faces and their transitive closure. All are
// stored at their appropriate strata within the label.
DMLabel boundary_label;
ierr = DMCreateLabel(tdy->dm, "boundary"); CHKERRQ(ierr);
ierr = DMGetLabel(tdy->dm, "boundary", &boundary_label); CHKERRQ(ierr);
ierr = DMPlexMarkBoundaryFaces(tdy->dm, 1, boundary_label); CHKERRQ(ierr);
ierr = DMPlexLabelComplete(tdy->dm, boundary_label); CHKERRQ(ierr);

// Compute/store plex geometry.
PetscLogEvent t1 = TDyGetTimer("ComputePlexGeometry");
TDyStartTimer(t1);
ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr);
Expand Down Expand Up @@ -470,13 +480,50 @@ PetscErrorCode TDyGetDimension(TDy tdy,PetscInt *dim) {
PetscFunctionReturn(0);
}

/// Retrieves the DM used by the dycore. This must be called after
/// TDySetDM or TDySetFromOptions.
/// @param dm A pointer that stores the DM in use by the dycore
PetscErrorCode TDyGetDM(TDy tdy,DM *dm) {
PetscFunctionBegin;
PetscValidHeaderSpecific(tdy,TDY_CLASSID,1);
*dm = tdy->dm;
PetscFunctionReturn(0);
}

/// Retrieves the indices of the faces belonging to the domain boundary
/// for the dycore. This must be called after TDySetFromOptions. Call
/// TDyRestoreBoundaryFaces when you're finished manipulating boundary
/// faces.
/// @param num_faces A pointer that stores the number of boundary faces
/// @param faces A pointer to an array that stores the indices of boundary faces
PetscErrorCode TDyGetBoundaryFaces(TDy tdy, PetscInt *num_faces,
const PetscInt **faces) {
PetscFunctionBegin;
PetscErrorCode ierr;
DMLabel label;
IS is;
ierr = DMGetLabel(tdy->dm, "boundary", &label); CHKERRQ(ierr);
ierr = DMLabelGetStratumSize(label, 1, num_faces); CHKERRQ(ierr);
ierr = DMLabelGetStratumIS(label, 1, &is); CHKERRQ(ierr);
ierr = ISGetIndices(is, faces); CHKERRQ(ierr);
PetscFunctionReturn(0);
}

/// Resets the pointers set by TDyGetBoundaryFaces.
/// @param num_faces A pointer that stores the number of boundary faces
/// @param faces A pointer to an array that stores the indices of boundary faces
PetscErrorCode TDyRestoreBoundaryFaces(TDy tdy, PetscInt *num_faces,
const PetscInt** faces) {
PetscFunctionBegin;
PetscErrorCode ierr;
DMLabel label;
IS is;
ierr = DMGetLabel(tdy->dm, "boundary", &label); CHKERRQ(ierr);
ierr = DMLabelGetStratumIS(label, 1, &is); CHKERRQ(ierr);
ierr = ISRestoreIndices(is, faces); CHKERRQ(ierr);
PetscFunctionReturn(0);
}

PetscErrorCode TDyGetCentroidArray(TDy tdy,PetscReal **X) {
PetscFunctionBegin;
PetscValidHeaderSpecific(tdy,TDY_CLASSID,1);
Expand Down Expand Up @@ -552,8 +599,10 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
PetscErrorCode ierr;
PetscFunctionBegin;

MPI_Comm comm = PETSC_COMM_WORLD;

if ((tdy->setupflags & TDySetupFinished) != 0) {
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetupNumericalMethods()");
SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetupNumericalMethods()");
}

// Collect options from command line arguments.
Expand All @@ -568,7 +617,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
PetscBool flag;

// Material property options
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Material property options",""); CHKERRQ(ierr);
ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Material property options",""); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_porosity", "Value of porosity", NULL, options->porosity, &options->porosity, NULL); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_permability", "Value of permeability", NULL, options->permeability, &options->permeability, NULL); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_soil_density", "Value of soil density", NULL, options->soil_density, &options->soil_density, NULL); CHKERRQ(ierr);
Expand All @@ -577,15 +626,15 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
ierr = PetscOptionsEnd(); CHKERRQ(ierr);

// Characteristic curve options
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Characteristic curve options",""); CHKERRQ(ierr);
ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Characteristic curve options",""); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_residual_satuaration", "Value of residual saturation", NULL, options->residual_saturation, &options->residual_saturation, NULL); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_gardner_param_n", "Value of Gardner n parameter", NULL, options->gardner_n, &options->gardner_n, NULL); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_vangenuchten_param_m", "Value of VanGenuchten m parameter", NULL, options->vangenuchten_m, &options->vangenuchten_m, NULL); CHKERRQ(ierr);
ierr = PetscOptionsReal("-tdy_vangenuchten_param_alpha", "Value of VanGenuchten alpha parameter", NULL, options->vangenuchten_alpha, &options->vangenuchten_alpha, NULL); CHKERRQ(ierr);
ierr = PetscOptionsEnd(); CHKERRQ(ierr);

// Model options
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Model options",""); CHKERRQ(ierr);
ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Model options",""); CHKERRQ(ierr);
ierr = PetscOptionsEnum("-tdy_mode","Flow mode",
"TDySetMode",TDyModes,(PetscEnum)options->mode,
(PetscEnum *)&options->mode, NULL); CHKERRQ(ierr);
Expand Down Expand Up @@ -632,7 +681,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
ierr = PetscOptionsEnd(); CHKERRQ(ierr);

// Numerics options
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr);
ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr);
ierr = PetscOptionsEnum("-tdy_method","Discretization method",
"TDySetDiscretizationMethod",TDyMethods,
(PetscEnum)options->method,(PetscEnum *)&options->method,
Expand Down Expand Up @@ -669,7 +718,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
&options->init_from_file); CHKERRQ(ierr);

if (options->init_from_file && options->init_with_random_field) {
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,
SETERRQ(comm,PETSC_ERR_USER,
"Only one of -tdy_init_from_file and -tdy_init_with_random_field can be specified");
}

Expand All @@ -682,16 +731,16 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
sizeof(options->mesh_file),
&options->read_mesh); CHKERRQ(ierr);
if (options->generate_mesh && options->read_mesh) {
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,
SETERRQ(comm,PETSC_ERR_USER,
"Only one of -tdy_generate_mesh and -tdy_read_mesh can be specified");
}
if ((tdy->dm != NULL) && (options->generate_mesh || options->read_mesh)) {
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,
SETERRQ(comm,PETSC_ERR_USER,
"TDySetDM was called before TDySetFromOptions: can't generate or "
"read a mesh");
}
if ((tdy->dm == NULL) && !(options->generate_mesh || options->read_mesh)) {
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,
SETERRQ(comm,PETSC_ERR_USER,
"No mesh is available for TDycore: please use TDySetDM, "
"-tdy_generate_mesh, or -tdy_read_mesh to specify a mesh");
}
Expand All @@ -709,8 +758,18 @@ PetscErrorCode TDySetFromOptions(TDy tdy) {
ierr = PetscOptionsEnd(); CHKERRQ(ierr);
tdy->setupflags |= TDyOptionsSet;

// Create our mesh.
ierr = TDyCreateGrid(tdy); CHKERRQ(ierr);

// If we have been instructed to initialize the solution from a file, do so.
if (options->init_from_file) {
ierr = TDyCreateVectors(tdy); CHKERRQ(ierr);
PetscViewer viewer;
ierr = PetscViewerBinaryOpen(comm, options->init_file, FILE_MODE_READ,
&viewer); CHKERRQ(ierr);
ierr = VecLoad(tdy->solution, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}

Expand Down Expand Up @@ -1561,24 +1620,33 @@ PetscErrorCode TDyPostSolveSNESSolver(TDy tdy,Vec U) {
PetscFunctionReturn(0);
}

/// Allocates storage for the vectors used by the dycore. Storage is allocated
/// the first time the function is called. Subsequent calls have no effect.
PetscErrorCode TDyCreateVectors(TDy tdy) {
PetscErrorCode ierr;
PetscFunctionBegin;
TDY_START_FUNCTION_TIMER()
ierr = TDyCreateGlobalVector(tdy,&tdy->solution); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->residual); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->accumulation_prev); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->soln_prev); CHKERRQ(ierr);
if (tdy->solution == NULL) {
ierr = TDyCreateGlobalVector(tdy,&tdy->solution); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->residual); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->accumulation_prev); CHKERRQ(ierr);
ierr = VecDuplicate(tdy->solution,&tdy->soln_prev); CHKERRQ(ierr);
}
TDY_STOP_FUNCTION_TIMER()
PetscFunctionReturn(0);
}

/// Allocates storage for the dycore's Jacobian and preconditioner matrices.
/// Storage is allocated the first time the function is called. Subsequent calls
/// have no effect.
PetscErrorCode TDyCreateJacobian(TDy tdy) {
PetscErrorCode ierr;
PetscFunctionBegin;
TDY_START_FUNCTION_TIMER()
ierr = TDyCreateJacobianMatrix(tdy,&tdy->J); CHKERRQ(ierr);
ierr = TDyCreateJacobianMatrix(tdy,&tdy->Jpre); CHKERRQ(ierr);
if (tdy->J == NULL) {
ierr = TDyCreateJacobianMatrix(tdy,&tdy->J); CHKERRQ(ierr);
ierr = TDyCreateJacobianMatrix(tdy,&tdy->Jpre); CHKERRQ(ierr);
}
TDY_STOP_FUNCTION_TIMER()
PetscFunctionReturn(0);
}
Expand Down
10 changes: 5 additions & 5 deletions src/tdyth.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,9 @@ PetscErrorCode TDyTHTSPostStep(TS ts) {
PetscFunctionBegin;
ierr = PetscObjectGetComm((PetscObject)ts,&comm); CHKERRQ(ierr);
ierr = MPI_Comm_rank(comm,&rank); CHKERRQ(ierr);
ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
ierr = TSGetTime(ts,&time); CHKERRQ(ierr);
ierr = TSGetStepNumber(ts,&istep); CHKERRQ(ierr);
ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
ierr = TSGetTime(ts,&time); CHKERRQ(ierr);
ierr = TSGetStepNumber(ts,&istep); CHKERRQ(ierr);
ierr = TSGetSNES(ts,&snes); CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&nit); CHKERRQ(ierr); ierr = SNESGetLinearSolveIterations(snes,&lit); CHKERRQ(ierr);
ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
Expand Down Expand Up @@ -108,8 +108,8 @@ PetscErrorCode TDyTHConvergenceTest(SNES snes, PetscInt it,
ierr = VecView(r,viewer);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
#endif
if (!rank)
if (!rank)
printf("%3d: %9.2e %9.2e %9.2e\n",it,fnorm,xnorm,unorm);

PetscFunctionReturn(0);
}
Loading

0 comments on commit 974f999

Please sign in to comment.