From f2667fb8d76382c6ee32674ad670e1c34148436c Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 30 Sep 2021 14:38:04 -0700 Subject: [PATCH 01/23] Simplifying some things. Not getting rid of TDySetDM after all. --- demo/transient/transient_mpfaof90.F90 | 22 +++---------- include/private/tdyoptions.h | 1 - src/tdycore.c | 46 +++++++++++++-------------- 3 files changed, 27 insertions(+), 42 deletions(-) diff --git a/demo/transient/transient_mpfaof90.F90 b/demo/transient/transient_mpfaof90.F90 index 7adf2790..8494adaa 100644 --- a/demo/transient/transient_mpfaof90.F90 +++ b/demo/transient/transient_mpfaof90.F90 @@ -71,7 +71,7 @@ program main implicit none TDy :: tdy - DM :: dm, dmDist + DM :: dm Vec :: U TS :: ts PetscInt :: rank, successful_exit_code @@ -109,15 +109,11 @@ program main call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); CHKERRA(ierr); - call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr); + call TDySetDM(tdy, dm, ierr); + CHKERRA(ierr); + + call TDySetFromOptions(tdy,ierr); CHKERRA(ierr); - if (dmDist /= PETSC_NULL_DM) then - call DMDestroy(dm, ierr); - CHKERRA(ierr); - dm = dmDist; - end if - call DMSetUp(dm,ierr); - CHKERRA(ierr) call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr); CHKERRA(ierr); @@ -138,14 +134,6 @@ program main enddo enddo - call DMSetFromOptions(dm, ierr); - CHKERRA(ierr); - - call TDySetDM(tdy, dm, ierr); - CHKERRA(ierr); - call TDySetFromOptions(tdy,ierr); - CHKERRA(ierr); - call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr); CHKERRA(ierr); diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index f6844a0d..f96ce75f 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -43,7 +43,6 @@ typedef struct { char init_file[PETSC_MAX_PATH_LEN]; // Mesh-related options - PetscBool generate_mesh; PetscBool read_mesh; char mesh_file[PETSC_MAX_PATH_LEN]; diff --git a/src/tdycore.c b/src/tdycore.c index a0a5a1b3..01e916da 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -264,19 +264,27 @@ PetscErrorCode TDyCreate(TDy *_tdy) { PetscFunctionReturn(0); } +/// Sets the DM to be used by the dycore. Useful for generating a specific +/// geometry to be used by a driver or demo. This function can only be called +/// after TDyCreate and before TDySetFromOptions. +/// @param [in] dm a DM object for the dycore to use. PetscErrorCode TDySetDM(TDy tdy, DM dm) { PetscFunctionBegin; + if (!(tdy->setup_flags | TDyParametersInitialized)) { + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + "Can't call TDySetDM before TDyCreate!"); + } if (tdy->dm != NULL) { - if (tdy->options.generate_mesh) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM: A DM has already been generated from supplied options."); - } else { // we read a mesh from a file + if (tdy->options.read_mesh) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM: A DM has already been read from a given file."); + "Can't call TDySetDM: -tdy_read_mesh was specified!"); + } else { + // throw out the existing DM. + DMDestroy(&tdy->dm); } } if (!dm) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A DM must be created prior to TDySetDM()"); + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Cannot pass an empty DM to TDySetDM()"); } tdy->dm = dm; PetscFunctionReturn(0); @@ -373,16 +381,17 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { if (!tdy->dm) { DM dm; - if (tdy->options.generate_mesh) { + if (tdy->options.read_mesh) { + ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, + PETSC_TRUE, &dm); CHKERRQ(ierr); + } else { // Here we lean on PETSc's DM* options. ierr = DMCreate(comm, &dm); CHKERRQ(ierr); ierr = DMSetType(dm, DMPLEX); CHKERRQ(ierr); ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); - } else { // if (tdy->options.read_mesh) - ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, - PETSC_TRUE, &dm); CHKERRQ(ierr); } + // Distribute the mesh, however we got it. + ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); tdy->dm = dm; } @@ -714,21 +723,10 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Mesh-related options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); - ierr = PetscOptionsBool("-tdy_generate_mesh","Generate a mesh using provided PETSc DM options","",options->generate_mesh,&(options->generate_mesh),NULL); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_mesh", options->mesh_file,sizeof(options->mesh_file),&options->read_mesh); CHKERRQ(ierr); - if (options->generate_mesh && options->read_mesh) { - 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(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)) { + if ((tdy->dm != NULL) && options->read_mesh) { 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"); + "TDySetDM was called before TDySetFromOptions: can't read a mesh"); } ierr = PetscOptionsBool("-tdy_output_mesh","Enable output of mesh attributes","",options->output_mesh,&(options->output_mesh),NULL); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_output_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->output_geom_attributes); CHKERRQ(ierr); From 7dd3814c4a96752151e84fb5ce7391d1334ad70c Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Wed, 13 Oct 2021 17:00:16 -0700 Subject: [PATCH 02/23] Interim commit -- in the middle of wrapping TDyDefineDefaultDM. --- include/private/tdycoreimpl.h | 9 ++++---- include/tdycore.h | 2 +- src/f90-mod/tdycoremod.F90 | 19 +++++++++++++++++ src/interface/ftn/tdycoref.c | 12 ----------- src/tdycore.c | 39 +++++++++++++---------------------- 5 files changed, 39 insertions(+), 42 deletions(-) diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index 89f08d74..a882c65c 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -28,15 +28,16 @@ struct _TDyOps { // viewer. PetscErrorCode (*view)(void*, PetscViewer); - //PetscErrorCode (*setup)(TDy); + // Configures a default DM in the absence of grid information. + PetscErrorCode (*config_default_dm)(void*, DM); + + // Called by TDySetup -- configures the DM for the dycore. + PetscErrorCode (*setup)(TDy, DM); // Called by TDySetFromOptions -- sets implementation-specific options // from command-line arguments. PetscErrorCode (*set_from_options)(void*); - // Called by TDySetup -- configures the DM for the dycore. - PetscErrorCode (*config_dm)(void*, DM); - // Called by TDyComputeErrorNorms -- computes error norms given a solution // vector. PetscErrorCode (*compute_error_norms)(void*,Vec,PetscReal*,PetscReal*); diff --git a/include/tdycore.h b/include/tdycore.h index 6d909c0b..7de55b13 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -155,7 +155,7 @@ PETSC_EXTERN PetscErrorCode TDySetMPFAOGmatrixMethod(TDy,TDyMPFAOGmatrixMethod); PETSC_EXTERN PetscErrorCode TDySetMPFAOBoundaryConditionType(TDy,TDyMPFAOBoundaryConditionType); // We will probably remove the following functions. -PETSC_EXTERN PetscErrorCode TDySetDM(TDy,DM); +PETSC_EXTERN PetscErrorCode TDyDefineDefaultDM(TDy,PetscErrorCode(*)(TDy, DM)); PETSC_EXTERN PetscErrorCode TDyComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscErrorCode TDySetIFunction(TS,TDy); PETSC_EXTERN PetscErrorCode TDySetIJacobian(TS,TDy); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index f217c9da..e0bc74e2 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -461,6 +461,25 @@ function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) end function end interface + abstract interface + subroutine TDyDefaultDMFunction(tdy, dm, ierr) + use tdycoredef + TDy :: tdy + DM :: dm + PetscErrorCode :: ierr + end subroutine + end interface + + interface + function DefineDefaultDM(name, c_func) bind (c, name="TDyDefineDefaultDM") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + TDy :: tdy + type(c_funptr) :: c_func + integer(c_int) :: ierr + end function + end interface + contains subroutine TDyInit(ierr) diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 4058835b..c2c8d107 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -10,7 +10,6 @@ #define tdyinitnoarguments_ TDYINITNOARGUMENTS #define tdyfinalize_ TDYFINALIZE #define tdycreate_ TDYCREATE -#define tdycreatesetdm_ TDYSETDM #define tdysetmode_ TDYSETMODE #define tdysetdiscretizationmethod_ TDYSETDISCRETIZATIONMETHOD #define tdysetfromoptions_ TDYSETFROMOPTIONS @@ -72,7 +71,6 @@ #define tdyinitnoarguments_ tdyinitnoarguments #define tdyfinalize_ tdyfinalize #define tdycreate_ tdycreate -#define tdysetdm_ tdysetdm #define tdysetmode_ tdysetmode #define tdysetdiscretizationmethod_ tdysetdiscretizationmethod #define tdysetfromoptions_ tdysetfromoptions @@ -192,16 +190,6 @@ PETSC_EXTERN void tdysetdiscretizationmethod_(TDy _tdy, PetscInt *method, int * } #endif -#if defined(__cplusplus) -extern "C" { -#endif -PETSC_EXTERN void tdysetdm_(TDy _tdy, DM dm, int *__ierr){ -*__ierr = TDySetDM((TDy)PetscToPointer((_tdy)), (DM)PetscToPointer((dm))); -} -#if defined(__cplusplus) -} -#endif - #if defined(__cplusplus) extern "C" { #endif diff --git a/src/tdycore.c b/src/tdycore.c index 01e916da..5f9ee994 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -264,29 +264,11 @@ PetscErrorCode TDyCreate(TDy *_tdy) { PetscFunctionReturn(0); } -/// Sets the DM to be used by the dycore. Useful for generating a specific -/// geometry to be used by a driver or demo. This function can only be called -/// after TDyCreate and before TDySetFromOptions. -/// @param [in] dm a DM object for the dycore to use. -PetscErrorCode TDySetDM(TDy tdy, DM dm) { - PetscFunctionBegin; - if (!(tdy->setup_flags | TDyParametersInitialized)) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM before TDyCreate!"); - } - if (tdy->dm != NULL) { - if (tdy->options.read_mesh) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM: -tdy_read_mesh was specified!"); - } else { - // throw out the existing DM. - DMDestroy(&tdy->dm); - } - } - if (!dm) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Cannot pass an empty DM to TDySetDM()"); - } - tdy->dm = dm; +/// Provides a function to be used to create a default DM in the case that no +/// grid-related command line arguments are provided to the dycore. +/// @param [in] dm_func A function that creates a DM and returns an error code +PetscErrorCode TDyDefineDefaultDM(TDy tdy, PetscErrorCode (*dm_func)(TDy, DM)) { + tdy->ops->config_default_dm = dm_func; PetscFunctionReturn(0); } @@ -385,9 +367,12 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, PETSC_TRUE, &dm); CHKERRQ(ierr); } else { + ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); + if (tdy->ops->config_default_dm) { + // We've been instructed to configure a default DM. + ierr = tdy->ops->config_default_dm(tdy->context, dm); + } // Here we lean on PETSc's DM* options. - ierr = DMCreate(comm, &dm); CHKERRQ(ierr); - ierr = DMSetType(dm, DMPLEX); CHKERRQ(ierr); ierr = DMSetFromOptions(dm); CHKERRQ(ierr); } // Distribute the mesh, however we got it. @@ -801,6 +786,10 @@ PetscErrorCode TDySetup(TDy tdy) { } TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); + + // Perform implementation-specific setup. + ierr = tdy->ops->setup(tdy->context, tdy->dm); CHKERRQ(ierr); + // TODO: Stick a call to tdy->ops->config_dm here to configure the DM for our // TODO: specific dycore setup. ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); From fe47ecc5881e3ff1176062378b783a5717bea1e0 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 14 Oct 2021 14:58:09 -0700 Subject: [PATCH 03/23] Got rid of TDySetDM in favor of TDySetDMConstructor. --- demo/steady/steady.c | 165 ++++++++++++++++++++-------------- include/private/tdycoreimpl.h | 10 +-- include/tdycore.h | 8 +- src/f90-mod/tdycoremod.F90 | 50 ++++++----- src/interface/ftn/tdycoref.c | 2 +- src/tdycore.c | 82 +++++++++++------ 6 files changed, 192 insertions(+), 125 deletions(-) diff --git a/demo/steady/steady.c b/demo/steady/steady.c index f67a4bea..ac3d7914 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -264,7 +264,6 @@ PetscErrorCode Forcing3(TDy tdy,PetscReal *x,PetscReal *f,void *ctx) { PetscFunctionReturn(0); } - PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { PetscErrorCode ierr; PetscFunctionBegin; @@ -275,7 +274,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { PetscInt v,vStart,vEnd,offset,value,dim; ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); /* this is the 'marker' label which marks boundary entities */ - ierr = DMGetLabelByNum(dm,2,&label); CHKERRQ(ierr); + ierr = DMGetLabel(dm, "boundary", &label); CHKERRQ(ierr); ierr = DMGetCoordinateSection(dm, &coordSection); CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); CHKERRQ(ierr); @@ -303,6 +302,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { ierr = VecRestoreArray(coordinates,&coords); CHKERRQ(ierr); PetscFunctionReturn(0); } + PetscErrorCode PerturbVerticesSmooth(DM dm) { PetscErrorCode ierr; PetscFunctionBegin; @@ -393,7 +393,9 @@ PetscErrorCode SaveVertices(DM dm, char filename[256]){ Vec coordinates; PetscViewer viewer; - ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); ierr = VecView(coordinates,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); @@ -414,7 +416,9 @@ PetscErrorCode SaveCentroids(DM dm, char filename[256]){ ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); - ierr = VecCreateMPI(PETSC_COMM_WORLD,(cEnd-cStart)*dim,PETSC_DETERMINE,¢roids); CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + ierr = VecCreateMPI(comm,(cEnd-cStart)*dim,PETSC_DETERMINE,¢roids); CHKERRQ(ierr); ierr = VecGetArray(centroids,¢roids_p); CHKERRQ(ierr); for(c=cStart; csetup_flags |= TDyCreated; @@ -264,11 +264,25 @@ PetscErrorCode TDyCreate(TDy *_tdy) { PetscFunctionReturn(0); } -/// Provides a function to be used to create a default DM in the case that no -/// grid-related command line arguments are provided to the dycore. -/// @param [in] dm_func A function that creates a DM and returns an error code -PetscErrorCode TDyDefineDefaultDM(TDy tdy, PetscErrorCode (*dm_func)(TDy, DM)) { - tdy->ops->config_default_dm = dm_func; +/// Provides a function to be used to create a DM in special cases where a +/// specific geometry is needed. After the function is executed on the DM, +/// DMSetFromOptions is called to apply overrides, and then the DM is +/// distributed appropriately. This function must be called before +/// TDySetFromOptions. +/// @param [in] dm_func A function that creates a given DM and returns an error +PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(void*, MPI_Comm, DM*)) { + PetscFunctionBegin; + MPI_Comm comm; + int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if (!tdy->context) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetMode and TDySetDiscretization before TDySetDMConstructor()"); + } + if (tdy->setup_flags & TDyOptionsSet) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDyDefineDM before TDySetFromOptions()"); + } + tdy->ops->create_dm = dm_func; PetscFunctionReturn(0); } @@ -276,8 +290,11 @@ PetscErrorCode TDyMalloc(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if (!tdy->dm) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"tdy->dm must be set prior to TDyMalloc()"); + SETERRQ(comm,PETSC_ERR_USER,"tdy->dm must be set prior to TDyMalloc()"); } PetscInt dim; @@ -356,7 +373,10 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { PetscScalar *coords; PetscErrorCode ierr; PetscFunctionBegin; - MPI_Comm comm = PETSC_COMM_WORLD; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if ((tdy->setup_flags & TDyOptionsSet) == 0) { SETERRQ(comm,PETSC_ERR_USER,"Options must be set prior to TDyCreateGrid()"); } @@ -367,12 +387,13 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, PETSC_TRUE, &dm); CHKERRQ(ierr); } else { - ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); - if (tdy->ops->config_default_dm) { - // We've been instructed to configure a default DM. - ierr = tdy->ops->config_default_dm(tdy->context, dm); + if (tdy->ops->create_dm) { + // We've been instructed to create a DM ourselves. + ierr = tdy->ops->create_dm(tdy->context, comm, &dm); + } else { + ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); } - // Here we lean on PETSc's DM* options. + // Here we lean on PETSc's DM* options for overrides. ierr = DMSetFromOptions(dm); CHKERRQ(ierr); } // Distribute the mesh, however we got it. @@ -545,6 +566,9 @@ PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { PetscFunctionBegin; PetscValidPointer(tdy,1); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); @@ -566,7 +590,7 @@ PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { case 3: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretizationMethod"); + SETERRQ(comm,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretizationMethod"); break; } // if (tdy->subc_Gmatrix) { ierr = TDyDeallocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, @@ -608,7 +632,8 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; - MPI_Comm comm = PETSC_COMM_WORLD; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); if ((tdy->setup_flags & TDySetupFinished) != 0) { SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); @@ -707,7 +732,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { } // Mesh-related options - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); + ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_mesh", options->mesh_file,sizeof(options->mesh_file),&options->read_mesh); CHKERRQ(ierr); if ((tdy->dm != NULL) && options->read_mesh) { SETERRQ(comm,PETSC_ERR_USER, @@ -717,7 +742,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsGetString(NULL,NULL,"-tdy_output_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->output_geom_attributes); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->read_geom_attributes); CHKERRQ(ierr); if (options->output_geom_attributes && options->read_geom_attributes){ - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Only one of -tdy_output_geom_attributes and -tdy_read_geom_attributes can be specified"); + SETERRQ(comm,PETSC_ERR_USER,"Only one of -tdy_output_geom_attributes and -tdy_read_geom_attributes can be specified"); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); @@ -777,12 +802,16 @@ PetscErrorCode TDySetupDiscretizationScheme(TDy tdy) { } PetscErrorCode TDySetup(TDy tdy) { - /* must follow TDySetFromOptions() is it relies upon options set by - TDySetFromOptions */ PetscErrorCode ierr; PetscFunctionBegin; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + /* must follow TDySetFromOptions() is it relies upon options set by + TDySetFromOptions */ if ((tdy->setup_flags & TDyOptionsSet) == 0) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); + SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); } TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); @@ -790,8 +819,6 @@ PetscErrorCode TDySetup(TDy tdy) { // Perform implementation-specific setup. ierr = tdy->ops->setup(tdy->context, tdy->dm); CHKERRQ(ierr); - // TODO: Stick a call to tdy->ops->config_dm here to configure the DM for our - // TODO: specific dycore setup. ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); if (tdy->options.regression_testing) { /* must come after Sections are set up in @@ -800,7 +827,7 @@ PetscErrorCode TDySetup(TDy tdy) { } if (tdy->options.output_mesh) { if (tdy->options.method != MPFA_O) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + SETERRQ(comm,PETSC_ERR_USER, "-tdy_output_mesh only supported for MPFA-O method"); } } @@ -1068,6 +1095,9 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&P);CHKERRQ(ierr); ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&temp);CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if (tdy->options.mode == TH) { for (PetscInt c=0;cvg_m[c],alpha,cc->sr[i],tdy->Pref-P[i],&(cc->S[i]),&cc->dS_dP[i],&(cc->d2S_dP2[i])); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unknown saturation function"); + SETERRQ(comm,PETSC_ERR_SUP,"Unknown saturation function"); break; } @@ -1110,7 +1140,7 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { RelativePermeability_Mualem(cc->mualem_m[c],cc->mualem_poly_low[c],cc->mualem_poly_coeffs[c],Se,&Kr,&dKr_dSe); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unknown relative permeability function"); + SETERRQ(comm,PETSC_ERR_SUP,"Unknown relative permeability function"); break; } cc->Kr[i] = Kr; From cfa5fa21e8808238963e4915ed574e69a05546bb Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 14 Oct 2021 15:19:30 -0700 Subject: [PATCH 04/23] Changed TDySetDiscretizationMethod -> TDySetDiscretization. --- include/tdycore.h | 5 ++--- src/interface/ftn/tdycoref.c | 8 ++++---- src/tdycore.c | 6 +++--- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/include/tdycore.h b/include/tdycore.h index 14c5510b..9155983c 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -81,7 +81,7 @@ PETSC_EXTERN PetscErrorCode TDyFinalize(void); PETSC_EXTERN PetscErrorCode TDyCreate(MPI_Comm, TDy*); PETSC_EXTERN PetscErrorCode TDySetMode(TDy,TDyMode); -PETSC_EXTERN PetscErrorCode TDySetDiscretizationMethod(TDy,TDyMethod); +PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyMethod); PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,PetscErrorCode(*)(void*, MPI_Comm, DM*)); PETSC_EXTERN PetscErrorCode TDySetFromOptions(TDy); PETSC_EXTERN PetscErrorCode TDySetup(TDy); @@ -147,8 +147,7 @@ PETSC_EXTERN PetscErrorCode TDyGetNumCellsLocal(TDy,PetscInt*); PETSC_EXTERN PetscErrorCode TDyGetCellNaturalIDsLocal(TDy,PetscInt*,PetscInt[]); PETSC_EXTERN PetscErrorCode TDyGetCellIsLocal(TDy,PetscInt*,PetscInt[]); - -PETSC_EXTERN PetscErrorCode TDyResetDiscretizationMethod(TDy); +PETSC_EXTERN PetscErrorCode TDyResetDiscretization(TDy); PETSC_EXTERN PetscErrorCode TDySetQuadratureType(TDy,TDyQuadratureType); PETSC_EXTERN PetscErrorCode TDySetWaterDensityType(TDy,TDyWaterDensityType); diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 9e322360..612a23a2 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -11,7 +11,7 @@ #define tdyfinalize_ TDYFINALIZE #define tdycreate_ TDYCREATE #define tdysetmode_ TDYSETMODE -#define tdysetdiscretizationmethod_ TDYSETDISCRETIZATIONMETHOD +#define tdysetdiscretization_ TDYSETDISCRETIZATION #define tdysetfromoptions_ TDYSETFROMOPTIONS #define tdydriverinitializettdy_ TDYDRIVERINITIALIZETDY #define tdydtimeintegratorruntotime_ TDYTIMEINTEGRATORRUNTOTIME @@ -72,7 +72,7 @@ #define tdyfinalize_ tdyfinalize #define tdycreate_ tdycreate #define tdysetmode_ tdysetmode -#define tdysetdiscretizationmethod_ tdysetdiscretizationmethod +#define tdysetdiscretization_ tdysetdiscretization #define tdysetfromoptions_ tdysetfromoptions #define tdydriverinitializettdy_ tdydriverinitializetdy #define tdydtimeintegratorruntotime_ tdydtimeintegratorruntotime @@ -183,8 +183,8 @@ PETSC_EXTERN void tdysetmode_(TDy _tdy, PetscInt *mode, int *__ierr){ #if defined(__cplusplus) extern "C" { #endif -PETSC_EXTERN void tdysetdiscretizationmethod_(TDy _tdy, PetscInt *method, int *__ierr){ -*__ierr = TDySetDiscretizationMethod((TDy)PetscToPointer((_tdy)), *method); +PETSC_EXTERN void tdysetdiscretization_(TDy _tdy, PetscInt *method, int *__ierr){ +*__ierr = TDySetDiscretization((TDy)PetscToPointer((_tdy)), *method); } #if defined(__cplusplus) } diff --git a/src/tdycore.c b/src/tdycore.c index 10f1c863..6e53c7ce 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -471,7 +471,7 @@ PetscErrorCode TDyDestroy(TDy *_tdy) { if (!tdy) PetscFunctionReturn(0); - ierr = TDyResetDiscretizationMethod(tdy); CHKERRQ(ierr); + ierr = TDyResetDiscretization(tdy); CHKERRQ(ierr); ierr = PetscFree(tdy->V); CHKERRQ(ierr); ierr = PetscFree(tdy->X); CHKERRQ(ierr); ierr = PetscFree(tdy->N); CHKERRQ(ierr); @@ -560,7 +560,7 @@ PetscErrorCode TDyGetCentroidArray(TDy tdy,PetscReal **X) { PetscFunctionReturn(0); } -PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { +PetscErrorCode TDyResetDiscretization(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; @@ -590,7 +590,7 @@ PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { case 3: break; default: - SETERRQ(comm,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretizationMethod"); + SETERRQ(comm,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretization"); break; } // if (tdy->subc_Gmatrix) { ierr = TDyDeallocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, From dff8509b1b54cbd77540d06e3132f11d2ab5d183 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 14 Oct 2021 15:26:43 -0700 Subject: [PATCH 05/23] Missed a few spots for TDySetDiscretization rename. --- demo/steady/steady.c | 10 ++++++---- src/tdycore.c | 4 ++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/demo/steady/steady.c b/demo/steady/steady.c index ac3d7914..4e8274c6 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -555,7 +555,7 @@ PetscErrorCode GeometryColumn(DM dm){ PetscFunctionReturn(0); } -// This set of globals are used by CreateDM below to create a DM for this demo. +// This set of globals is used by CreateDM below to create a DM for this demo. typedef struct DMOptions { PetscInt dim; // Dimension of DM (2 or 3) PetscInt N; // Number of cells on a side @@ -611,7 +611,7 @@ int main(int argc, char **argv) { MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy; ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-dim","Problem dimension","",dim,&dim,NULL); CHKERRQ(ierr); @@ -638,11 +638,13 @@ int main(int argc, char **argv) { dm_options_.exo = exo; if (exo) strcpy(dm_options_.exofile, exofile); - // Set up a default DM for this demo. + // Specify a special DM to be constructed for this demo. ierr = TDySetDMConstructor(tdy, CreateDM); CHKERRQ(ierr); + + // Apply overrides. ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); - /* Setup problem parameters */ + // Setup problem parameters if(wheeler2006){ if(dim != 2){ SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2006 is only for -dim 2 problems"); diff --git a/src/tdycore.c b/src/tdycore.c index 6e53c7ce..1878d45d 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -709,7 +709,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Numerics options ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_method","Discretization method", - "TDySetDiscretizationMethod",TDyMethods, + "TDySetDiscretization",TDyMethods, (PetscEnum)options->method,(PetscEnum *)&options->method, NULL); CHKERRQ(ierr); TDyQuadratureType qtype = FULL; @@ -837,7 +837,7 @@ PetscErrorCode TDySetup(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetDiscretizationMethod(TDy tdy,TDyMethod method) { +PetscErrorCode TDySetDiscretization(TDy tdy,TDyMethod method) { PetscFunctionBegin; tdy->options.method = method; PetscFunctionReturn(0); From 755633a038ea364fbd941e6d06f7b166cf4dd79a Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 14 Oct 2021 16:14:47 -0700 Subject: [PATCH 06/23] Made way for polymorphic setup function. --- include/tdycore.h | 6 ++-- src/tdycore.c | 90 +++++++++++++++++++++++++---------------------- 2 files changed, 52 insertions(+), 44 deletions(-) diff --git a/include/tdycore.h b/include/tdycore.h index 9155983c..d83f804a 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -12,7 +12,8 @@ typedef enum { MPFA_O_DAE, /* multipoint flux approximation - O method using DAE */ MPFA_O_TRANSIENTVAR, /* multipoint flux approximation - O method using TS transient (conservative) approach */ BDM, /* P0,BDM1 spaces, standard approach */ - WY /* P0,BDM1 spaces, vertex quadrature, statically condensed */ + WY, /* P0,BDM1 spaces, vertex quadrature, statically condensed */ + UNSPECIFIED_METHOD } TDyMethod; PETSC_EXTERN const char *const TDyMethods[]; @@ -41,7 +42,8 @@ PETSC_EXTERN const char *const TDyQuadratureTypes[]; typedef enum { RICHARDS=0, - TH + TH, + UNSPECIFIED_MODE } TDyMode; typedef enum { diff --git a/src/tdycore.c b/src/tdycore.c index 1878d45d..390424d8 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -27,6 +27,7 @@ const char *const TDyMethods[] = { "MPFA_O_TRANSIENTVAR", "BDM", "WY", + "UNSPECIFIED_METHOD", /* */ "TDyMethod","TDY_METHOD_",NULL }; @@ -49,6 +50,7 @@ const char *const TDyMPFAOBoundaryConditionTypes[] = { const char *const TDyModes[] = { "RICHARDS", "TH", + "UNSPECIFIED_MODE", /* */ "TDyMode","TDY_MODE_",NULL }; @@ -194,8 +196,8 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { PetscFunctionBegin; TDyOptions *options = &tdy->options; - options->mode = RICHARDS; - options->method = WY; + options->mode = UNSPECIFIED_MODE; + options->method = UNSPECIFIED_METHOD; options->gravity_constant = 9.8068; options->rho_type = WATER_DENSITY_CONSTANT; options->mu_type = WATER_VISCOSITY_CONSTANT; @@ -274,7 +276,8 @@ PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(void*, MPI PetscFunctionBegin; MPI_Comm comm; int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - if (!tdy->context) { + if ((tdy->options.mode == UNSPECIFIED_MODE) || + (tdy->options.method == UNSPECIFIED_METHOD)) { SETERRQ(comm,PETSC_ERR_USER, "You must call TDySetMode and TDySetDiscretization before TDySetDMConstructor()"); } @@ -636,7 +639,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); if ((tdy->setup_flags & TDySetupFinished) != 0) { - SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); + SETERRQ(comm,PETSC_ERR_USER,"You must call TDySetFromOptions before TDySetup()"); } // Collect options from command line arguments. @@ -768,39 +771,6 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetupDiscretizationScheme(TDy tdy) { - MPI_Comm comm; - PetscErrorCode ierr; - PetscValidPointer(tdy,1); - PetscFunctionBegin; - ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - ierr = TDyTPFInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O_DAE: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O_TRANSIENTVAR: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case BDM: - ierr = TDyBDMInitialize(tdy); CHKERRQ(ierr); - break; - case WY: - ierr = TDyWYInitialize(tdy); CHKERRQ(ierr); - break; - } - - // Record metadata for scaling studies. - TDySetTimingMetadata(tdy); - - PetscFunctionReturn(0); -} - PetscErrorCode TDySetup(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; @@ -808,10 +778,14 @@ PetscErrorCode TDySetup(TDy tdy) { MPI_Comm comm; ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - /* must follow TDySetFromOptions() is it relies upon options set by - TDySetFromOptions */ + if ((tdy->options.mode == UNSPECIFIED_MODE) || + (tdy->options.method == UNSPECIFIED_METHOD)) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetMode and TDySetDiscretization before TDySetup()"); + } + if ((tdy->setup_flags & TDyOptionsSet) == 0) { - SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); + SETERRQ(comm,PETSC_ERR_USER,"You must call TDySetFromOptions before TDySetup()"); } TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); @@ -819,10 +793,12 @@ PetscErrorCode TDySetup(TDy tdy) { // Perform implementation-specific setup. ierr = tdy->ops->setup(tdy->context, tdy->dm); CHKERRQ(ierr); - ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); + // Record metadata for scaling studies. + TDySetTimingMetadata(tdy); + if (tdy->options.regression_testing) { /* must come after Sections are set up in - TDySetupDiscretizationScheme->XXXInitialize */ + TDySetupDiscretization->XXXInitialize */ ierr = TDyRegressionInitialize(tdy); CHKERRQ(ierr); } if (tdy->options.output_mesh) { @@ -839,6 +815,36 @@ PetscErrorCode TDySetup(TDy tdy) { PetscErrorCode TDySetDiscretization(TDy tdy,TDyMethod method) { PetscFunctionBegin; + + MPI_Comm comm; + PetscErrorCode ierr; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + if (tdy->options.mode == UNSPECIFIED_MODE) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetMode before TDySetDiscretization"); + } + + // FIXME: All of these functions are current called TDy*Initialize, and + // FIXME: should be renamed and altered to accept a context pointer and a DM + // FIXME: as arguments. + if (tdy->options.mode == RICHARDS) { + if (method == TPF) { + tdy->ops->setup = TDyTPFRichardsSetup; + } else if (method == MPFA_O) { + tdy->ops->setup = TDyMPFAORichardsSetup; + } else if (method == MPFA_O_DAE) { + tdy->ops->setup = TDyMPFAORichardsSetup; + } else if (method == MPFA_O_TRANSIENTVAR) { + tdy->ops->setup = TDyMPFAORichardsSetup; + } else if (method == BDM) { + tdy->ops->setup = TDyBDMRichardsSetup; + } else if (method == WY) { + tdy->ops->setup = TDyWYRichardsSetup; + } + } else if (tdy->options.mode == TH) { + // How many methods do we support for TH? + } tdy->options.method = method; PetscFunctionReturn(0); } From c59b60df3362b330f28c0cb6f71eaec7a525ccec Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 08:42:05 -0700 Subject: [PATCH 07/23] Renaming "method" to "discretization" for clarity. --- include/private/tdycoreimpl.h | 2 + include/private/tdyoptions.h | 2 +- include/tdycore.h | 53 +++++++++++-------- src/mpfao/tdympfao.c | 2 +- src/tdycore.c | 95 +++++++++++++++++++---------------- 5 files changed, 89 insertions(+), 65 deletions(-) diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index 0d8d7585..68617afb 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -42,6 +42,8 @@ struct _TDyOps { // vector. PetscErrorCode (*compute_error_norms)(void*,Vec,PetscReal*,PetscReal*); + // Functions used to define solver behavior. + // Material and boundary condition functions--we'll sort these out later. PetscErrorCode (*computeporosity)(TDy,PetscReal*,PetscReal*,void*); PetscErrorCode (*computepermeability)(TDy,PetscReal*,PetscReal*,void*); diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index f96ce75f..2dbb688c 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -13,7 +13,7 @@ typedef struct { PetscInt enthalpy_type; // Numerics settings - TDyMethod method; + TDyDiscretization discretization; TDyQuadratureType qtype; PetscInt mpfao_gmatrix_method; PetscInt mpfao_bc_type; diff --git a/include/tdycore.h b/include/tdycore.h index d83f804a..df2d906b 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -6,17 +6,36 @@ /* ---------------------------------------------------------------- */ +/// This type lists modes that identify the set of governing equations solved +/// by the dycore. typedef enum { - TPF=0, /* two point flux, classic finite volumes */ - MPFA_O, /* multipoint flux approximation - O method */ - MPFA_O_DAE, /* multipoint flux approximation - O method using DAE */ - MPFA_O_TRANSIENTVAR, /* multipoint flux approximation - O method using TS transient (conservative) approach */ - BDM, /* P0,BDM1 spaces, standard approach */ - WY, /* P0,BDM1 spaces, vertex quadrature, statically condensed */ - UNSPECIFIED_METHOD -} TDyMethod; + /// Richards equation + RICHARDS=0, + /// Non-isothermal flows + TH +} TDyMode; + +PETSC_EXTERN const char *const TDyModes[]; -PETSC_EXTERN const char *const TDyMethods[]; +/// This type enumerates discretizations supported by the dycore. +typedef enum { + /// two-point flux, classic finite volumes + TPF=0, + /// multi-point flux approximation - O method + MPFA_O, + /// multi-point flux approximation - O method using DAE + MPFA_O_DAE, + /// multipoint flux approximation - O method using TS transient (conservative) + /// approach + MPFA_O_TRANSIENTVAR, + /// finite element using P0, BDM1 spaces, standard approach + BDM, + /// finite element using P0,BDM1 spaces, vertex quadrature, statically + /// condensed + WY +} TDyDiscretization; + +PETSC_EXTERN const char *const TDyDiscretizations[]; typedef enum { MPFAO_GMATRIX_DEFAULT=0, /* default method to compute gmatrix for MPFA-O method */ @@ -40,12 +59,6 @@ typedef enum { PETSC_EXTERN const char *const TDyQuadratureTypes[]; -typedef enum { - RICHARDS=0, - TH, - UNSPECIFIED_MODE -} TDyMode; - typedef enum { TDySNES=0, TDyTS @@ -54,12 +67,12 @@ typedef enum { typedef enum { TDyCreated=0x0, TDyParametersInitialized=0x1, - TDyOptionsSet=0x2, - TDySetupFinished=0x4, + TDyModeSet=0x1<<1, + TDyDiscretizationSet=0x1<<2, + TDyOptionsSet=0x1<<3, + TDySetupFinished=0x1<<4, } TDySetupFlags; -PETSC_EXTERN const char *const TDyModes[]; - typedef void (*SpatialFunction)(PetscReal *x,PetscReal *f); /* returns f(x) */ typedef struct _p_TDy *TDy; @@ -83,7 +96,7 @@ PETSC_EXTERN PetscErrorCode TDyFinalize(void); PETSC_EXTERN PetscErrorCode TDyCreate(MPI_Comm, TDy*); PETSC_EXTERN PetscErrorCode TDySetMode(TDy,TDyMode); -PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyMethod); +PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyDiscretization); PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,PetscErrorCode(*)(void*, MPI_Comm, DM*)); PETSC_EXTERN PetscErrorCode TDySetFromOptions(TDy); PETSC_EXTERN PetscErrorCode TDySetup(TDy); diff --git a/src/mpfao/tdympfao.c b/src/mpfao/tdympfao.c index 1169e75d..ee348716 100644 --- a/src/mpfao/tdympfao.c +++ b/src/mpfao/tdympfao.c @@ -359,7 +359,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { PetscInt p, pStart, pEnd; PetscBool use_dae; - use_dae = (tdy->options.method == MPFA_O_DAE); + use_dae = (tdy->options.discretization == MPFA_O_DAE); ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); if (!use_dae) { if (tdy->options.mode == TH){ diff --git a/src/tdycore.c b/src/tdycore.c index 390424d8..b511f947 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -20,16 +20,15 @@ static char help[] = "TDycore \n\ #include #include -const char *const TDyMethods[] = { +const char *const TDyDiscretizations[] = { "TPF", "MPFA_O", "MPFA_O_DAE", "MPFA_O_TRANSIENTVAR", "BDM", "WY", - "UNSPECIFIED_METHOD", /* */ - "TDyMethod","TDY_METHOD_",NULL + "TDyDiscretization","TDY_DISCRETIZATION_",NULL }; const char *const TDyMPFAOGmatrixMethods[] = { @@ -196,8 +195,6 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { PetscFunctionBegin; TDyOptions *options = &tdy->options; - options->mode = UNSPECIFIED_MODE; - options->method = UNSPECIFIED_METHOD; options->gravity_constant = 9.8068; options->rho_type = WATER_DENSITY_CONSTANT; options->mu_type = WATER_VISCOSITY_CONSTANT; @@ -276,10 +273,9 @@ PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(void*, MPI PetscFunctionBegin; MPI_Comm comm; int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - if ((tdy->options.mode == UNSPECIFIED_MODE) || - (tdy->options.method == UNSPECIFIED_METHOD)) { + if ((tdy->setup_flags & TDyDiscretizationSet) == 0) { SETERRQ(comm,PETSC_ERR_USER, - "You must call TDySetMode and TDySetDiscretization before TDySetDMConstructor()"); + "You must call TDySetDiscretization before TDySetDMConstructor()"); } if (tdy->setup_flags & TDyOptionsSet) { SETERRQ(comm,PETSC_ERR_USER, @@ -630,7 +626,7 @@ PetscErrorCode TDyView(TDy tdy,PetscViewer viewer) { /// Sets options for the dycore based on command line arguments supplied by a /// user. TDySetFromOptions must be called before TDySetup, /// since the latter uses options specified by the former. -/// @param tdy The dycore instance +/// @param [inout] tdy The dycore instance PetscErrorCode TDySetFromOptions(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; @@ -712,8 +708,8 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Numerics options ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_method","Discretization method", - "TDySetDiscretization",TDyMethods, - (PetscEnum)options->method,(PetscEnum *)&options->method, + "TDySetDiscretization",TDyDiscretizations, + (PetscEnum)options->discretization,(PetscEnum *)&options->discretization, NULL); CHKERRQ(ierr); TDyQuadratureType qtype = FULL; ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods","TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&options->qtype,NULL); CHKERRQ(ierr); @@ -771,6 +767,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { PetscFunctionReturn(0); } +/// Performs setup, including allocation and configuration of any bookkeeping +/// data structures and the configuration of the dycore's DM object, which can +/// subsequently be passed to a solver (e.g. SNES or TS). This function must be +/// called after TDySetOptions. +/// @param [inout] tdy the dycore instance PetscErrorCode TDySetup(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; @@ -778,10 +779,9 @@ PetscErrorCode TDySetup(TDy tdy) { MPI_Comm comm; ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - if ((tdy->options.mode == UNSPECIFIED_MODE) || - (tdy->options.method == UNSPECIFIED_METHOD)) { + if ((tdy->setup_flags & TDyDiscretizationSet) == 0) { SETERRQ(comm,PETSC_ERR_USER, - "You must call TDySetMode and TDySetDiscretization before TDySetup()"); + "You must call TDySetDiscretization before TDySetup()"); } if ((tdy->setup_flags & TDyOptionsSet) == 0) { @@ -802,9 +802,9 @@ PetscErrorCode TDySetup(TDy tdy) { ierr = TDyRegressionInitialize(tdy); CHKERRQ(ierr); } if (tdy->options.output_mesh) { - if (tdy->options.method != MPFA_O) { + if (tdy->options.discretization != MPFA_O) { SETERRQ(comm,PETSC_ERR_USER, - "-tdy_output_mesh only supported for MPFA-O method"); + "-tdy_output_mesh only supported for MPFA-O discretization"); } } TDyExitProfilingStage("TDycore Setup"); @@ -813,14 +813,29 @@ PetscErrorCode TDySetup(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetDiscretization(TDy tdy,TDyMethod method) { +/// Sets the dycore's "mode", specifying the governing equations it solves. +/// @param [inout] tdy the dycore +/// @param [in] mode the selected mode +PetscErrorCode TDySetMode(TDy tdy, TDyMode mode) { + PetscValidPointer(tdy,1); + PetscFunctionBegin; + tdy->options.mode = mode; + tdy->setup_flags |= TDyModeSet; + PetscFunctionReturn(0); +} + +/// Sets the discretization used by the dycore. This must be called after +/// TDySetMode. +/// @param [inout] tdy the dycore +/// @param [in] discretization the selected discretization +PetscErrorCode TDySetDiscretization(TDy tdy, TDyDiscretization discretization) { PetscFunctionBegin; MPI_Comm comm; PetscErrorCode ierr; ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - if (tdy->options.mode == UNSPECIFIED_MODE) { + if ((tdy->setup_flags & TDyModeSet) == 0) { SETERRQ(comm,PETSC_ERR_USER, "You must call TDySetMode before TDySetDiscretization"); } @@ -829,30 +844,24 @@ PetscErrorCode TDySetDiscretization(TDy tdy,TDyMethod method) { // FIXME: should be renamed and altered to accept a context pointer and a DM // FIXME: as arguments. if (tdy->options.mode == RICHARDS) { - if (method == TPF) { + if (discretization == TPF) { tdy->ops->setup = TDyTPFRichardsSetup; - } else if (method == MPFA_O) { + } else if (discretization == MPFA_O) { tdy->ops->setup = TDyMPFAORichardsSetup; - } else if (method == MPFA_O_DAE) { + } else if (discretization == MPFA_O_DAE) { tdy->ops->setup = TDyMPFAORichardsSetup; - } else if (method == MPFA_O_TRANSIENTVAR) { + } else if (discretization == MPFA_O_TRANSIENTVAR) { tdy->ops->setup = TDyMPFAORichardsSetup; - } else if (method == BDM) { + } else if (discretization == BDM) { tdy->ops->setup = TDyBDMRichardsSetup; - } else if (method == WY) { + } else if (discretization == WY) { tdy->ops->setup = TDyWYRichardsSetup; } } else if (tdy->options.mode == TH) { - // How many methods do we support for TH? + // How many discretizations do we support for TH? } - tdy->options.method = method; - PetscFunctionReturn(0); -} - -PetscErrorCode TDySetMode(TDy tdy,TDyMode mode) { - PetscValidPointer(tdy,1); - PetscFunctionBegin; - tdy->options.mode = mode; + tdy->options.discretization = discretization; + tdy->setup_flags |= TDyDiscretizationSet; PetscFunctionReturn(0); } @@ -907,7 +916,7 @@ PetscErrorCode TDySetIFunction(TS ts,TDy tdy) { ierr = PetscSectionGetNumFields(sec, &num_fields); ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"IFunction not implemented for TPF"); break; @@ -947,7 +956,7 @@ PetscErrorCode TDySetIJacobian(TS ts,TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)ts,&comm); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"IJacobian not implemented for TPF"); break; @@ -995,7 +1004,7 @@ PetscErrorCode TDySetSNESFunction(SNES snes,TDy tdy) { PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"SNESFunction not implemented for TPF"); break; @@ -1032,7 +1041,7 @@ PetscErrorCode TDySetSNESJacobian(SNES snes,TDy tdy) { PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"SNESJacobian not implemented for TPF"); break; @@ -1062,7 +1071,7 @@ PetscErrorCode TDyComputeSystem(TDy tdy,Mat K,Vec F) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: ierr = TDyTPFComputeSystem(tdy,K,F); CHKERRQ(ierr); break; @@ -1166,9 +1175,9 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { } } - if ( (tdy->options.method == MPFA_O || - tdy->options.method == MPFA_O_DAE || - tdy->options.method == MPFA_O_TRANSIENTVAR)) { + if ( (tdy->options.discretization == MPFA_O || + tdy->options.discretization == MPFA_O_DAE || + tdy->options.discretization == MPFA_O_TRANSIENTVAR)) { PetscReal *p_vec_ptr, gz; TDyMesh *mesh = tdy->mesh; TDyCell *cells = &mesh->cells; @@ -1467,7 +1476,7 @@ PetscErrorCode TDyComputeErrorNorms(TDy tdy,Vec U,PetscReal *normp, PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: if(normp != NULL) { *normp = TDyTPFPressureNorm(tdy,U); } if(normv != NULL) { *normv = TDyTPFVelocityNorm(tdy,U); } @@ -1558,7 +1567,7 @@ PetscErrorCode TDyPreSolveSNESSolver(TDy tdy) { ierr = PetscObjectGetComm((PetscObject)tdy->dm,&comm); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { + switch (tdy->options.discretization) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"TDyPreSolveSNESSolver not implemented for TPF"); break; From 09f3706b3040d94a04f5b9414a461e094e38ce9d Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 10:06:33 -0700 Subject: [PATCH 08/23] Factored out setup functions. --- include/private/tdycoreimpl.h | 14 +- include/private/tdympfaoimpl.h | 6 +- include/private/tdythimpl.h | 2 +- include/private/tdytpfimpl.h | 9 + include/tdycore.h | 5 +- src/f90-mod/tdycoremod.F90 | 3 +- src/mpfao/tdympfao.c | 311 +++++++++++++++++++++++++-------- src/tdybdm.c | 4 +- src/tdycore.c | 47 +++-- src/tdydriver.c | 6 +- src/tdytimers.c | 31 +--- src/tdytpf.c | 3 +- src/tdywy.c | 3 +- 13 files changed, 310 insertions(+), 134 deletions(-) create mode 100644 include/private/tdytpfimpl.h diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index 68617afb..0383c6d5 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -24,18 +24,24 @@ struct _TDyOps { // Called by TDyDestroy to free implementation-specific resources. PetscErrorCode (*destroy)(void*); + // Creates a DM for a particular simulation (optional). + // We pass the dycore as the first argument here because we don't expect + // the caller to know implementation details. + PetscErrorCode (*create_dm)(MPI_Comm, DM*); + // Implements the view operation for the TDy implementation with the given // viewer. PetscErrorCode (*view)(void*, PetscViewer); - // Creates a DM for a particular simulation (optional) - PetscErrorCode (*create_dm)(void*, MPI_Comm, DM*); - // Called by TDySetFromOptions -- sets implementation-specific options // from command-line arguments. - PetscErrorCode (*set_from_options)(void*); + // FIXME: convert the first arg here to void* when we've moved specific data + // FIXME: out of TDy. + PetscErrorCode (*set_from_options)(TDy); // Called by TDySetup -- configures the DM for solvers. + // FIXME: we should convert the first argument here to void* when we've moved + // FIXME: all discretization-specific data out of TDy itself. PetscErrorCode (*setup)(TDy, DM); // Called by TDyComputeErrorNorms -- computes error norms given a solution diff --git a/include/private/tdympfaoimpl.h b/include/private/tdympfaoimpl.h index c87e247a..1f70eacb 100644 --- a/include/private/tdympfaoimpl.h +++ b/include/private/tdympfaoimpl.h @@ -3,7 +3,9 @@ #include -PETSC_INTERN PetscErrorCode TDyMPFAOInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_Setup(TDy, DM); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_DAE_Setup(TDy, DM); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_TRANSIENTVAR_Setup(TDy, DM); PETSC_INTERN PetscErrorCode TDyComputeGMatrixMPFAO(TDy); PETSC_INTERN PetscErrorCode TDyComputeGMatrixTPF(TDy); PETSC_INTERN PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy); @@ -21,7 +23,7 @@ PETSC_INTERN PetscErrorCode TDyMPFAOTransientVariable(TS,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOIFunction_TransientVariable(TS,PetscReal,Vec,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOSNESFunction(SNES,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOSNESJacobian(SNES,Vec,Mat,Mat,void*); -PETSC_INTERN PetscErrorCode TDyMPFAOSetFromOptions(TDy); +PETSC_INTERN PetscErrorCode TDyMPFAO_SetFromOptions(TDy); PETSC_INTERN PetscErrorCode TDyMPFAOSNESPreSolve(TDy); #endif diff --git a/include/private/tdythimpl.h b/include/private/tdythimpl.h index d0171237..fe269cf5 100644 --- a/include/private/tdythimpl.h +++ b/include/private/tdythimpl.h @@ -4,7 +4,7 @@ #include #include -PETSC_INTERN PetscErrorCode TDyTHInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyTH_MPFA_O_Setup(TDy, DM); PETSC_INTERN PetscErrorCode TDyTHTSPostStep(TS); PETSC_INTERN PetscErrorCode TDyTHSNESPostCheck(SNESLineSearch,Vec,Vec,Vec, PetscBool*,PetscBool*, diff --git a/include/private/tdytpfimpl.h b/include/private/tdytpfimpl.h new file mode 100644 index 00000000..c0699474 --- /dev/null +++ b/include/private/tdytpfimpl.h @@ -0,0 +1,9 @@ +#if !defined(TDYTPFIMPL_H) +#define TDYTPFIMPL_H + +#include +#include + +PETSC_INTERN PetscErrorCode TDyTPF_Setup(TDy, DM); + +#endif diff --git a/include/tdycore.h b/include/tdycore.h index df2d906b..aefd6528 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -97,7 +97,7 @@ PETSC_EXTERN PetscErrorCode TDyFinalize(void); PETSC_EXTERN PetscErrorCode TDyCreate(MPI_Comm, TDy*); PETSC_EXTERN PetscErrorCode TDySetMode(TDy,TDyMode); PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyDiscretization); -PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,PetscErrorCode(*)(void*, MPI_Comm, DM*)); +PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,PetscErrorCode(*)(MPI_Comm, DM*)); PETSC_EXTERN PetscErrorCode TDySetFromOptions(TDy); PETSC_EXTERN PetscErrorCode TDySetup(TDy); PETSC_EXTERN PetscErrorCode TDyDestroy(TDy*); @@ -186,16 +186,13 @@ PETSC_EXTERN PetscErrorCode TDyPostSolveSNESSolver(TDy,Vec); PETSC_EXTERN PetscErrorCode TDyOutputRegression(TDy,Vec); -PETSC_EXTERN PetscErrorCode TDyTPFInitialize(TDy); PETSC_EXTERN PetscErrorCode TDyTPFComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscReal TDyTPFPressureNorm(TDy,Vec); PETSC_EXTERN PetscReal TDyTPFVelocityNorm(TDy,Vec); PETSC_EXTERN PetscErrorCode TDyTPFCheckMeshSuitability(TDy); -PETSC_EXTERN PetscErrorCode TDyWYInitialize(TDy); PETSC_EXTERN PetscErrorCode TDyWYComputeSystem(TDy,Mat,Vec); -PETSC_EXTERN PetscErrorCode TDyBDMInitialize(TDy); PETSC_EXTERN PetscErrorCode TDyBDMComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscReal TDyBDMPressureNorm(TDy,Vec); PETSC_EXTERN PetscReal TDyBDMVelocityNorm(TDy,Vec); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index a9fb4293..0c691021 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -453,10 +453,9 @@ function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) end interface abstract interface - subroutine TDyDMConstructor(tdy, comm, dm, ierr) + subroutine TDyDMConstructor(comm, dm, ierr) use tdycoredef use petscdm - TDy :: tdy MPI_Comm :: comm DM :: dm PetscErrorCode :: ierr diff --git a/src/mpfao/tdympfao.c b/src/mpfao/tdympfao.c index ee348716..06d4a20f 100644 --- a/src/mpfao/tdympfao.c +++ b/src/mpfao/tdympfao.c @@ -307,15 +307,217 @@ PetscErrorCode TDyMPFAO_AllocateMemoryForEnergySourceSinkValues(TDy tdy) { TDY_STOP_FUNCTION_TIMER() PetscFunctionReturn(0); } - /* -------------------------------------------------------------------------- */ -PetscErrorCode TDyMPFAOInitialize(TDy tdy) { +//----------------- +// Setup functions +//----------------- + +// There's a lot of duplicated code here for the moment. We'll factor out the +// repeated stuff when we move the data out of the TDy struct and into +// discretization-specific types. +// Setup function for Richards + TPF +PetscErrorCode TDyRichards_TPF_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +// Setup function for Richards + MPFA_O +PetscErrorCode TDyRichards_MPFA_O_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + + PetscInt dim; + ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); + if (dim == 2) { + SETERRQ(comm,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); + } + + // Allocate mesh data. + tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); + ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); + ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), + &(tdy->vel )); CHKERRQ(ierr); + ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscInt), + &(tdy->vel_count)); CHKERRQ(ierr); + ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; + ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); + + // Set up the section, 1 dof per cell + PetscSection sec; + PetscInt p, pStart, pEnd; + ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); + ierr = PetscSectionSetNumFields(sec, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + + ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); + ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); + for(p=pStart; pmesh->num_faces; + nLocalCells = TDyMeshGetNumberOfLocalCells(tdy->mesh); + nNonLocalFaces = TDyMeshGetNumberOfNonLocalFacess(tdy->mesh); + nNonInternalFaces = TDyMeshGetNumberOfNonInternalFacess(tdy->mesh); + + nrow = 4*nFaces; + ncol = nLocalCells + nNonLocalFaces + nNonInternalFaces; + nz = tdy->nfv; + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->TtimesP_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); + ierr = VecZeroEntries(tdy->GravDisVec); + } + + // Set up material models. + if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } + if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} + + // why must these be placed after SetPermeabilityFromFunction()? + ierr = TDyComputeGMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeTransmissibilityMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeGravityDiscretization(tdy); CHKERRQ(ierr); + + ierr = TDyMPFAO_AllocateMemoryForBoundaryValues(tdy); CHKERRQ(ierr); + ierr = TDyMPFAO_AllocateMemoryForSourceSinkValues(tdy); CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +// Setup function for Richards + MPFA_O_DAE +PetscErrorCode TDyRichards_MPFA_O_DAE_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + PetscFunctionReturn(0); + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + + PetscInt dim; + ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); + if (dim == 2) { + SETERRQ(comm,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); + } + + // Allocate mesh data. + tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); + ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); + ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), + &(tdy->vel )); CHKERRQ(ierr); + ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscInt), + &(tdy->vel_count)); CHKERRQ(ierr); + ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; + ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); + + // Set up the section, 1 dof per cell + PetscSection sec; + PetscInt p, pStart, pEnd; + ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); + ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 1, "LiquidMass"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); + + ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); + ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); + for(p=pStart; pmesh->num_faces; + nLocalCells = TDyMeshGetNumberOfLocalCells(tdy->mesh); + nNonLocalFaces = TDyMeshGetNumberOfNonLocalFacess(tdy->mesh); + nNonInternalFaces = TDyMeshGetNumberOfNonInternalFacess(tdy->mesh); + + nrow = 4*nFaces; + ncol = nLocalCells + nNonLocalFaces + nNonInternalFaces; + nz = tdy->nfv; + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->TtimesP_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); + ierr = VecZeroEntries(tdy->GravDisVec); + } + + // Set up material models. + if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } + if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} + + // why must these be placed after SetPermeabilityFromFunction()? + ierr = TDyComputeGMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeTransmissibilityMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeGravityDiscretization(tdy); CHKERRQ(ierr); + + ierr = TDyMPFAO_AllocateMemoryForBoundaryValues(tdy); CHKERRQ(ierr); + ierr = TDyMPFAO_AllocateMemoryForSourceSinkValues(tdy); CHKERRQ(ierr); + +} + +// Setup function for Richards + MPFA_O_TRANSIENTVAR +PetscErrorCode TDyRichards_MPFA_O_TRANSIENTVAR_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + // This is essentially the same as the MPFA-O one. + ierr = TDyRichards_MPFA_O_Setup(tdy, dm); + PetscFunctionReturn(0); +} + +// Setup function for TH + MPFA-O +PetscErrorCode TDyTH_MPFA_O_Setup(TDy tdy, DM dm) { PetscErrorCode ierr; - MPI_Comm comm; - DM dm = tdy->dm; - PetscInt nrow,ncol,nsubcells; PetscFunctionBegin; TDY_START_FUNCTION_TIMER() @@ -326,16 +528,15 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); } + MPI_Comm comm; ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + // Allocate mesh data. tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); - ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); - if (tdy->options.mode == TH) { - ierr = TDyAllocate_RealArray_3D(&tdy->Temp_Trans, tdy->mesh->num_vertices, - tdy->nfv, tdy->nfv); CHKERRQ(ierr); - } + ierr = TDyAllocate_RealArray_3D(&tdy->Temp_Trans, tdy->mesh->num_vertices, + tdy->nfv, tdy->nfv); CHKERRQ(ierr); ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), &(tdy->vel )); CHKERRQ(ierr); ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); @@ -343,74 +544,41 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { &(tdy->vel_count)); CHKERRQ(ierr); ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); - nsubcells = 8; - nrow = 3; - ncol = 3; + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, nsubcells, nrow, ncol); CHKERRQ(ierr); - if (tdy->options.mode == TH) { - ierr = TDyAllocate_RealArray_4D(&tdy->Temp_subc_Gmatrix, tdy->mesh->num_cells, - nsubcells, nrow, ncol); CHKERRQ(ierr); - } + ierr = TDyAllocate_RealArray_4D(&tdy->Temp_subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); - /* Set up the section, 1 dof per cell */ + // Set up the section, 1 dof per cell PetscSection sec; PetscInt p, pStart, pEnd; - PetscBool use_dae; - - use_dae = (tdy->options.discretization == MPFA_O_DAE); ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); - if (!use_dae) { - if (tdy->options.mode == TH){ - ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 1, "LiquidTemperature"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); - } else { - ierr = PetscSectionSetNumFields(sec, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - } - } else { - if (tdy->options.mode == TH) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_DAE"); - } - ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 1, "LiquidMass"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); - } + ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 1, "LiquidTemperature"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); for(p=pStart; poptions.mode == TH) { - ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); - ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); - } - if (use_dae) { - if (tdy->options.mode == TH) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_DAE"); - } - ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); - ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); - } + ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); + ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); } ierr = PetscSectionSetUp(sec); CHKERRQ(ierr); ierr = DMSetSection(dm,sec); CHKERRQ(ierr); ierr = PetscSectionViewFromOptions(sec, NULL, "-layout_view"); CHKERRQ(ierr); ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr); ierr = DMSetBasicAdjacency(dm,PETSC_TRUE,PETSC_TRUE); CHKERRQ(ierr); - //ierr = DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);CHKERRQ(ierr); - //ierr = DMPlexSetAdjacencyUseClosure(dm,PETSC_TRUE);CHKERRQ(ierr); + // Build the mesh. ierr = TDyBuildMesh(tdy); CHKERRQ(ierr); - { PetscInt nLocalCells, nFaces, nNonLocalFaces, nNonInternalFaces; PetscInt nrow, ncol, nz; @@ -430,25 +598,22 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); ierr = VecZeroEntries(tdy->GravDisVec); - if (tdy->options.mode == TH) { - ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Temp_Trans_mat); CHKERRQ(ierr); - ierr = MatSetOption(tdy->Temp_Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); - ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->Temp_P_vec); - ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->Temp_TtimesP_vec); - } + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Temp_Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Temp_Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->Temp_P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->Temp_TtimesP_vec); } + // Set up material models. if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} - if (tdy->options.mode == TH){ - if (tdy->ops->computethermalconductivity) {ierr = SetThermalConductivityFromFunction(tdy); CHKERRQ(ierr);} - if (tdy->ops->computesoildensity) { - ierr = SetSoilDensityFromFunction(tdy); CHKERRQ(ierr); - } - if (tdy->ops->computesoilspecificheat) { - printf("SetSoilSpecificHeatFromFunction\n"); - ierr = SetSoilSpecificHeatFromFunction(tdy); CHKERRQ(ierr); - } + if (tdy->ops->computethermalconductivity) {ierr = SetThermalConductivityFromFunction(tdy); CHKERRQ(ierr);} + if (tdy->ops->computesoildensity) { + ierr = SetSoilDensityFromFunction(tdy); CHKERRQ(ierr); + } + if (tdy->ops->computesoilspecificheat) { + printf("SetSoilSpecificHeatFromFunction\n"); + ierr = SetSoilSpecificHeatFromFunction(tdy); CHKERRQ(ierr); } // why must these be placed after SetPermeabilityFromFunction()? @@ -469,7 +634,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { } /* -------------------------------------------------------------------------- */ -PetscErrorCode TDyMPFAOSetFromOptions(TDy tdy) { +PetscErrorCode TDyMPFAO_SetFromOptions(TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() diff --git a/src/tdybdm.c b/src/tdybdm.c index 7b865a7f..156880f4 100644 --- a/src/tdybdm.c +++ b/src/tdybdm.c @@ -1,4 +1,5 @@ #include +#include #include #include @@ -133,7 +134,7 @@ void HdivBasisHex(const PetscReal *x,PetscReal *B,PetscReal *DF,PetscReal J) { TDY_STOP_FUNCTION_TIMER() } -PetscErrorCode TDyBDMInitialize(TDy tdy) { +PetscErrorCode TDyBDM_Setup(TDy tdy, DM dm) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() PetscErrorCode ierr; @@ -142,7 +143,6 @@ PetscErrorCode TDyBDMInitialize(TDy tdy) { PetscSection sec; PetscInt d,dim,dofs_per_face = 1; PetscBool found; - DM dm = tdy->dm; ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); /* Get plex limits */ diff --git a/src/tdycore.c b/src/tdycore.c index b511f947..d5f4ec5b 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -9,7 +9,11 @@ static char help[] = "TDycore \n\ #include #include #include +#include #include +#include +#include +#include #include #include #include @@ -269,7 +273,7 @@ PetscErrorCode TDyCreate(MPI_Comm comm, TDy *_tdy) { /// distributed appropriately. This function must be called before /// TDySetFromOptions. /// @param [in] dm_func A function that creates a given DM and returns an error -PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(void*, MPI_Comm, DM*)) { +PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(MPI_Comm, DM*)) { PetscFunctionBegin; MPI_Comm comm; int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); @@ -388,7 +392,7 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { } else { if (tdy->ops->create_dm) { // We've been instructed to create a DM ourselves. - ierr = tdy->ops->create_dm(tdy->context, comm, &dm); + ierr = tdy->ops->create_dm(comm, &dm); } else { ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); } @@ -748,6 +752,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Other options ierr = PetscOptionsBool("-tdy_regression_test","Enable output of a regression file","",options->regression_testing,&(options->regression_testing),NULL); CHKERRQ(ierr); + // Mode/discretization-specific options. + if (tdy->ops->set_from_options) { + tdy->ops->set_from_options(tdy); + } + // Wrap up and indicate that options are set. ierr = PetscOptionsEnd(); CHKERRQ(ierr); tdy->setup_flags |= TDyOptionsSet; @@ -775,6 +784,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { PetscErrorCode TDySetup(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; + TDY_START_FUNCTION_TIMER() MPI_Comm comm; ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); @@ -787,11 +797,11 @@ PetscErrorCode TDySetup(TDy tdy) { if ((tdy->setup_flags & TDyOptionsSet) == 0) { SETERRQ(comm,PETSC_ERR_USER,"You must call TDySetFromOptions before TDySetup()"); } - TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); // Perform implementation-specific setup. - ierr = tdy->ops->setup(tdy->context, tdy->dm); CHKERRQ(ierr); + // FIXME: Pass tdy->context as first argument when we break data out of TDy. + ierr = tdy->ops->setup(tdy, tdy->dm); CHKERRQ(ierr); // Record metadata for scaling studies. TDySetTimingMetadata(tdy); @@ -840,25 +850,34 @@ PetscErrorCode TDySetDiscretization(TDy tdy, TDyDiscretization discretization) { "You must call TDySetMode before TDySetDiscretization"); } - // FIXME: All of these functions are current called TDy*Initialize, and - // FIXME: should be renamed and altered to accept a context pointer and a DM - // FIXME: as arguments. + // Set function pointers for operations. if (tdy->options.mode == RICHARDS) { if (discretization == TPF) { - tdy->ops->setup = TDyTPFRichardsSetup; + tdy->ops->set_from_options = NULL; // FIXME: define this and move things here! + tdy->ops->setup = TDyTPF_Setup; } else if (discretization == MPFA_O) { - tdy->ops->setup = TDyMPFAORichardsSetup; + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_Setup; } else if (discretization == MPFA_O_DAE) { - tdy->ops->setup = TDyMPFAORichardsSetup; + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_DAE_Setup; } else if (discretization == MPFA_O_TRANSIENTVAR) { - tdy->ops->setup = TDyMPFAORichardsSetup; + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_TRANSIENTVAR_Setup; } else if (discretization == BDM) { - tdy->ops->setup = TDyBDMRichardsSetup; + tdy->ops->set_from_options = NULL; + tdy->ops->setup = TDyBDM_Setup; } else if (discretization == WY) { - tdy->ops->setup = TDyWYRichardsSetup; + tdy->ops->set_from_options = NULL; + tdy->ops->setup = TDyWY_Setup; } } else if (tdy->options.mode == TH) { - // How many discretizations do we support for TH? + if (discretization == MPFA_O) { + tdy->ops->setup = TDyTH_MPFA_O_Setup; + } else { + SETERRQ(comm,PETSC_ERR_USER, + "The TH mode does not support the selected discretization!"); + } } tdy->options.discretization = discretization; tdy->setup_flags |= TDyDiscretizationSet; diff --git a/src/tdydriver.c b/src/tdydriver.c index 74f713b1..d025a848 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -22,7 +22,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver currently only supports 3D"); } - switch(tdy->options.method) { + switch(tdy->options.discretization) { case TPF: case MPFA_O: break; @@ -30,7 +30,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case MPFA_O_TRANSIENTVAR: case BDM: case WY: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver not supported for specified method."); + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver not supported for specified discretization."); break; } @@ -161,7 +161,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case TDyTS: ierr = TSSetPostStep(ts,TDyTHTSPostStep); CHKERRQ(ierr); } - ierr = TDyTHInitialize(tdy); CHKERRQ(ierr); + ierr = TDySetup(tdy); CHKERRQ(ierr); break; } PetscPrintf(PETSC_COMM_WORLD,"tdy->ti->time_integration_method = %d\n", diff --git a/src/tdytimers.c b/src/tdytimers.c index b50c5c8f..2b0e5eb8 100644 --- a/src/tdytimers.c +++ b/src/tdytimers.c @@ -11,8 +11,8 @@ khash_t(TDY_PROFILING_STAGE_MAP)* TDY_PROFILING_STAGES = NULL; // Timing metadata used by tdyperfplot and other tools. typedef struct { - TDyMethod method; TDyMode mode; + TDyDiscretization discretization; int num_cells; int num_proc; } TimingMetadata; @@ -91,8 +91,8 @@ PetscErrorCode TDySetTimingMetadata(TDy tdy) { iter = kh_put(TDY_PROFILING_MD_MAP, TDY_PROFILING_METADATA, tdy_addr, &retval); kh_val(TDY_PROFILING_METADATA, iter) = md; } - md->method = tdy->options.method; md->mode = tdy->options.mode; + md->discretization = tdy->options.discretization; if (tdy->mesh != NULL) { md->num_cells = TDyMeshGetNumberOfLocalCells(tdy->mesh); } else { @@ -136,31 +136,12 @@ PetscErrorCode TDyWriteTimingProfile(const char* filename) { int num_cells; MPI_Reduce(&(md->num_cells), &num_cells, 1, MPI_INT, MPI_SUM, 0, PETSC_COMM_WORLD); - const char* method_name; - if (md->method == TPF) { - method_name = "TPF"; - } else if (md->method == MPFA_O) { - method_name = "MPFA_O"; - } else if (md->method == MPFA_O_DAE) { - method_name = "MPFA_O_DAE"; - } else if (md->method == MPFA_O_TRANSIENTVAR) { - method_name = "MPFA_O_TRANSIENTVAR"; - } else if (md->method == BDM) { - method_name = "BDM"; - } else { // (md->method == BDM) - method_name = "WY"; - } - const char* mode_name; - if (md->mode == RICHARDS) { - mode_name = "RICHARDS"; - } else { // (md->mode == TH) - mode_name = "TH"; - } + const char* mode_name = TDyModes[md->mode]; + const char* disc_name = TDyDiscretizations[md->discretization]; FILE* f = fopen("tdycore_profile.csv", "a"); fprintf(f, "METADATA\n"); - fprintf(f, "Method,Mode,NumProc,NumCells\n"); - fprintf(f, "%s,%s,%d,%d", method_name, mode_name, - md->num_proc, num_cells); + fprintf(f, "Mode,Discretiztion,NumProc,NumCells\n"); + fprintf(f, "%s,%s,%d,%d", mode_name, disc_name, md->num_proc, num_cells); fclose(f); } else { // rank > 0 MPI_Reduce(&md->num_cells, NULL, 1, MPI_INT, MPI_SUM, diff --git a/src/tdytpf.c b/src/tdytpf.c index da0509bf..964e1348 100644 --- a/src/tdytpf.c +++ b/src/tdytpf.c @@ -7,12 +7,11 @@ PETSC_STATIC_INLINE PetscScalar Dot(PetscInt dim, const PetscScalar *x, const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;} PETSC_STATIC_INLINE PetscReal Norm(PetscInt dim, const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(Dot(dim,x,x)));} -PetscErrorCode TDyTPFInitialize(TDy tdy) { +PetscErrorCode TDyTPF_Setup(TDy tdy, DM dm) { PetscErrorCode ierr; MPI_Comm comm; PetscInt dim,c,cStart,cEnd,pStart,pEnd; PetscSection sec; - DM dm = tdy->dm; PetscFunctionBegin; TDY_START_FUNCTION_TIMER() diff --git a/src/tdywy.c b/src/tdywy.c index 7b80db00..aa9d3802 100644 --- a/src/tdywy.c +++ b/src/tdywy.c @@ -217,7 +217,7 @@ PetscErrorCode TDyWYLocalElementCompute(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDyWYInitialize(TDy tdy) { +PetscErrorCode TDyWY_Setup(TDy tdy, DM dm) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() PetscErrorCode ierr; @@ -225,7 +225,6 @@ PetscErrorCode TDyWYInitialize(TDy tdy) { PetscInt i,dim,ncv,nfv,c,cStart,cEnd,f,fStart,fEnd,vStart,vEnd,p,pStart,pEnd; PetscInt closureSize, *closure; PetscSection sec; - DM dm = tdy->dm; MaterialProp *matprop = tdy->matprop; ierr = PetscObjectGetComm((PetscObject)dm,&comm); CHKERRQ(ierr); From 7d50c9072fc3c1181f74743d6283ec0964888fe2 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 10:07:14 -0700 Subject: [PATCH 09/23] Forgot to add a couple of files. --- include/private/tdybdmimpl.h | 9 +++++++++ include/private/tdywyimpl.h | 9 +++++++++ 2 files changed, 18 insertions(+) create mode 100644 include/private/tdybdmimpl.h create mode 100644 include/private/tdywyimpl.h diff --git a/include/private/tdybdmimpl.h b/include/private/tdybdmimpl.h new file mode 100644 index 00000000..39e8328c --- /dev/null +++ b/include/private/tdybdmimpl.h @@ -0,0 +1,9 @@ +#if !defined(TDYBDMIMPL_H) +#define TDYBDMIMPL_H + +#include +#include + +PETSC_INTERN PetscErrorCode TDyBDM_Setup(TDy, DM); + +#endif diff --git a/include/private/tdywyimpl.h b/include/private/tdywyimpl.h new file mode 100644 index 00000000..b24f2c52 --- /dev/null +++ b/include/private/tdywyimpl.h @@ -0,0 +1,9 @@ +#if !defined(TDYWYIMPL_H) +#define TDYWYIMPL_H + +#include +#include + +PETSC_INTERN PetscErrorCode TDyWY_Setup(TDy, DM); + +#endif From 5db4d16a448181e96b0d20e2166c64dc4a727a0a Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 10:50:05 -0700 Subject: [PATCH 10/23] Managed to get the steady demo working without TDySetDM. --- demo/steady/steady.c | 44 ++++++++++++++++------------------- demo/steady/steady.cfg | 30 ++++++++++++------------ include/private/tdycoreimpl.h | 5 +++- include/tdycore.h | 2 +- src/f90-mod/tdycoremod.F90 | 17 ++++++++------ src/tdycore.c | 29 ++++++++++++++++++----- 6 files changed, 73 insertions(+), 54 deletions(-) diff --git a/demo/steady/steady.c b/demo/steady/steady.c index 4e8274c6..ca2634cc 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -557,36 +557,35 @@ PetscErrorCode GeometryColumn(DM dm){ // This set of globals is used by CreateDM below to create a DM for this demo. typedef struct DMOptions { - PetscInt dim; // Dimension of DM (2 or 3) - PetscInt N; // Number of cells on a side - PetscBool perturb; // whether to perturb randomly (as opposed to smoothly) - PetscBool exo; // whether to load a named exodus file - PetscBool column; // column mesh? - char exofile[256]; // name of the exodus file to load + PetscInt dim; // Dimension of DM (2 or 3) + PetscInt N; // Number of cells on a side + PetscBool perturb; // whether to perturb randomly (as opposed to smoothly) + PetscBool exo; // whether to load a named exodus file + PetscBool column; // column mesh? + const char* exofile; // name of the exodus file to load } DMOptions; -static DMOptions dm_options_; - // This function creates a DM specifically for this demo. Overrides are applied // to the resulting DM with TDySetFromOptions. -PetscErrorCode CreateDM(void* context, MPI_Comm comm, DM* dm) { +PetscErrorCode CreateDM(MPI_Comm comm, void* context, DM* dm) { int ierr; + DMOptions* options = context; - PetscInt N = dm_options_.N; - if(dm_options_.exo) { - ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD, dm_options_.exofile, + PetscInt N = options->N; + if(options->exo) { + ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD, options->exofile, PETSC_TRUE,dm); CHKERRQ(ierr); } else { PetscInt Nx=N,Ny=N,Nz=N; PetscReal Lx=1,Ly=1,Lz=1; - if(dm_options_.column){ + if(options->column){ Nx = 1; Ny = 1; Nz = N; Lx = 10; Ly = 10; Lz = 1; } const PetscInt faces[3] = {Nx ,Ny ,Nz }; const PetscReal lower[3] = {0.0,0.0,0.0}; const PetscReal upper[3] = {Lx ,Ly ,Lz }; - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dm_options_.dim,PETSC_FALSE, + ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,options->dim,PETSC_FALSE, faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); } } @@ -611,10 +610,12 @@ int main(int argc, char **argv) { MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy; ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); + ierr = TDySetMode(tdy,RICHARDS); CHKERRQ(ierr); ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-dim","Problem dimension","",dim,&dim,NULL); CHKERRQ(ierr); + ierr = PetscOptionsInt("-N","Number of elements in 1D","",N,&N,NULL); CHKERRQ(ierr); ierr = PetscOptionsInt("-problem","Problem number","",problem,&problem,NULL); CHKERRQ(ierr); ierr = PetscOptionsReal("-alpha","Permeability scaling","",alpha,&alpha,NULL); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code","Code passed on successful completion","",successful_exit_code,&successful_exit_code,NULL); @@ -630,16 +631,11 @@ int main(int argc, char **argv) { ierr = PetscStrcasecmp(paper,"wheeler2012",&wheeler2012); CHKERRQ(ierr); ierr = PetscStrcasecmp(paper,"column",&column); CHKERRQ(ierr); - // Copy DM-related globals into place for use with CreateDM. - dm_options_.N = N; - dm_options_.dim = dim; - dm_options_.perturb = perturb; - dm_options_.column = column; - dm_options_.exo = exo; - if (exo) strcpy(dm_options_.exofile, exofile); - - // Specify a special DM to be constructed for this demo. - ierr = TDySetDMConstructor(tdy, CreateDM); CHKERRQ(ierr); + // Specify a special DM to be constructed for this demo, and pass it the + // relevant options. + DMOptions dm_options = {.N = N, .dim = dim, .perturb = perturb, + .exo = exo, .exofile = exofile}; + ierr = TDySetDMConstructor(tdy, &dm_options, CreateDM); CHKERRQ(ierr); // Apply overrides. ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); diff --git a/demo/steady/steady.cfg b/demo/steady/steady.cfg index c203d220..e107310a 100644 --- a/demo/steady/steady.cfg +++ b/demo/steady/steady.cfg @@ -20,48 +20,48 @@ standard= pressure = 1.0e-12 relative [steady-wy-prob1] -input_arguments=-tdy_method wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 +input_arguments=-tdy_disc wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 [steady-wy-prob2] -input_arguments=-tdy_method wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 +input_arguments=-tdy_disc wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 [steady-wy-prob3] -input_arguments=-tdy_method wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 +input_arguments=-tdy_disc wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 [steady-wy-prob4] -input_arguments=-tdy_method wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 +input_arguments=-tdy_disc wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 [steady-wy-prob1-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D +input_arguments=-tdy_disc wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D [steady-wy-prob2-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D +input_arguments=-tdy_disc wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D [steady-wy-prob3-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D +input_arguments=-tdy_disc wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D [steady-tpf-prob1] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 +input_arguments=-tdy_disc tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 [steady-tpf-prob2] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 +input_arguments=-tdy_disc tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 [steady-tpf-prob3] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 +input_arguments=-tdy_disc tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 [steady-tpf-prob4] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 +input_arguments=-tdy_disc tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 [steady-bdm-prob1] pressure = 1.0e-12 absolute -input_arguments=-tdy_method bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu +input_arguments=-tdy_disc bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu [steady-bdm-prob2] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu +input_arguments=-tdy_disc bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu [steady-bdm-prob3] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu +input_arguments=-tdy_disc bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu [steady-bdm-prob4] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu +input_arguments=-tdy_disc bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index 0383c6d5..d2df81e2 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -27,7 +27,7 @@ struct _TDyOps { // Creates a DM for a particular simulation (optional). // We pass the dycore as the first argument here because we don't expect // the caller to know implementation details. - PetscErrorCode (*create_dm)(MPI_Comm, DM*); + PetscErrorCode (*create_dm)(MPI_Comm, void*, DM*); // Implements the view operation for the TDy implementation with the given // viewer. @@ -77,6 +77,9 @@ struct _p_TDy { // configured DM dm; + // Contextual information passed to create_dm (if given). + void* create_dm_context; + // We'll likely get rid of this. TDyTimeIntegrator ti; diff --git a/include/tdycore.h b/include/tdycore.h index aefd6528..4c43d6ed 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -97,7 +97,7 @@ PETSC_EXTERN PetscErrorCode TDyFinalize(void); PETSC_EXTERN PetscErrorCode TDyCreate(MPI_Comm, TDy*); PETSC_EXTERN PetscErrorCode TDySetMode(TDy,TDyMode); PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyDiscretization); -PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,PetscErrorCode(*)(MPI_Comm, DM*)); +PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,void*,PetscErrorCode(*)(MPI_Comm, void*, DM*)); PETSC_EXTERN PetscErrorCode TDySetFromOptions(TDy); PETSC_EXTERN PetscErrorCode TDySetup(TDy); PETSC_EXTERN PetscErrorCode TDyDestroy(TDy*); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index 0c691021..3e35b1ad 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -453,12 +453,14 @@ function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) end interface abstract interface - subroutine TDyDMConstructor(comm, dm, ierr) + subroutine TDyDMConstructor(comm, context, dm, ierr) + use, intrinsic :: iso_c_binding use tdycoredef use petscdm - MPI_Comm :: comm - DM :: dm - PetscErrorCode :: ierr + MPI_Comm :: comm + type(c_ptr), value :: context + DM :: dm + PetscErrorCode :: ierr end subroutine end interface @@ -504,16 +506,17 @@ subroutine TDySetDMConstructor(tdy, dm_ctor, ierr) PetscErrorCode :: ierr interface - function SetDMConstructor(tdy, c_func) bind (c, name="TDySetDMConstructor") result(ierr) + function SetDMConstructor(tdy, context, dm_ctor) bind (c, name="TDySetDMConstructor") result(ierr) use, intrinsic :: iso_c_binding implicit none type(c_ptr) :: tdy - type(c_funptr) :: c_func + type(c_ptr) :: context + type(c_funptr) :: dm_ctor integer(c_int) :: ierr end function end interface - ierr = SetDMConstructor(c_loc(tdy), c_funloc(dm_ctor)) + ierr = SetDMConstructor(c_loc(tdy), c_null_ptr, c_funloc(dm_ctor)) end subroutine subroutine TDySelectPorosityFunction(tdy, name, ierr) diff --git a/src/tdycore.c b/src/tdycore.c index d5f4ec5b..d03a2d55 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -272,8 +272,14 @@ PetscErrorCode TDyCreate(MPI_Comm comm, TDy *_tdy) { /// DMSetFromOptions is called to apply overrides, and then the DM is /// distributed appropriately. This function must be called before /// TDySetFromOptions. -/// @param [in] dm_func A function that creates a given DM and returns an error -PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(MPI_Comm, DM*)) { +/// @param [inout] tdy the dycore +/// @param [in] context a pointer to contextual information that can be used by +/// dm_func to create a DM. This pointer is not managed +/// by the dycore. +/// @param [in] dm_func A function that, given an MPI communicator and a context +/// pointer, creates a given DM and returns an error +PetscErrorCode TDySetDMConstructor(TDy tdy, void* context, + PetscErrorCode (*dm_func)(MPI_Comm, void*, DM*)) { PetscFunctionBegin; MPI_Comm comm; int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); @@ -285,6 +291,7 @@ PetscErrorCode TDySetDMConstructor(TDy tdy, PetscErrorCode (*dm_func)(MPI_Comm, SETERRQ(comm,PETSC_ERR_USER, "You must call TDyDefineDM before TDySetFromOptions()"); } + tdy->create_dm_context = context; tdy->ops->create_dm = dm_func; PetscFunctionReturn(0); } @@ -392,7 +399,7 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { } else { if (tdy->ops->create_dm) { // We've been instructed to create a DM ourselves. - ierr = tdy->ops->create_dm(comm, &dm); + ierr = tdy->ops->create_dm(comm, tdy->create_dm_context, &dm); } else { ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); } @@ -671,10 +678,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Model options + TDyMode mode; 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); + (PetscEnum *)&mode, NULL); CHKERRQ(ierr); ierr = PetscOptionsReal("-tdy_gravity", "Magnitude of gravity vector", NULL, options->gravity_constant, &options->gravity_constant, NULL); CHKERRQ(ierr); @@ -710,10 +718,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options + TDyDiscretization discretization; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); - ierr = PetscOptionsEnum("-tdy_method","Discretization method", + ierr = PetscOptionsEnum("-tdy_disc","Discretization", "TDySetDiscretization",TDyDiscretizations, - (PetscEnum)options->discretization,(PetscEnum *)&options->discretization, + (PetscEnum)options->discretization,(PetscEnum *)&discretization, NULL); CHKERRQ(ierr); TDyQuadratureType qtype = FULL; ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods","TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&options->qtype,NULL); CHKERRQ(ierr); @@ -752,6 +761,14 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Other options ierr = PetscOptionsBool("-tdy_regression_test","Enable output of a regression file","",options->regression_testing,&(options->regression_testing),NULL); CHKERRQ(ierr); + // Override the mode and/or discretization if needed. + if (options->mode != mode) { + TDySetMode(tdy, mode); + } + if (options->discretization != discretization) { + TDySetDiscretization(tdy, discretization); + } + // Mode/discretization-specific options. if (tdy->ops->set_from_options) { tdy->ops->set_from_options(tdy); From faba44db7d393bacb1e0e4ade0642f41fa5b5b81 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 11:13:32 -0700 Subject: [PATCH 11/23] Adjustments to comments. --- demo/steady/steady.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/demo/steady/steady.c b/demo/steady/steady.c index ca2634cc..c4a43096 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -273,7 +273,6 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { PetscScalar *coords; PetscInt v,vStart,vEnd,offset,value,dim; ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - /* this is the 'marker' label which marks boundary entities */ ierr = DMGetLabel(dm, "boundary", &label); CHKERRQ(ierr); ierr = DMGetCoordinateSection(dm, &coordSection); CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); @@ -291,7 +290,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { coords[offset+1] += r*PetscSinReal(t); } } else { - /* This is because 'marker' is broken in 3D */ + /* This is because 'boundary' is broken in 3D */ if(coords[offset] > 0 && coords[offset] < 1 && coords[offset+1] > 0 && coords[offset+1] < 1 && coords[offset+2] > 0 && coords[offset+2] < 1) { @@ -555,7 +554,7 @@ PetscErrorCode GeometryColumn(DM dm){ PetscFunctionReturn(0); } -// This set of globals is used by CreateDM below to create a DM for this demo. +// This data is used by CreateDM below to create a DM for this demo. typedef struct DMOptions { PetscInt dim; // Dimension of DM (2 or 3) PetscInt N; // Number of cells on a side From bad0066a96c1e490b18568605be0360502899d28 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 11:26:01 -0700 Subject: [PATCH 12/23] Changed -tdy_disc -> -tdy_discretization --- demo/steady/steady.cfg | 30 +++++++++++++++--------------- src/tdycore.c | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/demo/steady/steady.cfg b/demo/steady/steady.cfg index e107310a..6db2cc9f 100644 --- a/demo/steady/steady.cfg +++ b/demo/steady/steady.cfg @@ -20,48 +20,48 @@ standard= pressure = 1.0e-12 relative [steady-wy-prob1] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 [steady-wy-prob2] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 [steady-wy-prob3] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 [steady-wy-prob4] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 [steady-wy-prob1-3D] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D [steady-wy-prob2-3D] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D [steady-wy-prob3-3D] -input_arguments=-tdy_disc wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D [steady-tpf-prob1] -input_arguments=-tdy_disc tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 +input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 [steady-tpf-prob2] -input_arguments=-tdy_disc tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 +input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 [steady-tpf-prob3] -input_arguments=-tdy_disc tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 +input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 [steady-tpf-prob4] -input_arguments=-tdy_disc tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 +input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 [steady-bdm-prob1] pressure = 1.0e-12 absolute -input_arguments=-tdy_disc bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu [steady-bdm-prob2] -input_arguments=-tdy_disc bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu [steady-bdm-prob3] -input_arguments=-tdy_disc bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu [steady-bdm-prob4] -input_arguments=-tdy_disc bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu diff --git a/src/tdycore.c b/src/tdycore.c index d03a2d55..ce52d245 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -720,7 +720,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Numerics options TDyDiscretization discretization; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); - ierr = PetscOptionsEnum("-tdy_disc","Discretization", + ierr = PetscOptionsEnum("-tdy_discretization","Discretization", "TDySetDiscretization",TDyDiscretizations, (PetscEnum)options->discretization,(PetscEnum *)&discretization, NULL); CHKERRQ(ierr); From 15c22de108ba8821b5c824c8cc5729f729a22671 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 12:55:11 -0700 Subject: [PATCH 13/23] Worked through the transient demo (and fixed a few bugs). --- demo/steady/steady.c | 9 +++-- demo/transient/transient.c | 68 +++++++++++++++++++++++--------------- src/tdycore.c | 13 ++++++-- 3 files changed, 57 insertions(+), 33 deletions(-) diff --git a/demo/steady/steady.c b/demo/steady/steady.c index c4a43096..fab4c1c4 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -558,7 +558,6 @@ PetscErrorCode GeometryColumn(DM dm){ typedef struct DMOptions { PetscInt dim; // Dimension of DM (2 or 3) PetscInt N; // Number of cells on a side - PetscBool perturb; // whether to perturb randomly (as opposed to smoothly) PetscBool exo; // whether to load a named exodus file PetscBool column; // column mesh? const char* exofile; // name of the exodus file to load @@ -584,9 +583,10 @@ PetscErrorCode CreateDM(MPI_Comm comm, void* context, DM* dm) { const PetscInt faces[3] = {Nx ,Ny ,Nz }; const PetscReal lower[3] = {0.0,0.0,0.0}; const PetscReal upper[3] = {Lx ,Ly ,Lz }; - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,options->dim,PETSC_FALSE, - faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); + ierr = DMPlexCreateBoxMesh(comm, options->dim, PETSC_FALSE, + faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); } + PetscFunctionReturn(0); } int main(int argc, char **argv) { @@ -632,8 +632,7 @@ int main(int argc, char **argv) { // Specify a special DM to be constructed for this demo, and pass it the // relevant options. - DMOptions dm_options = {.N = N, .dim = dim, .perturb = perturb, - .exo = exo, .exofile = exofile}; + DMOptions dm_options = {.N = N, .dim = dim, .exo = exo, .exofile = exofile}; ierr = TDySetDMConstructor(tdy, &dm_options, CreateDM); CHKERRQ(ierr); // Apply overrides. diff --git a/demo/transient/transient.c b/demo/transient/transient.c index ac70aa4c..b0c7febb 100644 --- a/demo/transient/transient.c +++ b/demo/transient/transient.c @@ -58,6 +58,34 @@ PetscErrorCode PerturbInteriorVertices(DM dm,PetscReal h) { PetscFunctionReturn(0); } +// This data is used by CreateDM below to create a DM for this demo. +typedef struct DMOptions { + PetscInt dim; // Dimension of DM (2 or 3) + PetscInt N; // Number of cells on a side + PetscBool exo; // whether to load a named exodus file + const char* exofile; // name of the exodus file to load +} DMOptions; + +// This function creates a DM specifically for this demo. Overrides are applied +// to the resulting DM with TDySetFromOptions. +PetscErrorCode CreateDM(MPI_Comm comm, void* context, DM* dm) { + int ierr; + DMOptions* options = context; + + PetscInt N = options->N; + if(options->exo){ + ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,options->exofile, + PETSC_TRUE,dm); CHKERRQ(ierr); + } else { + const PetscInt faces[3] = {N,N,N }; + const PetscReal lower[3] = {-0.5,-0.5,-0.5}; + const PetscReal upper[3] = {+0.5,+0.5,+0.5}; + ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,options->dim,PETSC_FALSE, + faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); + ierr = PerturbInteriorVertices(*dm,1./N); CHKERRQ(ierr); + } +} + int main(int argc, char **argv) { /* Initialize */ PetscErrorCode ierr; @@ -67,9 +95,11 @@ int main(int argc, char **argv) { PetscBool exo = PETSC_FALSE; ierr = TDyInit(argc, argv); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy; - ierr = TDyCreate(&tdy); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,WY); CHKERRQ(ierr); + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); + ierr = TDySetMode(tdy, RICHARDS); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,WY); CHKERRQ(ierr); ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL, "Transient Options",""); CHKERRQ(ierr); @@ -81,33 +111,19 @@ int main(int argc, char **argv) { exofile,exofile,256,&exo); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); - /* Create and distribute the mesh */ - DM dm, dmDist = NULL; - DMLabel marker; - if(exo){ - ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,exofile, - PETSC_TRUE,&dm); CHKERRQ(ierr); - //ierr = DMPlexOrient(dm); CHKERRQ(ierr); - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMCreateLabel(dm,"marker"); CHKERRQ(ierr); - ierr = DMGetLabel(dm,"marker",&marker); CHKERRQ(ierr); - ierr = DMPlexMarkBoundaryFaces(dm,1,marker); CHKERRQ(ierr); - }else{ - const PetscInt faces[3] = {N,N,N }; - const PetscReal lower[3] = {-0.5,-0.5,-0.5}; - const PetscReal upper[3] = {+0.5,+0.5,+0.5}; - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dim,PETSC_FALSE,faces,lower,upper, - NULL,PETSC_TRUE,&dm); CHKERRQ(ierr); - ierr = PerturbInteriorVertices(dm,1./N); CHKERRQ(ierr); - } - ierr = DMPlexDistribute(dm, 1, NULL, &dmDist); - if (dmDist) {DMDestroy(&dm); dm = dmDist;} - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + // Specify a special DM to be constructed for this demo, and pass it the + // relevant options. + DMOptions dm_options = {.N = N, .dim = dim, .exo = exo, .exofile = exofile}; + ierr = TDySetDMConstructor(tdy, &dm_options, CreateDM); CHKERRQ(ierr); - ierr = TDySetDM(tdy,dm); CHKERRQ(ierr); + // Apply overrides. ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); + // View the configured DM. + DM dm; + TDyGetDM(tdy, &dm); + ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + /* Setup problem parameters */ ierr = TDySetPorosity(tdy,Porosity); CHKERRQ(ierr); ierr = TDySetPermeabilityScalar(tdy,Permeability); CHKERRQ(ierr); diff --git a/src/tdycore.c b/src/tdycore.c index ce52d245..34374ab7 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -199,6 +199,8 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { PetscFunctionBegin; TDyOptions *options = &tdy->options; + options->mode = RICHARDS; + options->discretization = MPFA_O; options->gravity_constant = 9.8068; options->rho_type = WATER_DENSITY_CONSTANT; options->mu_type = WATER_VISCOSITY_CONSTANT; @@ -678,7 +680,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Model options - TDyMode mode; + TDyMode mode = options->mode; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Model options",""); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_mode","Flow mode", "TDySetMode",TDyModes,(PetscEnum)options->mode, @@ -718,7 +720,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options - TDyDiscretization discretization; + TDyDiscretization discretization = options->discretization; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_discretization","Discretization", "TDySetDiscretization",TDyDiscretizations, @@ -848,6 +850,13 @@ PetscErrorCode TDySetMode(TDy tdy, TDyMode mode) { PetscFunctionBegin; tdy->options.mode = mode; tdy->setup_flags |= TDyModeSet; + + // If we are resetting the mode and have already set the discretization, + // we call TDySetDiscretization again. + if (tdy->setup_flags & TDyDiscretizationSet) { + PetscInt ierr = TDySetDiscretization(tdy, tdy->options.discretization); CHKERRQ(ierr); + } + PetscFunctionReturn(0); } From f743f30d9f978e48ff7c8410bca518ba334dcc70 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 15:30:58 -0700 Subject: [PATCH 14/23] Worked through transient Fortran demo and fixed up some issues with DM constructors. --- demo/transient/transient_mpfaof90.F90 | 44 ++++++++++++++++++--------- include/finclude/tdycore.h | 2 +- src/f90-mod/tdycoremod.F90 | 16 +++++----- src/interface/ftn/tdycoref.c | 12 ++++++++ src/tdycore.c | 16 +++++++++- 5 files changed, 65 insertions(+), 25 deletions(-) diff --git a/demo/transient/transient_mpfaof90.F90 b/demo/transient/transient_mpfaof90.F90 index 8494adaa..cc88cc58 100644 --- a/demo/transient/transient_mpfaof90.F90 +++ b/demo/transient/transient_mpfaof90.F90 @@ -4,6 +4,9 @@ module mpfaof90mod #include #include + ! Module variables + PetscInt :: dim = 3 + contains subroutine PorosityFunction(tdy,x,theta,dummy,ierr) @@ -53,6 +56,27 @@ subroutine PressureFunction(tdy,x,pressure,dummy,ierr) ierr = 0 end subroutine PressureFunction + subroutine CreateDM(comm, dm, ierr) + use iso_c_binding + use petscdm + implicit none + + MPI_Comm :: comm + DM :: dm + PetscErrorCode :: ierr + + PetscInt :: faces(3), nx, ny, nz + PetscReal :: lower(3), upper(3) + + nx = 3; ny = 3; nz = 3; + faces(1) = nx; faces(2) = ny; faces(3) = nz; + lower(:) = 0.d0; + upper(:) = 1.d0; + + call DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, faces, lower, upper, & + PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); + CHKERRA(ierr) + end subroutine CreateDM end module mpfaof90mod program main @@ -76,10 +100,7 @@ program main TS :: ts PetscInt :: rank, successful_exit_code PetscBool :: flg - PetscInt :: dim, faces(3) - PetscReal :: lower(3), upper(3) PetscErrorCode :: ierr - PetscInt :: nx, ny, nz PetscInt, pointer :: index(:) PetscReal, pointer :: residualSat(:), blockPerm(:), liquid_sat(:), liquid_mass(:) PetscReal, pointer :: p_loc(:) @@ -90,11 +111,11 @@ program main CHKERRA(ierr); call TDyCreate(tdy, ierr); CHKERRA(ierr); - call TDySetDiscretizationMethod(tdy,MPFA_O,ierr); + call TDySetMode(tdy,RICHARDS,ierr); + CHKERRA(ierr); + call TDySetDiscretization(tdy,MPFA_O,ierr); CHKERRA(ierr); - nx = 3; ny = 3; nz = 3; - dim = 3 successful_exit_code= 0 call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-successful_exit_code',successful_exit_code,flg,ierr); @@ -102,19 +123,14 @@ program main call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr); CHKERRA(ierr) - faces(1) = nx; faces(2) = ny; faces(3) = nz; - lower(:) = 0.d0; - upper(:) = 1.d0; - - call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & - PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); - CHKERRA(ierr); - call TDySetDM(tdy, dm, ierr); + ! Set a constructor for a DM. + call TDySetDMConstructor(tdy, CreateDM,ierr); CHKERRA(ierr); call TDySetFromOptions(tdy,ierr); CHKERRA(ierr); + call TDyGetDM(tdy, dm, ierr); call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr); CHKERRA(ierr); diff --git a/include/finclude/tdycore.h b/include/finclude/tdycore.h index f4b01478..fdc9d1e7 100644 --- a/include/finclude/tdycore.h +++ b/include/finclude/tdycore.h @@ -2,7 +2,7 @@ #define TDYCORECOREDEF_H #define TDy type(tTDy) -#define TDyMethod PetscEnum +#define TDyDiscretization PetscEnum #define MPFAO_DIRICHLET_BC 0 #define MPFAO_NEUMANN_BC 1 diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index 3e35b1ad..98cfe8ac 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -453,12 +453,11 @@ function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) end interface abstract interface - subroutine TDyDMConstructor(comm, context, dm, ierr) + subroutine TDyDMConstructor(comm, dm, ierr) use, intrinsic :: iso_c_binding use tdycoredef use petscdm MPI_Comm :: comm - type(c_ptr), value :: context DM :: dm PetscErrorCode :: ierr end subroutine @@ -501,22 +500,21 @@ subroutine TDySetDMConstructor(tdy, dm_ctor, ierr) use, intrinsic :: iso_c_binding use tdycoredef implicit none - TDy, target :: tdy + TDy, target :: tdy procedure(TDyDMConstructor) :: dm_ctor PetscErrorCode :: ierr interface - function SetDMConstructor(tdy, context, dm_ctor) bind (c, name="TDySetDMConstructor") result(ierr) + function SetDMConstructor(tdy, dm_ctor) bind (c, name="TDySetDMConstructorF90") result(ierr) use, intrinsic :: iso_c_binding implicit none - type(c_ptr) :: tdy - type(c_ptr) :: context - type(c_funptr) :: dm_ctor - integer(c_int) :: ierr + type(c_ptr), value :: tdy + type(c_funptr), value :: dm_ctor + integer(c_int) :: ierr end function end interface - ierr = SetDMConstructor(c_loc(tdy), c_null_ptr, c_funloc(dm_ctor)) + ierr = SetDMConstructor(c_loc(tdy), c_funloc(dm_ctor)) end subroutine subroutine TDySelectPorosityFunction(tdy, name, ierr) diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 612a23a2..68104f07 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -18,6 +18,7 @@ #define tdydtimeintegratorsettimestep_ TDYTIMEINTEGRATORSETTIMESTEP #define tdydtimeintegratoroutputregression_ TDYTIMEINTEGRATOROUTPUTREGRESSION #define tdysetup_ TDYSETUP +#define tdygetdm_ TDYGETDM #define tdysetwaterdensitytype_ TDYSETWATERDENSITYTYPE #define tdysetmpfaogmatrixmethod_ TDYSETMPFAOGMATRIXMETHOD #define tdysetmpfaoboundaryconditiontype_ TDYSETMPFAOGBOUNDARYCONDITIONTYPE @@ -79,6 +80,7 @@ #define tdydtimeintegratorsettimestep_ tdydtimeintegratorsettimestep #define tdydtimeintegratoroutputregression_ tdydtimeintegratoroutputregression #define tdysetup_ tdysetup +#define tdygetdm_ tdygetdm #define tdysetwaterdensitytype_ tdysetwaterdensitytype #define tdysetmpfaogmatrixmethod_ tdysetmpfaogmatrixmethod #define tdysetmpfaoboundaryconditiontype_ tdysetmpfaoboundaryconditiontype @@ -250,6 +252,16 @@ PETSC_EXTERN void tdysetup_(TDy _tdy, int *__ierr){ } #endif +#if defined(__cplusplus) +extern "C" { +#endif +PETSC_EXTERN void tdygetdm_(TDy _tdy, DM *dm, int *__ierr){ +*__ierr = TDyGetDM((TDy)PetscToPointer((_tdy)), dm); +} +#if defined(__cplusplus) +} +#endif + #if defined(__cplusplus) extern "C" { #endif diff --git a/src/tdycore.c b/src/tdycore.c index 34374ab7..4082d35b 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -291,13 +291,23 @@ PetscErrorCode TDySetDMConstructor(TDy tdy, void* context, } if (tdy->setup_flags & TDyOptionsSet) { SETERRQ(comm,PETSC_ERR_USER, - "You must call TDyDefineDM before TDySetFromOptions()"); + "You must call TDySetDMConstructor before TDySetFromOptions()"); } tdy->create_dm_context = context; tdy->ops->create_dm = dm_func; PetscFunctionReturn(0); } +// This function is a wrapper used to eliminate the context pointer argument +// from TDySetDMConstructor so the function can be called from Fortran. +static void (*create_dm_f90_)(MPI_Fint*, DM*, PetscErrorCode*) = NULL; +PetscErrorCode TDySetDMConstructorF90(TDy tdy, + void (*dm_func)(MPI_Fint*, DM*, PetscErrorCode*)) { + PetscFunctionBegin; + create_dm_f90_ = dm_func; + PetscFunctionReturn(0); +} + PetscErrorCode TDyMalloc(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; @@ -402,6 +412,10 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { if (tdy->ops->create_dm) { // We've been instructed to create a DM ourselves. ierr = tdy->ops->create_dm(comm, tdy->create_dm_context, &dm); + } else if (create_dm_f90_) { + // Create a DM all-Fortran-like. + MPI_Fint comm_f = MPI_Comm_c2f(comm); + create_dm_f90_(&comm_f, &dm, &ierr); } else { ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); } From 2feb082dd3493b0187b61eeb5a2681174d84ffae Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 16:35:31 -0700 Subject: [PATCH 15/23] Converted transient SNES demo. --- demo/transient/transient_mpfaof90.F90 | 1 - demo/transient/transient_snes_mpfaof90.F90 | 113 +++++++++++---------- demo/transient/transient_snes_mpfaof90.cfg | 2 +- 3 files changed, 63 insertions(+), 53 deletions(-) diff --git a/demo/transient/transient_mpfaof90.F90 b/demo/transient/transient_mpfaof90.F90 index cc88cc58..343c9b99 100644 --- a/demo/transient/transient_mpfaof90.F90 +++ b/demo/transient/transient_mpfaof90.F90 @@ -57,7 +57,6 @@ subroutine PressureFunction(tdy,x,pressure,dummy,ierr) end subroutine PressureFunction subroutine CreateDM(comm, dm, ierr) - use iso_c_binding use petscdm implicit none diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90 index 14113b83..6d74f0f3 100644 --- a/demo/transient/transient_snes_mpfaof90.F90 +++ b/demo/transient/transient_snes_mpfaof90.F90 @@ -4,6 +4,13 @@ module snes_mpfaof90mod #include #include + ! DM-related variables for CreateDM. + character (len=256) :: mesh_filename + PetscBool :: mesh_file_flg + PetscInt :: nx, ny, nz + PetscInt :: dim + PetscInt :: dm_plex_extrude_layers + contains subroutine PorosityFunction(tdy,x,theta,dummy,ierr) @@ -113,6 +120,40 @@ subroutine PressureFunction(tdy,x,pressure,dummy,ierr) ierr = 0 end subroutine PressureFunction + subroutine CreateDM(comm, dm, ierr) + use petscdm + implicit none + + MPI_Comm :: comm + DM :: dm, edm + PetscErrorCode :: ierr + + PetscInt :: faces(3) + PetscReal :: lower(3), upper(3) + + if (.not.mesh_file_flg) then + faces(1) = nx; faces(2) = ny; faces(3) = nz; + lower(:) = 0.d0; + upper(:) = 1.d0; + + call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & + PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); + CHKERRA(ierr); + else + call DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh_filename, PETSC_TRUE, dm, ierr); + CHKERRA(ierr); + call DMGetDimension(dm, dim, ierr); + CHKERRA(ierr); + if (dm_plex_extrude_layers > 0) then + call DMPlexExtrude(dm, PETSC_DETERMINE, -1.d0, PETSC_TRUE, PETSC_NULL_REAL, PETSC_TRUE, edm, ierr); + CHKERRA(ierr); + call DMDestroy(dm ,ierr); + dm = edm + end if + endif + + end subroutine + end module snes_mpfaof90mod program main @@ -134,30 +175,31 @@ program main implicit none TDy :: tdy - DM :: dm, edm, dmDist + DM :: dm Vec :: U !TS :: ts SNES :: snes PetscInt :: rank, successful_exit_code PetscBool :: flg - PetscInt :: dim, faces(3) - PetscReal :: lower(3), upper(3) PetscErrorCode :: ierr - PetscInt :: nx, ny, nz, ncell, bc_type, dm_plex_extrude_layers + PetscInt :: ncell, bc_type PetscInt , pointer :: index(:) PetscReal , pointer :: residualSat(:), blockPerm(:), liquid_sat(:), liquid_mass(:) PetscReal , pointer :: alpha(:), m(:) PetscReal :: perm(9), resSat PetscInt :: c, cStart, cEnd, j, nvalues,g, max_steps, step PetscReal :: dtime, mass_pre, mass_post, ic_value - character (len=256) :: mesh_filename, ic_filename + character (len=256) :: ic_filename character(len=256) :: string, bc_type_name - PetscBool :: mesh_file_flg, ic_file_flg, pflotran_consistent, use_tdydriver + PetscBool :: ic_file_flg, pflotran_consistent, use_tdydriver PetscViewer :: viewer PetscInt :: step_mod PetscFE :: fe SNESConvergedReason :: reason + dim = 3 + nx = 1; ny = 1; nz = 15; + call TDyInit(ierr); CHKERRA(ierr); @@ -169,11 +211,11 @@ program main call TDyCreate(tdy, ierr); CHKERRA(ierr); - call TDySetDiscretizationMethod(tdy,MPFA_O,ierr); + call TDySetMode(tdy,RICHARDS,ierr); + CHKERRA(ierr); + call TDySetDiscretization(tdy,MPFA_O,ierr); CHKERRA(ierr); - nx = 1; ny = 1; nz = 15; - dim = 3 successful_exit_code= 0 max_steps = 2 dtime = 1800.d0 @@ -224,27 +266,17 @@ program main bc_type = MPFAO_SEEPAGE_BC endif - if (.not.mesh_file_flg) then - faces(1) = nx; faces(2) = ny; faces(3) = nz; - lower(:) = 0.d0; - upper(:) = 1.d0; + ! Set a constructor for a DM. + call TDySetDMConstructor(tdy, CreateDM,ierr); + CHKERRA(ierr); - call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & - PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); - CHKERRA(ierr); - else - call DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh_filename, PETSC_TRUE, dm, ierr); - CHKERRA(ierr); - call DMGetDimension(dm, dim, ierr); - CHKERRA(ierr); - if (dm_plex_extrude_layers > 0) then - call DMPlexExtrude(dm, PETSC_DETERMINE, -1.d0, PETSC_TRUE, PETSC_NULL_REAL, PETSC_TRUE, edm, ierr); - CHKERRA(ierr); - call DMDestroy(dm ,ierr); - dm = edm - end if - endif + ! Apply overrides. + call TDySetFromOptions(tdy,ierr); + CHKERRA(ierr); + ! Set up the discretization. + call TDyGetDM(tdy, dm, ierr); + CHKERRA(ierr); call DMGetDimension(dm, dim, ierr); CHKERRA(ierr) call PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, "p_", -1, fe, ierr); @@ -259,25 +291,6 @@ program main CHKERRA(ierr) call DMSetUseNatural(dm, PETSC_TRUE, ierr); CHKERRA(ierr) - - call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr); - CHKERRA(ierr); - if (dmDist /= PETSC_NULL_DM) then - call DMDestroy(dm, ierr); - CHKERRA(ierr); - dm = dmDist; - end if - call DMSetUp(dm,ierr); - CHKERRA(ierr) - - - call DMSetFromOptions(dm, ierr); - CHKERRA(ierr); - - - call TDySetWaterDensityType(tdy,WATER_DENSITY_EXPONENTIAL,ierr); - CHKERRA(ierr) - call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr); CHKERRA(ierr); @@ -306,10 +319,8 @@ program main enddo enddo - call TDySetDM(tdy, dm, ierr); - CHKERRA(ierr); - call TDySetFromOptions(tdy,ierr); - CHKERRA(ierr); + call TDySetWaterDensityType(tdy,WATER_DENSITY_EXPONENTIAL,ierr); + CHKERRA(ierr) if (pflotran_consistent) then ! call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr); diff --git a/demo/transient/transient_snes_mpfaof90.cfg b/demo/transient/transient_snes_mpfaof90.cfg index 4339405f..424d79d8 100644 --- a/demo/transient/transient_snes_mpfaof90.cfg +++ b/demo/transient/transient_snes_mpfaof90.cfg @@ -26,7 +26,7 @@ input_arguments=-max_steps 48 -ic_value 90325 -nx 5 -ny 5 -nz 5 -tdy_water_densi input_arguments=-max_steps 10 -use_tdydriver -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-tdydriver -pc_type lu [transient-snes-mpfaof90-transect-neumann] -input_arguments=-mesh_filename ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -ic_value 57180 -dtime 1800 -tdy_mpfao_boundary_condition_type MPFAO_NEUMANN_BC -max_steps 10 -snes_rtol 1e-30 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-transect-neumann +input_arguments=-tdy_read_mesh ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -ic_value 57180 -dtime 1800 -tdy_mpfao_boundary_condition_type MPFAO_NEUMANN_BC -max_steps 10 -snes_rtol 1e-30 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-transect-neumann [transient-snes-mpfaof90-region] input_arguments=-max_steps 10 -tdy_regression_test -tdy_regression_test_num_cells_per_process 10 -tdy_regression_test_filename transient-snes-mpfaof90-region -pc_type lu -tdy_connected_region ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh_regions.bin -mesh_filename ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -tdy_mpfao_gmatrix_method MPFAO_GMATRIX_TPF -tdy_mpfao_boundary_condition_type MPFAO_DIRICHLET_BC -tdy_water_density EXPONENTIAL -pflotran_consistent From bb5ef4630f17e78cc6a061e6833b4f617925ad20 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Fri, 15 Oct 2021 16:59:32 -0700 Subject: [PATCH 16/23] Removed TPF. --- demo/steady/steady.cfg | 16 -- include/private/tdytpfimpl.h | 9 -- include/tdycore.h | 7 - src/tdycore.c | 31 +--- src/tdydriver.c | 1 - src/tdytpf.c | 280 ----------------------------------- 6 files changed, 3 insertions(+), 341 deletions(-) delete mode 100644 include/private/tdytpfimpl.h delete mode 100644 src/tdytpf.c diff --git a/demo/steady/steady.cfg b/demo/steady/steady.cfg index 6db2cc9f..a463e9df 100644 --- a/demo/steady/steady.cfg +++ b/demo/steady/steady.cfg @@ -7,10 +7,6 @@ standard= steady-wy-prob1-3D steady-wy-prob2-3D steady-wy-prob3-3D - steady-tpf-prob1 - steady-tpf-prob2 - steady-tpf-prob3 - steady-tpf-prob4 steady-bdm-prob1 steady-bdm-prob2 steady-bdm-prob3 @@ -40,18 +36,6 @@ input_arguments=-tdy_discretization wy -tdy_regression_test -problem 2 -dim 3 -N [steady-wy-prob3-3D] input_arguments=-tdy_discretization wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D -[steady-tpf-prob1] -input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 - -[steady-tpf-prob2] -input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 - -[steady-tpf-prob3] -input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 - -[steady-tpf-prob4] -input_arguments=-tdy_discretization tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 - [steady-bdm-prob1] pressure = 1.0e-12 absolute input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu diff --git a/include/private/tdytpfimpl.h b/include/private/tdytpfimpl.h deleted file mode 100644 index c0699474..00000000 --- a/include/private/tdytpfimpl.h +++ /dev/null @@ -1,9 +0,0 @@ -#if !defined(TDYTPFIMPL_H) -#define TDYTPFIMPL_H - -#include -#include - -PETSC_INTERN PetscErrorCode TDyTPF_Setup(TDy, DM); - -#endif diff --git a/include/tdycore.h b/include/tdycore.h index 4c43d6ed..8b10f391 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -19,8 +19,6 @@ PETSC_EXTERN const char *const TDyModes[]; /// This type enumerates discretizations supported by the dycore. typedef enum { - /// two-point flux, classic finite volumes - TPF=0, /// multi-point flux approximation - O method MPFA_O, /// multi-point flux approximation - O method using DAE @@ -186,11 +184,6 @@ PETSC_EXTERN PetscErrorCode TDyPostSolveSNESSolver(TDy,Vec); PETSC_EXTERN PetscErrorCode TDyOutputRegression(TDy,Vec); -PETSC_EXTERN PetscErrorCode TDyTPFComputeSystem(TDy,Mat,Vec); -PETSC_EXTERN PetscReal TDyTPFPressureNorm(TDy,Vec); -PETSC_EXTERN PetscReal TDyTPFVelocityNorm(TDy,Vec); -PETSC_EXTERN PetscErrorCode TDyTPFCheckMeshSuitability(TDy); - PETSC_EXTERN PetscErrorCode TDyWYComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscErrorCode TDyBDMComputeSystem(TDy,Mat,Vec); diff --git a/src/tdycore.c b/src/tdycore.c index 4082d35b..1cfefefb 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -9,7 +9,6 @@ static char help[] = "TDycore \n\ #include #include #include -#include #include #include #include @@ -25,7 +24,6 @@ static char help[] = "TDycore \n\ #include const char *const TDyDiscretizations[] = { - "TPF", "MPFA_O", "MPFA_O_DAE", "MPFA_O_TRANSIENTVAR", @@ -892,10 +890,7 @@ PetscErrorCode TDySetDiscretization(TDy tdy, TDyDiscretization discretization) { // Set function pointers for operations. if (tdy->options.mode == RICHARDS) { - if (discretization == TPF) { - tdy->ops->set_from_options = NULL; // FIXME: define this and move things here! - tdy->ops->setup = TDyTPF_Setup; - } else if (discretization == MPFA_O) { + if (discretization == MPFA_O) { tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; tdy->ops->setup = TDyRichards_MPFA_O_Setup; } else if (discretization == MPFA_O_DAE) { @@ -910,6 +905,8 @@ PetscErrorCode TDySetDiscretization(TDy tdy, TDyDiscretization discretization) { } else if (discretization == WY) { tdy->ops->set_from_options = NULL; tdy->ops->setup = TDyWY_Setup; + } else { + SETERRQ(comm,PETSC_ERR_USER, "Invalid discretization given!"); } } else if (tdy->options.mode == TH) { if (discretization == MPFA_O) { @@ -976,9 +973,6 @@ PetscErrorCode TDySetIFunction(TS ts,TDy tdy) { ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"IFunction not implemented for TPF"); - break; case MPFA_O: switch (tdy->options.mode) { case RICHARDS: @@ -1016,9 +1010,6 @@ PetscErrorCode TDySetIJacobian(TS ts,TDy tdy) { TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)ts,&comm); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"IJacobian not implemented for TPF"); - break; case MPFA_O: ierr = TDyCreateJacobian(tdy); CHKERRQ(ierr); switch (tdy->options.mode) { @@ -1064,9 +1055,6 @@ PetscErrorCode TDySetSNESFunction(SNES snes,TDy tdy) { ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"SNESFunction not implemented for TPF"); - break; case MPFA_O: ierr = SNESSetFunction(snes,tdy->residual,TDyMPFAOSNESFunction,tdy); CHKERRQ(ierr); break; @@ -1101,9 +1089,6 @@ PetscErrorCode TDySetSNESJacobian(SNES snes,TDy tdy) { ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"SNESJacobian not implemented for TPF"); - break; case MPFA_O: ierr = SNESSetJacobian(snes,tdy->J,tdy->J,TDyMPFAOSNESJacobian,tdy); CHKERRQ(ierr); break; @@ -1131,9 +1116,6 @@ PetscErrorCode TDyComputeSystem(TDy tdy,Mat K,Vec F) { TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - ierr = TDyTPFComputeSystem(tdy,K,F); CHKERRQ(ierr); - break; case MPFA_O: ierr = TDyMPFAOComputeSystem(tdy,K,F); CHKERRQ(ierr); break; @@ -1536,10 +1518,6 @@ PetscErrorCode TDyComputeErrorNorms(TDy tdy,Vec U,PetscReal *normp, TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - if(normp != NULL) { *normp = TDyTPFPressureNorm(tdy,U); } - if(normv != NULL) { *normv = TDyTPFVelocityNorm(tdy,U); } - break; case MPFA_O: if(normv) { ierr = TDyMPFAORecoverVelocity(tdy,U); CHKERRQ(ierr); @@ -1627,9 +1605,6 @@ PetscErrorCode TDyPreSolveSNESSolver(TDy tdy) { ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); switch (tdy->options.discretization) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"TDyPreSolveSNESSolver not implemented for TPF"); - break; case MPFA_O: ierr = TDyMPFAOSNESPreSolve(tdy); CHKERRQ(ierr); break; diff --git a/src/tdydriver.c b/src/tdydriver.c index d025a848..47a0d584 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -23,7 +23,6 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { } switch(tdy->options.discretization) { - case TPF: case MPFA_O: break; case MPFA_O_DAE: diff --git a/src/tdytpf.c b/src/tdytpf.c deleted file mode 100644 index 964e1348..00000000 --- a/src/tdytpf.c +++ /dev/null @@ -1,280 +0,0 @@ -#include -#include - -PETSC_STATIC_INLINE void Waxpy(PetscInt dim, PetscScalar a, - const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];} -PETSC_STATIC_INLINE PetscScalar Dot(PetscInt dim, const PetscScalar *x, - const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;} -PETSC_STATIC_INLINE PetscReal Norm(PetscInt dim, const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(Dot(dim,x,x)));} - -PetscErrorCode TDyTPF_Setup(TDy tdy, DM dm) { - PetscErrorCode ierr; - MPI_Comm comm; - PetscInt dim,c,cStart,cEnd,pStart,pEnd; - PetscSection sec; - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - - ierr = PetscObjectGetComm((PetscObject)dm,&comm); CHKERRQ(ierr); - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - ierr = DMPlexGetChart(dm,&pStart,&pEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - - /* Setup the section, 1 dof per cell */ - ierr = PetscSectionCreate(comm,&sec); CHKERRQ(ierr); - ierr = PetscSectionSetNumFields(sec,1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec,0,"Pressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec,0,1); CHKERRQ(ierr); - ierr = DMPlexGetChart(dm,&pStart,&pEnd); CHKERRQ(ierr); - ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); - for(c=cStart; cdm; - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - if(!tdy->options.tpf_allow_all_meshes) { - ierr = TDyTPFCheckMeshSuitability(tdy); CHKERRQ(ierr); - } - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - dim2 = dim*dim; - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - MaterialProp *matprop = tdy->matprop; - - // Here's our domain boundary label. - DMLabel bnd_label; - ierr = DMGetLabel(tdy->dm, "boundary", &bnd_label); CHKERRQ(ierr); - - // Loop over all faces and compute fluxes. - for (PetscInt f = fStart; f < fEnd; ++f) { - - const PetscInt *supp; - ierr = DMPlexGetSupport(dm,f,&supp); CHKERRQ(ierr); - - // wrt - // v = -Ki (pd-p0)/dist - // 00 Ki/dist - // 0 -Ki pd / dist - - 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); - dist = Norm(dim,pnt2pnt); - Ki = 0.5*(matprop->K[(supp[0]-cStart)*dim2]+matprop->K[(supp[1]-cStart)*dim2]); - - //wrt 0 - // v = -Ki (p1-p0)/dist - // 00 Ki/dist - // 01 -Ki/dist - - //wrt 1 - // v = -Ki (p0-p1)/dist - // 11 Ki/dist - // 10 -Ki/dist - - ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); - ierr = DMPlexGetPointGlobal(dm,supp[1],&col,&junk); CHKERRQ(ierr); - ierr = MatSetValue(K,row,row, Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,row,col,-Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,col,col, Ki/dist*tdy->V[f]/tdy->V[supp[1]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,col,row,-Ki/dist*tdy->V[f]/tdy->V[supp[1]],ADD_VALUES); - CHKERRQ(ierr); - } - - if(tdy->ops->computeforcing) { - for(c=cStart; cops->computeforcing)(tdy,&(tdy->X[c*dim]),&force,tdy->forcingctx);CHKERRQ(ierr);; - ierr = VecSetValue(F,row,force,ADD_VALUES); CHKERRQ(ierr); - } - } - - // max = 3.63988 - // min = 3.15550 - - ierr = VecAssemblyBegin(F); CHKERRQ(ierr); - ierr = VecAssemblyEnd (F); CHKERRQ(ierr); - ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - ierr = MatAssemblyEnd (K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(0); -} - -/* - norm = sqrt( sum_i^n V_i * ( p(X_i) - P_i )^2 ) - - where n is the number of cells. - */ -PetscReal TDyTPFPressureNorm(TDy tdy,Vec U) { - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - PetscErrorCode ierr; - PetscSection sec; - PetscInt c,cStart,cEnd,offset,dim,gref,junk; - PetscReal p,*u,norm,norm_sum; - DM dm = tdy->dm; - if(!(tdy->ops->compute_boundary_pressure)) { - SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_USER, - "You must set the pressure function with TDySetBoundaryPressureFn"); - } - norm = 0; - ierr = VecGetArray(U,&u); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - ierr = DMGetSection(dm,&sec); CHKERRQ(ierr); - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - for(c=cStart; cops->compute_boundary_pressure)(tdy,&(tdy->X[c*dim]),&p,tdy->boundary_pressure_ctx);CHKERRQ(ierr); - norm += tdy->V[c]*PetscSqr(u[offset]-p); - } - ierr = MPI_Allreduce(&norm,&norm_sum,1,MPIU_REAL,MPI_SUM, - PetscObjectComm((PetscObject)U)); CHKERRQ(ierr); - norm_sum = PetscSqrtReal(norm_sum); - ierr = VecRestoreArray(U,&u); CHKERRQ(ierr); - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(norm_sum); -} - -PetscReal TDyTPFVelocityNorm(TDy tdy,Vec U) { - PetscErrorCode ierr; - 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; - TDY_START_FUNCTION_TIMER() - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - dim2 = dim*dim; - ierr = VecGetArray(U,&u); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - face_error = 0; - norm = 0; - MaterialProp *matprop = tdy->matprop; - - // Here's our domain boundary label. - DMLabel bnd_label; - ierr = DMGetLabel(tdy->dm, "boundary", &bnd_label); CHKERRQ(ierr); - - for(PetscInt c=cStart; cops->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); - Ki = matprop->K[(supp[0]-cStart)*dim2]; - ierr = (*tdy->ops->compute_boundary_pressure)(tdy, &(tdy->X[f*dim]),&p, tdy->boundary_pressure_ctx);CHKERRQ(ierr); - sign = PetscSign(TDyADotBMinusC(&(tdy->N[dim*f]),&(tdy->X[dim*f]), - &(tdy->X[dim*supp[0]]),dim)); - va = sign* -Ki*(p-u[(supp[0]-cStart)])/dist; - ierr = (*tdy->ops->compute_boundary_velocity)(tdy,&(tdy->X[f*dim]),&(vel[0]),tdy->boundary_velocity_ctx);CHKERRQ(ierr); - ve = TDyADotB(vel,&(tdy->N[dim*f]),dim); - face_error = PetscSqr((va-ve)/tdy->V[f]); - } else { - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[supp[1]*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = 0.5*(matprop->K[(supp[0]-cStart)*dim2]+matprop->K[(supp[1]-cStart)*dim2]); - sign = PetscSign(TDyADotBMinusC(&(tdy->N[dim*f]),&(tdy->X[dim*f]), - &(tdy->X[dim*c]),dim)); - va = sign* -Ki*(u[(supp[1]-cStart)]-u[(supp[0]-cStart)])/dist *tdy->V[f]; - if(c==supp[1]) va *= -1; - ierr = (*tdy->ops->compute_boundary_velocity)(tdy,&(tdy->X[f*dim]),&(vel[0]),tdy->boundary_velocity_ctx);CHKERRQ(ierr); - ve = TDyADotB(vel,&(tdy->N[dim*f]),dim)*tdy->V[f]; - face_error = PetscSqr((va-ve)/tdy->V[f]); - } - norm += face_error*tdy->V[c]; - } - } - ierr = MPI_Allreduce(&norm,&norm_sum,1,MPIU_REAL,MPI_SUM, - PetscObjectComm((PetscObject)dm)); CHKERRQ(ierr); - norm_sum = PetscSqrtReal(norm_sum); - ierr = VecRestoreArray(U,&u); CHKERRQ(ierr); - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(norm_sum); -} - -PetscErrorCode TDyTPFCheckMeshSuitability(TDy tdy) { - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - PetscErrorCode ierr; - PetscInt dim,f,fStart,fEnd; - PetscReal diff,dist,pnt2pnt[3]; - DM dm = tdy->dm; - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - for(f=fStart; fX[supp[1]*dim]),&(tdy->X[supp[0]*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - diff = PetscAbsReal(TDyADotB(&(tdy->N[f*dim]),pnt2pnt,dim)/dist); - if(PetscAbsReal(diff-1) > 10*PETSC_MACHINE_EPSILON) { - SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP, - "Mesh is unsuitable for a consistent two point flux method. To force rerun with -tdy_tpf_allow_unsuitable_mesh"); - } - } - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(0); -} From 904d8724b9870db054fc4a2041c23b888a3f0f8f Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Sat, 16 Oct 2021 12:49:27 -0700 Subject: [PATCH 17/23] Fixing minor glitches and cleaning things up. --- include/tdycore.h | 2 +- src/f90-mod/tdycoremod.F90 | 4 ++-- src/mpfao/tdympfao.c | 11 ----------- src/tdycore.c | 7 +------ 4 files changed, 4 insertions(+), 20 deletions(-) diff --git a/include/tdycore.h b/include/tdycore.h index 8b10f391..39dc62b9 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -20,7 +20,7 @@ PETSC_EXTERN const char *const TDyModes[]; /// This type enumerates discretizations supported by the dycore. typedef enum { /// multi-point flux approximation - O method - MPFA_O, + MPFA_O=0, /// multi-point flux approximation - O method using DAE MPFA_O_DAE, /// multipoint flux approximation - O method using TS transient (conservative) diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index 98cfe8ac..6a9df75b 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -31,12 +31,12 @@ subroutine TDySetMode(a,b,z) end subroutine TDySetMode end interface interface - subroutine TDySetDiscretizationMethod(a,b,z) + subroutine TDySetDiscretization(a,b,z) use tdycoredef TDy a PetscInt b integer z - end subroutine TDySetDiscretizationMethod + end subroutine TDySetDiscretization end interface interface subroutine TDySetFromOptions(a,z) diff --git a/src/mpfao/tdympfao.c b/src/mpfao/tdympfao.c index 06d4a20f..6d81c2bb 100644 --- a/src/mpfao/tdympfao.c +++ b/src/mpfao/tdympfao.c @@ -317,17 +317,6 @@ PetscErrorCode TDyMPFAO_AllocateMemoryForEnergySourceSinkValues(TDy tdy) { // repeated stuff when we move the data out of the TDy struct and into // discretization-specific types. -// Setup function for Richards + TPF -PetscErrorCode TDyRichards_TPF_Setup(TDy tdy, DM dm) { - PetscErrorCode ierr; - PetscFunctionBegin; - - MPI_Comm comm; - ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - // Setup function for Richards + MPFA_O PetscErrorCode TDyRichards_MPFA_O_Setup(TDy tdy, DM dm) { PetscErrorCode ierr; diff --git a/src/tdycore.c b/src/tdycore.c index 1cfefefb..9aa1acd8 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -51,7 +51,6 @@ const char *const TDyMPFAOBoundaryConditionTypes[] = { const char *const TDyModes[] = { "RICHARDS", "TH", - "UNSPECIFIED_MODE", /* */ "TDyMode","TDY_MODE_",NULL }; @@ -534,7 +533,7 @@ PetscErrorCode TDyGetDimension(TDy tdy,PetscInt *dim) { } /// Retrieves the DM used by the dycore. This must be called after -/// TDySetDM or TDySetFromOptions. +/// TDySetFromOptions. /// @param dm A pointer that stores the DM in use by the dycore PetscErrorCode TDyGetDM(TDy tdy,DM *dm) { PetscFunctionBegin; @@ -760,10 +759,6 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Mesh-related options ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_mesh", options->mesh_file,sizeof(options->mesh_file),&options->read_mesh); CHKERRQ(ierr); - if ((tdy->dm != NULL) && options->read_mesh) { - SETERRQ(comm,PETSC_ERR_USER, - "TDySetDM was called before TDySetFromOptions: can't read a mesh"); - } ierr = PetscOptionsBool("-tdy_output_mesh","Enable output of mesh attributes","",options->output_mesh,&(options->output_mesh),NULL); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_output_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->output_geom_attributes); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->read_geom_attributes); CHKERRQ(ierr); From 2791c3deceec7d2ea425cb81823794ed88cf8903 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 09:19:10 -0700 Subject: [PATCH 18/23] Fixed a bug caused by removing TPF discretization. --- src/f90-mod/tdycore.h | 5 ++--- src/interface/ftn/tdycoref.c | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/f90-mod/tdycore.h b/src/f90-mod/tdycore.h index 46f62e7c..9fc37a81 100644 --- a/src/f90-mod/tdycore.h +++ b/src/f90-mod/tdycore.h @@ -10,7 +10,6 @@ ! ! Types of TDy methods ! - PetscEnum TPF PetscEnum MPFA_O PetscEnum MPFA_O_DAE PetscEnum MPFA_O_TRANSIENT_VAR @@ -19,8 +18,8 @@ ! The parameters values need to match those defined in ! the C code (i.e. include/tdycore.h) - parameter (TPF=0,MPFA_O=1,MPFA_O_DAE=2,MPFA_O_TRANSIENT_VAR=3) - parameter (BDM=4,WY=5) + parameter (MPFA_O=0,MPFA_O_DAE=1,MPFA_O_TRANSIENT_VAR=2) + parameter (BDM=3,WY=4) ! ! Types of TDy modes diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 68104f07..b250b3a1 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -185,8 +185,8 @@ PETSC_EXTERN void tdysetmode_(TDy _tdy, PetscInt *mode, int *__ierr){ #if defined(__cplusplus) extern "C" { #endif -PETSC_EXTERN void tdysetdiscretization_(TDy _tdy, PetscInt *method, int *__ierr){ -*__ierr = TDySetDiscretization((TDy)PetscToPointer((_tdy)), *method); +PETSC_EXTERN void tdysetdiscretization_(TDy _tdy, PetscInt *discretization, int *__ierr){ +*__ierr = TDySetDiscretization((TDy)PetscToPointer((_tdy)), *discretization); } #if defined(__cplusplus) } From 7a6e9fe68f849eb04a5de911398c0f47236643c7 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 09:24:48 -0700 Subject: [PATCH 19/23] Updated richards and th drivers. --- demo/richards/richards_driver.c | 19 ++++++++++--------- demo/th/th_driver.c | 23 ++++++++++++----------- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/demo/richards/richards_driver.c b/demo/richards/richards_driver.c index 35c2fa98..96e9f9f2 100644 --- a/demo/richards/richards_driver.c +++ b/demo/richards/richards_driver.c @@ -7,18 +7,19 @@ int main(int argc, char **argv) { PetscErrorCode ierr; PetscInt successful_exit_code = 0; PetscMPIInt rank, size; + MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy = PETSC_NULL; TDyIOFormat format = HDF5Format; ierr = TDyInit(argc, argv); CHKERRQ(ierr); - ierr = TDyCreate(&tdy); CHKERRQ(ierr); + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,RICHARDS); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); - ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr); - ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"Beginning Richards Driver simulation.\n"); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); + ierr = MPI_Comm_rank(comm,&rank); CHKERRQ(ierr); + ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr); + PetscPrintf(comm,"Beginning Richards Driver simulation.\n"); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code", "Code passed on successful completion","", @@ -33,7 +34,7 @@ int main(int argc, char **argv) { if (!rank) { ierr = TDyIOSetIOProcess(tdy->io, PETSC_TRUE); CHKERRQ(ierr); } - PetscPrintf(PETSC_COMM_WORLD,"--\n"); + PetscPrintf(comm,"--\n"); if (size == 1) { ierr = TDyIOSetMode(tdy,format);CHKERRQ(ierr); ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr); @@ -44,8 +45,8 @@ int main(int argc, char **argv) { ierr = TDyOutputRegression(tdy,tdy->solution); CHKERRQ(ierr); ierr = TDyDestroy(&tdy); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"--\n"); - PetscPrintf(PETSC_COMM_WORLD,"Simulation complete.\n"); + PetscPrintf(comm,"--\n"); + PetscPrintf(comm,"Simulation complete.\n"); ierr = TDyFinalize(); CHKERRQ(ierr); return(successful_exit_code); } diff --git a/demo/th/th_driver.c b/demo/th/th_driver.c index 239232b0..74459aec 100644 --- a/demo/th/th_driver.c +++ b/demo/th/th_driver.c @@ -7,18 +7,19 @@ int main(int argc, char **argv) { PetscErrorCode ierr; PetscInt successful_exit_code; PetscMPIInt rank, size; + MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy = PETSC_NULL; - TDyIOFormat format = PetscViewerASCIIFormat; + TDyIOFormat format = PetscViewerASCIIFormat; ierr = TDyInit(argc, argv); CHKERRQ(ierr); - ierr = TDyCreate(&tdy); CHKERRQ(ierr); + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,TH); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); - ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr); - ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"Beginning TH Driver simulation.\n"); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); + ierr = MPI_Comm_rank(comm,&rank); CHKERRQ(ierr); + ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr); + PetscPrintf(comm,"Beginning TH Driver simulation.\n"); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code", "Code passed on successful completion","", @@ -33,19 +34,19 @@ int main(int argc, char **argv) { if (!rank) { ierr = TDyIOSetIOProcess(tdy->io, PETSC_TRUE); CHKERRQ(ierr); } - PetscPrintf(PETSC_COMM_WORLD,"--\n"); + PetscPrintf(comm,"--\n"); if (size == 1) { ierr = TDyIOSetMode(tdy,format);CHKERRQ(ierr); ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr); } - ierr = TDyTimeIntegratorRunToTime(tdy,tdy->ti->final_time); + ierr = TDyTimeIntegratorRunToTime(tdy,tdy->ti->final_time); CHKERRQ(ierr); if (size == 1) {ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr);} ierr = TDyOutputRegression(tdy,tdy->solution); CHKERRQ(ierr); ierr = TDyDestroy(&tdy); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"--\n"); - PetscPrintf(PETSC_COMM_WORLD,"Simulation complete.\n"); + PetscPrintf(comm,"--\n"); + PetscPrintf(comm,"Simulation complete.\n"); ierr = TDyFinalize(); CHKERRQ(ierr); return(successful_exit_code); } From 18bd8beec3e3753039b2985f4a8213227d45a0d5 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 09:55:26 -0700 Subject: [PATCH 20/23] Fixed a couple of bugs in TH and Richards demos. --- demo/richards/richards.cfg | 16 ++++++++-------- demo/richards/richards_driver.c | 2 +- demo/th/th.cfg | 8 ++++---- demo/th/th_driver.c | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/demo/richards/richards.cfg b/demo/richards/richards.cfg index d251ddb6..46756e7f 100644 --- a/demo/richards/richards.cfg +++ b/demo/richards/richards.cfg @@ -15,25 +15,25 @@ standard_parallel= pressure = 1.0e-12 relative [richards-driver-snes-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES [richards-driver-ts-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS [richards-driver-snes-prob1-np4] np=4 -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES [richards-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS [richards-driver-snes-initial-cond-isotropic-k] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-isotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_permeability_file initial.h5 -ic_file initial.h5 +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-isotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_permeability_file initial.h5 -ic_file initial.h5 [richards-driver-snes-initial-cond-anisotropic-k] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-anisotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_porosity_file initial.h5 -init_permeability_file initial.h5 -anisotropic_perm -ic_file initial.h5 -ic_dataset fields/IC +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-anisotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_porosity_file initial.h5 -init_permeability_file initial.h5 -anisotropic_perm -ic_file initial.h5 -ic_dataset fields/IC [richards-driver-64xy-3z-wedge] input_arguments=-problem 4 -dm_plex_simplex 0 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-64xy-3z-wedge -dm_plex_filename ../../share/meshes/64xy_3z_wedge.exo @@ -48,7 +48,7 @@ input_arguments=-snes_linesearch_basic -ts_dt 100 -ts_type beuler -ts_max_time 1 input_arguments=-snes_linesearch_type basic -ts_dt 100 -ts_type bdf -ts_adapt_type none -ts_max_time 1000 -ts_max_steps 10 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-transientvar -tdy_method MPFA_O_TRANSIENTVAR -pc_type lu [richards-driver-snes-checkpoint-write] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-write -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -enable_checkpoint +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-write -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -enable_checkpoint [richards-driver-snes-checkpoint-read] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-read -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -ic_file 3.15360e+03_chk.h5 -ic_dataset fields/IC \ No newline at end of file +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-read -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -ic_file 3.15360e+03_chk.h5 -ic_dataset fields/IC diff --git a/demo/richards/richards_driver.c b/demo/richards/richards_driver.c index 96e9f9f2..67efcbc3 100644 --- a/demo/richards/richards_driver.c +++ b/demo/richards/richards_driver.c @@ -7,11 +7,11 @@ int main(int argc, char **argv) { PetscErrorCode ierr; PetscInt successful_exit_code = 0; PetscMPIInt rank, size; - MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy = PETSC_NULL; TDyIOFormat format = HDF5Format; ierr = TDyInit(argc, argv); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,RICHARDS); CHKERRQ(ierr); ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); diff --git a/demo/th/th.cfg b/demo/th/th.cfg index de037b6a..98c3f2a0 100644 --- a/demo/th/th.cfg +++ b/demo/th/th.cfg @@ -9,17 +9,17 @@ standard_parallel= pressure = 1.0e-12 relative #[th-driver-snes-prob1] -#input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method SNES +#input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method SNES [th-driver-ts-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS #[th-driver-snes-prob1-np4] #np=4 -#input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES +#input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES [th-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS diff --git a/demo/th/th_driver.c b/demo/th/th_driver.c index 74459aec..5b4b271d 100644 --- a/demo/th/th_driver.c +++ b/demo/th/th_driver.c @@ -7,11 +7,11 @@ int main(int argc, char **argv) { PetscErrorCode ierr; PetscInt successful_exit_code; PetscMPIInt rank, size; - MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy = PETSC_NULL; TDyIOFormat format = PetscViewerASCIIFormat; ierr = TDyInit(argc, argv); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,TH); CHKERRQ(ierr); ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); From fcca22b0f697c41dd5a2174862eb999f4142c49b Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 11:12:06 -0700 Subject: [PATCH 21/23] Using a proper communicator in TDyDriver. --- src/tdydriver.c | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/tdydriver.c b/src/tdydriver.c index 47a0d584..036205cf 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -16,10 +16,17 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { SNES snes; SNESLineSearch linesearch; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + DM dm; + ierr = TDyGetDM(tdy,&dm); CHKERRQ(ierr); + PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); + if (dim != 3) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver currently only supports 3D"); + SETERRQ(comm,PETSC_ERR_USER,"Driver currently only supports 3D"); } switch(tdy->options.discretization) { @@ -29,7 +36,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case MPFA_O_TRANSIENTVAR: case BDM: case WY: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver not supported for specified discretization."); + SETERRQ(comm,PETSC_ERR_USER,"Driver not supported for specified discretization."); break; } @@ -81,7 +88,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case TH: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unrecognized flow mode."); + SETERRQ(comm,PETSC_ERR_USER,"Unrecognized flow mode."); } // check for unsupported time integration methods switch(tdy->ti->time_integration_method) { @@ -90,20 +97,20 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case RICHARDS: break; case TH: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"SNES not supported for TH mode."); + SETERRQ(comm,PETSC_ERR_USER,"SNES not supported for TH mode."); break; } break; case TDyTS: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unrecognized time integration method."); + SETERRQ(comm,PETSC_ERR_USER,"Unrecognized time integration method."); } // create time integrator switch(tdy->ti->time_integration_method) { case TDySNES: - ierr = SNESCreate(PETSC_COMM_WORLD,&snes); + ierr = SNESCreate(comm,&snes); CHKERRQ(ierr); ierr = TDySetSNESFunction(snes,tdy); CHKERRQ(ierr); ierr = TDySetSNESJacobian(snes,tdy); CHKERRQ(ierr); @@ -112,7 +119,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { tdy->ti->snes = snes; break; case TDyTS: - ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr); + ierr = TSCreate(comm,&ts); CHKERRQ(ierr); // ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr); // ierr = TSSetType(ts,TSPSEUDO); CHKERRQ(ierr); // ierr = TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,NULL); CHKERRQ(ierr); @@ -163,7 +170,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { ierr = TDySetup(tdy); CHKERRQ(ierr); break; } - PetscPrintf(PETSC_COMM_WORLD,"tdy->ti->time_integration_method = %d\n", + PetscPrintf(comm,"tdy->ti->time_integration_method = %d\n", tdy->ti->time_integration_method); size_t len; ierr = PetscStrlen(tdy->io->ic_filename, &len); CHKERRQ(ierr); From 7413c429c960840420951ad6464bff82c0c35ce8 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 11:33:57 -0700 Subject: [PATCH 22/23] Discovered some differences in the way TDyDriver sets things up. --- include/private/tdycoreimpl.h | 7 ++++++- include/private/tdyrichardsimpl.h | 2 +- include/private/tdythimpl.h | 3 ++- src/tdydriver.c | 4 +++- src/tdyth.c | 7 +++++-- 5 files changed, 17 insertions(+), 6 deletions(-) diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index d2df81e2..1c32167b 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -25,6 +25,11 @@ struct _TDyOps { PetscErrorCode (*destroy)(void*); // Creates a DM for a particular simulation (optional). + // Arguments: + // 1. A valid MPI communicator + // 2. A context pointer storing data specifically for the constructor + // function + // 3. A location for the newly created DM. // We pass the dycore as the first argument here because we don't expect // the caller to know implementation details. PetscErrorCode (*create_dm)(MPI_Comm, void*, DM*); @@ -35,7 +40,7 @@ struct _TDyOps { // Called by TDySetFromOptions -- sets implementation-specific options // from command-line arguments. - // FIXME: convert the first arg here to void* when we've moved specific data + // FIXME: convert the arg here to void* when we've moved specific data // FIXME: out of TDy. PetscErrorCode (*set_from_options)(TDy); diff --git a/include/private/tdyrichardsimpl.h b/include/private/tdyrichardsimpl.h index 67608fe7..753716f7 100644 --- a/include/private/tdyrichardsimpl.h +++ b/include/private/tdyrichardsimpl.h @@ -4,7 +4,7 @@ #include #include -PETSC_INTERN PetscErrorCode TDyRichardsInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyRichardsInitialize(TDy); // TDyDriver setup fn PETSC_INTERN PetscErrorCode TDyRichardsTSPostStep(TS); PETSC_INTERN PetscErrorCode TDyRichardsSNESPostCheck(SNESLineSearch,Vec,Vec,Vec, PetscBool*,PetscBool*, diff --git a/include/private/tdythimpl.h b/include/private/tdythimpl.h index fe269cf5..fb8ac240 100644 --- a/include/private/tdythimpl.h +++ b/include/private/tdythimpl.h @@ -4,7 +4,8 @@ #include #include -PETSC_INTERN PetscErrorCode TDyTH_MPFA_O_Setup(TDy, DM); +PETSC_INTERN PetscErrorCode TDyTH_MPFA_O_Setup(TDy tdy, DM dm); +PETSC_INTERN PetscErrorCode TDyTHInitialize(TDy); // TDyDriver setup fn PETSC_INTERN PetscErrorCode TDyTHTSPostStep(TS); PETSC_INTERN PetscErrorCode TDyTHSNESPostCheck(SNESLineSearch,Vec,Vec,Vec, PetscBool*,PetscBool*, diff --git a/src/tdydriver.c b/src/tdydriver.c index 036205cf..b5d3ecd5 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -154,6 +154,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { ierr = TSSetPostStep(ts,TDyRichardsTSPostStep); CHKERRQ(ierr); break; } + // FIXME: This is a different path from the one used by TDycore proper. ierr = TDyRichardsInitialize(tdy); CHKERRQ(ierr); break; case TH: @@ -167,7 +168,8 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case TDyTS: ierr = TSSetPostStep(ts,TDyTHTSPostStep); CHKERRQ(ierr); } - ierr = TDySetup(tdy); CHKERRQ(ierr); + // FIXME: This is a different path from the one used by TDycore proper. + ierr = TDyTHInitialize(tdy); CHKERRQ(ierr); break; } PetscPrintf(comm,"tdy->ti->time_integration_method = %d\n", diff --git a/src/tdyth.c b/src/tdyth.c index 09ed65ec..512d6357 100644 --- a/src/tdyth.c +++ b/src/tdyth.c @@ -9,11 +9,14 @@ PetscErrorCode TDyTHInitialize(TDy tdy) { PetscInt local_size; PetscFunctionBegin; - PetscPrintf(PETSC_COMM_WORLD,"Running TH mode.\n"); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + PetscPrintf(comm,"Running TH mode.\n"); if (tdy->options.init_with_random_field) { PetscRandom rand; - ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand); CHKERRQ(ierr); + ierr = PetscRandomCreate(comm,&rand); CHKERRQ(ierr); ierr = VecGetLocalSize(tdy->solution,&local_size); CHKERRQ(ierr); ierr = DMCreateGlobalVector(tdy->dm,&temp_vec); CHKERRQ(ierr); // pressure From 8ccbf056162faced249c0f0bbf902598692919df Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 18 Oct 2021 14:44:31 -0700 Subject: [PATCH 23/23] Fixed a spelling error in the timer report. --- src/tdytimers.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tdytimers.c b/src/tdytimers.c index 2b0e5eb8..5bcfd8c3 100644 --- a/src/tdytimers.c +++ b/src/tdytimers.c @@ -140,7 +140,7 @@ PetscErrorCode TDyWriteTimingProfile(const char* filename) { const char* disc_name = TDyDiscretizations[md->discretization]; FILE* f = fopen("tdycore_profile.csv", "a"); fprintf(f, "METADATA\n"); - fprintf(f, "Mode,Discretiztion,NumProc,NumCells\n"); + fprintf(f, "Mode,Discretization,NumProc,NumCells\n"); fprintf(f, "%s,%s,%d,%d", mode_name, disc_name, md->num_proc, num_cells); fclose(f); } else { // rank > 0