diff --git a/include/tdycore.h b/include/tdycore.h index cee1a769..d384241c 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -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*); diff --git a/src/mesh/tdycoremesh.c b/src/mesh/tdycoremesh.c index a04bcaab..7e599921 100644 --- a/src/mesh/tdycoremesh.c +++ b/src/mesh/tdycoremesh.c @@ -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 { @@ -737,19 +742,19 @@ PetscErrorCode SaveMeshConnectivityInfo(TDy tdy) { for (PetscInt s=0; sface_offset[icell]; - for (PetscInt ii=0; iinum_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; iinum_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 diff --git a/src/tdycore.c b/src/tdycore.c index 73ff6ab9..de1d9176 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -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; @@ -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); @@ -470,6 +480,9 @@ 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); @@ -477,6 +490,40 @@ PetscErrorCode TDyGetDM(TDy tdy,DM *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); @@ -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. @@ -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); @@ -577,7 +626,7 @@ 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); @@ -585,7 +634,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { 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); @@ -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, @@ -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"); } @@ -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"); } @@ -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); } @@ -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); } diff --git a/src/tdyth.c b/src/tdyth.c index 21eac4fa..09ed65ec 100644 --- a/src/tdyth.c +++ b/src/tdyth.c @@ -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); @@ -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); } diff --git a/src/tdytpf.c b/src/tdytpf.c index a2e83d2c..da0509bf 100644 --- a/src/tdytpf.c +++ b/src/tdytpf.c @@ -47,7 +47,7 @@ PetscErrorCode TDyTPFInitialize(TDy tdy) { PetscErrorCode TDyTPFComputeSystem(TDy tdy,Mat K,Vec F) { PetscErrorCode ierr; - PetscInt dim,dim2,f,fStart,fEnd,c,cStart,cEnd,row,col,junk,ss; + PetscInt dim,dim2,fStart,fEnd,c,cStart,cEnd,row,col,junk; PetscReal pnt2pnt[3],dist,Ki,p,force; DM dm = tdy->dm; PetscFunctionBegin; @@ -61,28 +61,40 @@ PetscErrorCode TDyTPFComputeSystem(TDy tdy,Mat K,Vec F) { ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); MaterialProp *matprop = tdy->matprop; - for(f=fStart; fdm, "boundary", &bnd_label); CHKERRQ(ierr); + + // Loop over all faces and compute fluxes. + for (PetscInt f = fStart; f < fEnd; ++f) { const PetscInt *supp; - ierr = DMPlexGetSupportSize(dm,f,&ss ); CHKERRQ(ierr); - ierr = DMPlexGetSupport (dm,f,&supp); CHKERRQ(ierr); + ierr = DMPlexGetSupport(dm,f,&supp); CHKERRQ(ierr); // wrt // v = -Ki (pd-p0)/dist // 00 Ki/dist // 0 -Ki pd / dist - if(ss==1 && tdy->ops->compute_boundary_pressure) { - ierr = (*tdy->ops->compute_boundary_pressure)(tdy, &(tdy->X[f*dim]),&p, tdy->boundary_pressure_ctx);CHKERRQ(ierr); - if (fabs(p+999.)<1.e-40) continue; - ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[f*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = matprop->K[(supp[0]-cStart)*dim2]; - ierr = MatSetValue(K,row,row,Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = VecSetValue(F,row,p*Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - continue; + + if (tdy->ops->compute_boundary_pressure) { + // Enforce boundary pressures if needed. Boundary faces have a label + // value of 1. + PetscInt bnd_label_val; + ierr = DMLabelGetValue(bnd_label, f, &bnd_label_val); + if (bnd_label_val == 1) { + + ierr = (*tdy->ops->compute_boundary_pressure)(tdy, &(tdy->X[f*dim]),&p, tdy->boundary_pressure_ctx);CHKERRQ(ierr); + if (fabs(p+999.)<1.e-40) continue; + ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); + Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[f*dim]),pnt2pnt); + dist = Norm(dim,pnt2pnt); + Ki = matprop->K[(supp[0]-cStart)*dim2]; + ierr = MatSetValue(K,row,row,Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); + CHKERRQ(ierr); + ierr = VecSetValue(F,row,p*Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); + CHKERRQ(ierr); + continue; + } } Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[supp[1]*dim]),pnt2pnt); @@ -147,7 +159,7 @@ PetscReal TDyTPFPressureNorm(TDy tdy,Vec U) { DM dm = tdy->dm; if(!(tdy->ops->compute_boundary_pressure)) { SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_USER, - "Must set the pressure function with TDySetDirichletValueFunction"); + "You must set the pressure function with TDySetBoundaryPressureFn"); } norm = 0; ierr = VecGetArray(U,&u); CHKERRQ(ierr); @@ -171,7 +183,7 @@ PetscReal TDyTPFPressureNorm(TDy tdy,Vec U) { PetscReal TDyTPFVelocityNorm(TDy tdy,Vec U) { PetscErrorCode ierr; - PetscInt dim,dim2,i,f,fStart,fEnd,c,cStart,cEnd,row,junk; + PetscInt dim,dim2,fStart,fEnd,cStart,cEnd,row,junk; PetscReal pnt2pnt[3],dist,Ki,p,vel[3],va,ve,*u,sign,face_error,norm,norm_sum; DM dm = tdy->dm; PetscFunctionBegin; @@ -185,21 +197,28 @@ PetscReal TDyTPFVelocityNorm(TDy tdy,Vec U) { norm = 0; MaterialProp *matprop = tdy->matprop; - for(c=cStart; cdm, "boundary", &bnd_label); CHKERRQ(ierr); + + for(PetscInt c=cStart; cops->compute_boundary_pressure) { + if((bnd_label_val == 1) && tdy->ops->compute_boundary_pressure) { ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[f*dim]),pnt2pnt); dist = Norm(dim,pnt2pnt);