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

We now construct a DMLabel for boundary faces. #199

Merged
merged 5 commits into from
Aug 14, 2021
Merged

We now construct a DMLabel for boundary faces. #199

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
129357c
Adding initialization to TDySetFromOptions.
jeff-cohere Aug 12, 2021
524e4d6
Rewiring TPF method to use our boundary label.
jeff-cohere Aug 12, 2021
5b9d6fe
Not touching MPFAO boundary face def just yet.
jeff-cohere Aug 12, 2021
dd022ff
Backed out initialization with random field (doesn't work with TH).
jeff-cohere Aug 12, 2021
16c8173
Elaborated on comments for boundary label.
jeff-cohere Aug 13, 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
2 changes: 2 additions & 0 deletions include/tdycore.h
View file
Open in desktop
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
Open in desktop
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
99 changes: 83 additions & 16 deletions src/tdycore.c
View file
Open in desktop
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,14 @@ PetscErrorCode TDyCreateGrid(TDy tdy) {
tdy->dm = dm;
}

/* compute/store plex geometry */
// Mark the grid's boundary faces.
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);
Copy link

Choose a reason for hiding this comment

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

Note that this Complete() marks the closure of every face, so we might change the comment above

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Good point. Are there cases where this isn't desirable? Just curious (I took this code from an example).

Copy link

Choose a reason for hiding this comment

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

Well, sometimes you want to do a face-specific thing, so you would need to filter out any non-faces, but that is usually easy. FE typically works with the whole boundary whereas FV just wants faces.


// Compute/store plex geometry.
PetscLogEvent t1 = TDyGetTimer("ComputePlexGeometry");
TDyStartTimer(t1);
ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr);
Expand Down Expand Up @@ -470,13 +479,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 +598,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 +616,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 +625,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 +680,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 +717,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 +730,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 +757,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 +1619,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
Open in desktop
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