Merge branch 'rosie/add-input' (PR #189)

Adds support for reading initial condition and permeability dataset from HDF5.
TDycores-Project · Aug 13, 2021 · 4586201 · 4586201
2 parents a585b2c + a9078df
commit 4586201
Show file tree

Hide file tree

Showing 10 changed files with 323 additions and 10 deletions.
diff --git a/demo/richards/initial.h5 b/demo/richards/initial.h5
diff --git a/demo/richards/richards-driver-snes-initial-cond-anisotropic-k.regression.gold b/demo/richards/richards-driver-snes-initial-cond-anisotropic-k.regression.gold
@@ -0,0 +1,5 @@
+-- PRESSURE: Liquid Pressure --
+      Max:   9.9000539026076e+04
+      Min:   9.9000488451729e+04
+     Mean:   9.9000512636697e+04
+        0:   9.9000488451729e+04
diff --git a/demo/richards/richards-driver-snes-initial-cond-isotropic-k.regression.gold b/demo/richards/richards-driver-snes-initial-cond-isotropic-k.regression.gold
@@ -0,0 +1,5 @@
+-- PRESSURE: Liquid Pressure --
+      Max:   9.9000581976854e+04
+      Min:   9.9000528668226e+04
+     Mean:   9.9000555027803e+04
+        0:   9.9000528668226e+04
diff --git a/demo/richards/richards.cfg b/demo/richards/richards.cfg
@@ -1,6 +1,8 @@
 [suites]
 standard=
 	richards-driver-snes-prob1
+	richards-driver-snes-initial-cond-isotropic-k
+	richards-driver-snes-initial-cond-anisotropic-k
 	#richards-driver-ts-prob1
 
 standard_parallel=
@@ -25,3 +27,8 @@ 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
 
+[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
+
+[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
diff --git a/include/private/tdydiscretization.h b/include/private/tdydiscretization.h
@@ -10,6 +10,7 @@ PETSC_INTERN PetscErrorCode TDyCreateJacobianMatrix(TDy,Mat*);
 PETSC_INTERN PetscErrorCode TDyGlobalToNatural(TDy,Vec,Vec);
 PETSC_INTERN PetscErrorCode TDyGlobalToLocal(TDy,Vec,Vec);
 PETSC_INTERN PetscErrorCode TDyNaturalToGlobal(TDy,Vec,Vec);
+PETSC_INTERN PetscErrorCode TDyNaturaltoLocal(TDy,Vec,Vec*);
 
 #endif
 
diff --git a/include/private/tdyioimpl.h b/include/private/tdyioimpl.h
@@ -6,23 +6,34 @@
 struct _p_TDyIO {
   PetscBool io_process;
   PetscBool print_intermediate;
+  PetscBool anisotropic_permeability;
   char filename[PETSC_MAX_PATH_LEN];
   char zonalVarNames[2][PETSC_MAX_PATH_LEN];
   int num_vars;
   TDyIOFormat format;
   int num_times;
+  char permeability_filename[PETSC_MAX_PATH_LEN];
+  char permeability_dataset[PETSC_MAX_PATH_LEN];
+  char porosity_filename[PETSC_MAX_PATH_LEN];
+  char porosity_dataset[PETSC_MAX_PATH_LEN];
+  char ic_filename[PETSC_MAX_PATH_LEN];
+  char ic_dataset[PETSC_MAX_PATH_LEN];
 };
 
 PETSC_INTERN PetscErrorCode TDyIOCreate(TDyIO*);
+PETSC_INTERN PetscErrorCode TDyIOReadPermeability(TDy);
+PETSC_INTERN PetscErrorCode TDyIOReadPorosity(TDy);
+PETSC_INTERN PetscErrorCode TDyIOReadIC(TDy);
+PETSC_INTERN PetscErrorCode TDyIOReadVariable(TDy,char*,char*,PetscReal**);
 PETSC_INTERN PetscErrorCode TDyIOInitializeExodus(char*,char*[],DM ,int);
 PETSC_INTERN PetscErrorCode TDyIOAddExodusTime(char*,PetscReal,DM,TDyIO);
 PETSC_INTERN PetscErrorCode TDyIOWriteExodusVar(char*,Vec,char*,TDyIO,PetscReal);
-PETSC_EXTERN PetscErrorCode TDyIOInitializeHDF5(char*,DM);
-PETSC_EXTERN PetscErrorCode TDyIOWriteHDF5Var(char*,DM,Vec,char*,PetscReal);
-PETSC_EXTERN PetscErrorCode TDyIOWriteAsciiViewer(Vec,PetscReal,char*);
-PETSC_EXTERN PetscErrorCode TDyIOWriteXMFHeader(PetscInt,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode TDyIOWriteXMFAttribute(char*,char*,PetscInt);
-PETSC_EXTERN PetscErrorCode TDyIOWriteXMFFooter();
+PETSC_INTERN PetscErrorCode TDyIOInitializeHDF5(char*,DM);
+PETSC_INTERN PetscErrorCode TDyIOWriteHDF5Var(char*,DM,Vec,char*,PetscReal);
+PETSC_INTERN PetscErrorCode TDyIOWriteAsciiViewer(Vec,PetscReal,char*);
+PETSC_INTERN PetscErrorCode TDyIOWriteXMFHeader(PetscInt,PetscInt,PetscInt,PetscInt);
+PETSC_INTERN PetscErrorCode TDyIOWriteXMFAttribute(char*,char*,PetscInt);
+PETSC_INTERN PetscErrorCode TDyIOWriteXMFFooter();
 PETSC_INTERN PetscErrorCode TDyIODestroy(TDyIO*);
 
 #endif

diff --git a/src/tdydiscretization.c b/src/tdydiscretization.c
@@ -146,3 +146,28 @@ PetscErrorCode TDyNaturalToGlobal(TDy tdy, Vec natural, Vec global){
 
   PetscFunctionReturn(0);
 }
+
+/* -------------------------------------------------------------------------- */
+/// Performs scatter of a natural vector to a local vector
+///
+/// @param [in] tdy     A TDy struct
+/// @param [in] natural A PETSc vector
+/// @param [out] local   A PETSc vector
+/// @returns 0 on success, or a non-zero error code on failure
+PetscErrorCode TDyNaturaltoLocal(TDy tdy,Vec natural, Vec *local) {
+
+  PetscFunctionBegin;
+
+  PetscErrorCode ierr;
+  Vec global;
+
+  ierr = TDyCreateGlobalVector(tdy, &global);CHKERRQ(ierr);
+
+  ierr = TDyNaturalToGlobal(tdy, natural, global);CHKERRQ(ierr);
+  ierr = TDyGlobalToLocal(tdy, global, *local); CHKERRQ(ierr);
+
+  ierr = VecDestroy(&global); CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+
+}
diff --git a/src/tdydriver.c b/src/tdydriver.c
@@ -3,6 +3,7 @@
 #include <private/tdyrichardsimpl.h>
 #include <private/tdymaterialpropertiesimpl.h>
 #include <private/tdythimpl.h>
+#include <private/tdyioimpl.h>
 #include <tdytimers.h>
 #include <private/tdycharacteristiccurvesimpl.h>
 
@@ -34,11 +35,23 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) {
   }
 
   if (!TDyIsPorositySet(tdy)) {
-    ierr = TDySetPorosityFunction(tdy,TDyConstantPorosityFunction,PETSC_NULL); CHKERRQ(ierr);
+    size_t len;
+    ierr = PetscStrlen(tdy->io->porosity_filename, &len); CHKERRQ(ierr);
+    if (!len){
+      ierr = TDySetPorosityFunction(tdy,TDyConstantPorosityFunction,PETSC_NULL);CHKERRQ(ierr);
+    } else {
+      ierr = TDyIOReadPorosity(tdy);CHKERRQ(ierr);
+    }
   }
 
   if (!TDyIsPermeabilitySet(tdy)){
-    ierr = TDySetPermeabilityFunction(tdy,TDyConstantPermeabilityFunction,PETSC_NULL); CHKERRQ(ierr);
+    size_t len;
+    ierr = PetscStrlen(tdy->io->permeability_filename, &len); CHKERRQ(ierr);
+    if (!len){
+      ierr = TDySetPermeabilityFunction(tdy,TDyConstantPermeabilityFunction,PETSC_NULL); CHKERRQ(ierr);     
+    } else {
+      ierr = TDyIOReadPermeability(tdy);CHKERRQ(ierr);
+    }
   }
 
   if (tdy->options.mode == TH) {
@@ -153,6 +166,11 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) {
   }
   PetscPrintf(PETSC_COMM_WORLD,"tdy->ti->time_integration_method = %d\n",
               tdy->ti->time_integration_method);
+   size_t len;
+   ierr = PetscStrlen(tdy->io->ic_filename, &len); CHKERRQ(ierr);
+   if (len){
+     TDyIOReadIC(tdy); CHKERRQ(ierr);
+   }
   // finish set of time integrators
   switch(tdy->ti->time_integration_method) {
     case TDySNES: