From 9cd1391958c36186193e1d90856ce38131f043f0 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 9 Aug 2021 12:32:47 -0700 Subject: [PATCH 01/10] Pull all options into TDySetFromOptions. --- src/tdycore.c | 129 +++++++++++++++++++++++++++++--------------------- 1 file changed, 75 insertions(+), 54 deletions(-) diff --git a/src/tdycore.c b/src/tdycore.c index 8f7ce54f..ff2bb2c6 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -181,7 +181,7 @@ PetscErrorCode TDyFinalize() { PetscFunctionReturn(0); } -PetscErrorCode TDyProcessOptions(TDyOptions *options){ +PetscErrorCode TDySetOptionDefaults(TDyOptions *options){ PetscFunctionBegin; PetscErrorCode ierr; @@ -198,25 +198,6 @@ PetscErrorCode TDyProcessOptions(TDyOptions *options){ options->default_vangenuchten_m=0.8; options->default_vangenutchen_alpha=1.e-4; - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: General options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_gravity", "Value of gravity", NULL, options->gravity_constant, &options->gravity_constant, NULL); CHKERRQ(ierr); - ierr = PetscOptionsEnd(); CHKERRQ(ierr); - - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Material propertiy options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_porosity", "Value of porosity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_permability", "Value of permeability", NULL, options->default_permeability, &options->default_permeability, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_soil_density", "Value of soil density", NULL, options->default_soil_density, &options->default_soil_density, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_soil_specific_heat", "Value of soil specific heat", NULL, options->default_soil_specific_heat, &options->default_soil_specific_heat, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_thermal_conductivity", "Value of thermal conductivity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); - ierr = PetscOptionsEnd(); CHKERRQ(ierr); - - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Characteristic curve options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_residual_satuaration", "Value of residual saturation", NULL, options->default_residual_saturation, &options->default_residual_saturation, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_gardner_param_n", "Value of Gardner n parameter", NULL, options->default_gardner_n, &options->default_gardner_n, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_vangenuchten_param_m", "Value of VanGenuchten m parameter", NULL, options->default_vangenuchten_m, &options->default_vangenuchten_m, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_vangenuchten_param_alpha", "Value of VanGenuchten alpha parameter", NULL, options->default_vangenutchen_alpha, &options->default_vangenutchen_alpha, NULL); CHKERRQ(ierr); - ierr = PetscOptionsEnd(); CHKERRQ(ierr); - PetscFunctionReturn(0); } @@ -236,7 +217,7 @@ PetscErrorCode TDyCreate(TDy *_tdy) { *_tdy = tdy; tdy->setupflags |= TDyCreated; - ierr = TDyProcessOptions(&tdy->options); + ierr = TDySetOptionDefaults(&tdy->options); ierr = TDyIOCreate(&tdy->io); CHKERRQ(ierr); @@ -534,8 +515,11 @@ PetscErrorCode TDyView(TDy tdy,PetscViewer viewer) { PetscFunctionReturn(0); } +/// Sets options for the dycore based on command line arguments supplied by a +/// user. TDySetFromOptions must be called before TDySetupNumericalMethods, +/// since the latter uses options specified by the former. +/// @param tdy The dycore instance PetscErrorCode TDySetFromOptions(TDy tdy) { - // must preceed TDySetupNumericalMethods() as it sets options used in TDySetupNumericalMethods() PetscErrorCode ierr; PetscBool flg; TDyMethod method = WY; @@ -545,20 +529,82 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { TDyMPFAOGmatrixMethod gmatrixmethod = MPFAO_GMATRIX_DEFAULT; TDyMPFAOBoundaryConditionType bctype = MPFAO_DIRICHLET_BC; PetscFunctionBegin; + if ((tdy->setupflags & TDySetupFinished) != 0) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetupNumericalMethods()"); } + + TDyOptions *options = &tdy->options; PetscValidHeaderSpecific(tdy,TDY_CLASSID,1); ierr = PetscObjectOptionsBegin((PetscObject)tdy); CHKERRQ(ierr); - ierr = PetscOptionsEnum("-tdy_method","Discretization method", - "TDySetDiscretizationMethod",TDyMethods,(PetscEnum)method,(PetscEnum *)&method, + + // Material property options + ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Material property options",""); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_porosity", "Value of porosity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_permability", "Value of permeability", NULL, options->default_permeability, &options->default_permeability, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_soil_density", "Value of soil density", NULL, options->default_soil_density, &options->default_soil_density, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_soil_specific_heat", "Value of soil specific heat", NULL, options->default_soil_specific_heat, &options->default_soil_specific_heat, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_thermal_conductivity", "Value of thermal conductivity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); + ierr = PetscOptionsEnd(); CHKERRQ(ierr); + + // Characteristic curve options + ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Characteristic curve options",""); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_residual_satuaration", "Value of residual saturation", NULL, options->default_residual_saturation, &options->default_residual_saturation, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_gardner_param_n", "Value of Gardner n parameter", NULL, options->default_gardner_n, &options->default_gardner_n, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_vangenuchten_param_m", "Value of VanGenuchten m parameter", NULL, options->default_vangenuchten_m, &options->default_vangenuchten_m, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_vangenuchten_param_alpha", "Value of VanGenuchten alpha parameter", NULL, options->default_vangenutchen_alpha, &options->default_vangenutchen_alpha, NULL); CHKERRQ(ierr); + ierr = PetscOptionsEnd(); CHKERRQ(ierr); + + // Model options + ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Model options",""); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_gravity", "Magntitude of gravity vector", NULL, + options->gravity_constant, &options->gravity_constant, + NULL); CHKERRQ(ierr); + ierr = PetscOptionsEnum("-tdy_water_density","Water density vertical profile", + "TDySetWaterDensityType",TDyWaterDensityTypes,(PetscEnum)densitytype,(PetscEnum *)&densitytype, &flg); CHKERRQ(ierr); + if (flg) { + ierr = TDySetWaterDensityType(tdy,densitytype); CHKERRQ(ierr); + } + + ierr = PetscOptionsEnum("-tdy_mode","Flow mode", + "TDySetMode",TDyModes,(PetscEnum)mode,(PetscEnum *)&mode, + &flg); CHKERRQ(ierr); + + if (flg && (mode != tdy->mode)) { + ierr = TDySetMode(tdy,mode); CHKERRQ(ierr); + } + ierr = PetscOptionsEnd(); CHKERRQ(ierr); + + // Numerics options + ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); + ierr = PetscOptionsEnum("-tdy_method","Discretization method", + "TDySetDiscretizationMethod",TDyMethods, + (PetscEnum)method,(PetscEnum *)&method, &flg); CHKERRQ(ierr); if (flg && (method != tdy->method)) { ierr = TDySetDiscretizationMethod(tdy,method); CHKERRQ(ierr); } - ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type", + ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods", "TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&qtype, &flg); CHKERRQ(ierr); + ierr = PetscOptionsBool("-tdy_tpf_allow_all_meshes", + "Enable to allow non-orthgonal meshes in finite volume TPF method", + "",tdy->allow_unsuitable_mesh, + &(tdy->allow_unsuitable_mesh),NULL); CHKERRQ(ierr); + ierr = PetscOptionsEnum("-tdy_mpfao_gmatrix_method","MPFA-O gmatrix method", + "TDySetMPFAOGmatrixMethod",TDyMPFAOGmatrixMethods,(PetscEnum)gmatrixmethod,(PetscEnum *)&gmatrixmethod, + &flg); CHKERRQ(ierr); + if (flg && (gmatrixmethod != tdy->mpfao_gmatrix_method)) { + ierr = TDySetMPFAOGmatrixMethod(tdy,gmatrixmethod); CHKERRQ(ierr); + } + ierr = PetscOptionsEnum("-tdy_mpfao_boundary_condition_type","MPFA-O boundary condition type", + "TDySetMPFAOBoundaryConditionType",TDyMPFAOBoundaryConditionTypes,(PetscEnum)bctype,(PetscEnum *)&bctype, + &flg); CHKERRQ(ierr); + if (flg && (bctype != tdy->mpfao_bc_type)) { + ierr = TDySetMPFAOBoundaryConditionType(tdy,bctype); CHKERRQ(ierr); + } + ierr = PetscOptionsEnd(); CHKERRQ(ierr); + if (flg && (qtype != tdy->qtype)) { ierr = TDySetQuadratureType(tdy,qtype); CHKERRQ(ierr); } ierr = PetscOptionsBool("-tdy_init_with_random_field", "Initialize solution with a random field","", @@ -567,10 +613,6 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsGetString(NULL,NULL,"-tdy_init_file", tdy->init_file, sizeof(tdy->init_file), &tdy->init_from_file); CHKERRQ(ierr); - ierr = PetscOptionsBool("-tdy_tpf_allow_unsuitable_mesh", - "Enable to allow non-orthgonal meshes in tpf","",tdy->allow_unsuitable_mesh, - &(tdy->allow_unsuitable_mesh),NULL); CHKERRQ(ierr); - ierr = PetscOptionsBool("-tdy_regression_test", "Enable output of a regression file","",tdy->regression_testing, &(tdy->regression_testing),NULL); CHKERRQ(ierr); @@ -579,27 +621,7 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { "Enable output of mesh attributes","",tdy->output_mesh, &(tdy->output_mesh),NULL); CHKERRQ(ierr); - ierr = PetscOptionsEnum("-tdy_water_density","Water density type", - "TDySetWaterDensityType",TDyWaterDensityTypes,(PetscEnum)densitytype,(PetscEnum *)&densitytype, - &flg); CHKERRQ(ierr); - if (flg) {ierr = TDySetWaterDensityType(tdy,densitytype); CHKERRQ(ierr);} - - ierr = PetscOptionsEnum("-tdy_mode","Flow mode", - "TDySetMode",TDyModes,(PetscEnum)mode,(PetscEnum *)&mode, - &flg); CHKERRQ(ierr); - - if (flg && (mode != tdy->mode)) { ierr = TDySetMode(tdy,mode); CHKERRQ(ierr); } - - ierr = PetscOptionsEnum("-tdy_mpfao_gmatrix_method","MPFA-O gmatrix method", - "TDySetMPFAOGmatrixMethod",TDyMPFAOGmatrixMethods,(PetscEnum)gmatrixmethod,(PetscEnum *)&gmatrixmethod, - &flg); CHKERRQ(ierr); - if (flg && (gmatrixmethod != tdy->mpfao_gmatrix_method)) { ierr = TDySetMPFAOGmatrixMethod(tdy,gmatrixmethod); CHKERRQ(ierr); } - - ierr = PetscOptionsEnum("-tdy_mpfao_boundary_condition_type","MPFA-O boundary condition type", - "TDySetMPFAOBoundaryConditionType",TDyMPFAOBoundaryConditionTypes,(PetscEnum)bctype,(PetscEnum *)&bctype, - &flg); CHKERRQ(ierr); - if (flg && (bctype != tdy->mpfao_bc_type)) { ierr = TDySetMPFAOBoundaryConditionType(tdy,bctype); CHKERRQ(ierr); } - + // Wrap up and indicate that options are set. ierr = PetscOptionsEnd(); CHKERRQ(ierr); tdy->setupflags |= TDyOptionsSet; @@ -646,7 +668,7 @@ PetscErrorCode TDySetupDiscretizationScheme(TDy tdy) { } PetscErrorCode TDySetupNumericalMethods(TDy tdy) { - /* must follow TDySetFromOptions() is it relies upon options set by + /* must follow TDySetFromOptions() is it relies upon options set by TDySetFromOptions */ PetscErrorCode ierr; PetscFunctionBegin; @@ -655,9 +677,9 @@ PetscErrorCode TDySetupNumericalMethods(TDy tdy) { } TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); - ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); + ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); if (tdy->regression_testing) { - /* must come after Sections are set up in + /* must come after Sections are set up in TDySetupDiscretizationScheme->XXXInitialize */ ierr = TDyRegressionInitialize(tdy); CHKERRQ(ierr); } @@ -1528,4 +1550,3 @@ const char* NewCString(char* f_str_ptr, int f_str_len) { return str; } } - From 58c3d968776fc0de2d283841961837f0bffab496 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 9 Aug 2021 13:37:24 -0700 Subject: [PATCH 02/10] Removed default_ prefix from options. --- include/private/tdyoptions.h | 24 ++++---- src/tdycore.c | 104 +++++++++++++++++++---------------- src/tdymaterialproperties.c | 14 ++--- 3 files changed, 75 insertions(+), 67 deletions(-) diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index 30362a00..931306e2 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -6,19 +6,19 @@ typedef struct{ PetscReal gravity_constant; - // Default values for material properties - PetscReal default_porosity; - PetscReal default_permeability; - PetscReal default_soil_density; - PetscReal default_soil_specific_heat; - PetscReal default_thermal_conductivity; + // Material properties + PetscReal porosity; + PetscReal permeability; + PetscReal soil_density; + PetscReal soil_specific_heat; + PetscReal thermal_conductivity; - // Default values for characteristic curves - PetscReal default_residual_saturation; - PetscReal default_gardner_n; - PetscReal default_vangenuchten_m; - PetscReal default_vangenutchen_alpha; + // Characteristic curves + PetscReal residual_saturation; + PetscReal gardner_n; + PetscReal vangenuchten_m; + PetscReal vangenuchten_alpha; } TDyOptions; -#endif \ No newline at end of file +#endif diff --git a/src/tdycore.c b/src/tdycore.c index ff2bb2c6..00943cf2 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -187,16 +187,16 @@ PetscErrorCode TDySetOptionDefaults(TDyOptions *options){ options->gravity_constant = 9.8068; - options->default_porosity=0.25; - options->default_permeability=1.e-12; - options->default_soil_density=2650.; - options->default_soil_specific_heat=1000.0; - options->default_thermal_conductivity=1.0; + options->porosity=0.25; + options->permeability=1.e-12; + options->soil_density=2650.; + options->soil_specific_heat=1000.0; + options->thermal_conductivity=1.0; - options->default_residual_saturation=0.15; - options->default_gardner_n=0.5; - options->default_vangenuchten_m=0.8; - options->default_vangenutchen_alpha=1.e-4; + options->residual_saturation=0.15; + options->gardner_n=0.5; + options->vangenuchten_m=0.8; + options->vangenuchten_alpha=1.e-4; PetscFunctionReturn(0); } @@ -294,14 +294,14 @@ PetscErrorCode TDyMalloc(TDy tdy) { PetscReal mualem_poly_low = 0.99; for (c=0; csr[c] = options->default_residual_saturation; - cc->gardner_n[c] = options->default_gardner_n; - cc->gardner_m[c] = options->default_vangenuchten_m; - cc->vg_m[c] = options->default_vangenuchten_m; + cc->sr[c] = options->residual_saturation; + cc->gardner_n[c] = options->gardner_n; + cc->gardner_m[c] = options->vangenuchten_m; + cc->vg_m[c] = options->vangenuchten_m; cc->mualem_poly_low[c] = mualem_poly_low; - cc->mualem_m[c] = options->default_vangenuchten_m; - cc->irmay_m[c] = options->default_vangenuchten_m; - cc->vg_alpha[c] = options->default_vangenutchen_alpha; + cc->mualem_m[c] = options->vangenuchten_m; + cc->irmay_m[c] = options->vangenuchten_m; + cc->vg_alpha[c] = options->vangenuchten_alpha; cc->SatFuncType[c] = SAT_FUNC_VAN_GENUCHTEN; cc->RelPermFuncType[c] = REL_PERM_FUNC_MUALEM; cc->Kr[c] = 0.0; @@ -521,13 +521,6 @@ PetscErrorCode TDyView(TDy tdy,PetscViewer viewer) { /// @param tdy The dycore instance PetscErrorCode TDySetFromOptions(TDy tdy) { PetscErrorCode ierr; - PetscBool flg; - TDyMethod method = WY; - TDyMode mode = RICHARDS; - TDyQuadratureType qtype = FULL; - TDyWaterDensityType densitytype = WATER_DENSITY_CONSTANT; - TDyMPFAOGmatrixMethod gmatrixmethod = MPFAO_GMATRIX_DEFAULT; - TDyMPFAOBoundaryConditionType bctype = MPFAO_DIRICHLET_BC; PetscFunctionBegin; if ((tdy->setupflags & TDySetupFinished) != 0) { @@ -537,81 +530,96 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { TDyOptions *options = &tdy->options; PetscValidHeaderSpecific(tdy,TDY_CLASSID,1); ierr = PetscObjectOptionsBegin((PetscObject)tdy); CHKERRQ(ierr); + PetscBool flag; // Material property options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Material property options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_porosity", "Value of porosity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_permability", "Value of permeability", NULL, options->default_permeability, &options->default_permeability, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_soil_density", "Value of soil density", NULL, options->default_soil_density, &options->default_soil_density, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_soil_specific_heat", "Value of soil specific heat", NULL, options->default_soil_specific_heat, &options->default_soil_specific_heat, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_thermal_conductivity", "Value of thermal conductivity", NULL, options->default_porosity, &options->default_porosity, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_porosity", "Value of porosity", NULL, options->porosity, &options->porosity, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_permability", "Value of permeability", NULL, options->permeability, &options->permeability, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_soil_density", "Value of soil density", NULL, options->soil_density, &options->soil_density, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_soil_specific_heat", "Value of soil specific heat", NULL, options->soil_specific_heat, &options->soil_specific_heat, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_thermal_conductivity", "Value of thermal conductivity", NULL, options->porosity, &options->porosity, NULL); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Characteristic curve options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Characteristic curve options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_residual_satuaration", "Value of residual saturation", NULL, options->default_residual_saturation, &options->default_residual_saturation, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_gardner_param_n", "Value of Gardner n parameter", NULL, options->default_gardner_n, &options->default_gardner_n, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_vangenuchten_param_m", "Value of VanGenuchten m parameter", NULL, options->default_vangenuchten_m, &options->default_vangenuchten_m, NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_vangenuchten_param_alpha", "Value of VanGenuchten alpha parameter", NULL, options->default_vangenutchen_alpha, &options->default_vangenutchen_alpha, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_residual_satuaration", "Value of residual saturation", NULL, options->residual_saturation, &options->residual_saturation, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_gardner_param_n", "Value of Gardner n parameter", NULL, options->gardner_n, &options->gardner_n, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_vangenuchten_param_m", "Value of VanGenuchten m parameter", NULL, options->vangenuchten_m, &options->vangenuchten_m, NULL); CHKERRQ(ierr); + ierr = PetscOptionsReal("-tdy_vangenuchten_param_alpha", "Value of VanGenuchten alpha parameter", NULL, options->vangenuchten_alpha, &options->vangenuchten_alpha, NULL); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Model options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Model options",""); CHKERRQ(ierr); - ierr = PetscOptionsReal("-tdy_gravity", "Magntitude of gravity vector", NULL, + ierr = PetscOptionsReal("-tdy_gravity", "Magnitude of gravity vector", NULL, options->gravity_constant, &options->gravity_constant, NULL); CHKERRQ(ierr); + TDyWaterDensityType densitytype = WATER_DENSITY_CONSTANT; ierr = PetscOptionsEnum("-tdy_water_density","Water density vertical profile", - "TDySetWaterDensityType",TDyWaterDensityTypes,(PetscEnum)densitytype,(PetscEnum *)&densitytype, - &flg); CHKERRQ(ierr); - if (flg) { + "TDySetWaterDensityType",TDyWaterDensityTypes, + (PetscEnum)densitytype,(PetscEnum *)&densitytype, + &flag); CHKERRQ(ierr); + if (flag) { ierr = TDySetWaterDensityType(tdy,densitytype); CHKERRQ(ierr); } - + TDyMode mode = RICHARDS; ierr = PetscOptionsEnum("-tdy_mode","Flow mode", "TDySetMode",TDyModes,(PetscEnum)mode,(PetscEnum *)&mode, - &flg); CHKERRQ(ierr); + &flag); CHKERRQ(ierr); - if (flg && (mode != tdy->mode)) { + if (flag && (mode != tdy->mode)) { ierr = TDySetMode(tdy,mode); CHKERRQ(ierr); } + ierr = PetscOptionsReal("-tdy_pressure_bc", + "Assigns a named pressure function to the domain boundary", + NULL, options->gravity_constant, &options->gravity_constant, + NULL); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); + TDyMethod method = WY; ierr = PetscOptionsEnum("-tdy_method","Discretization method", "TDySetDiscretizationMethod",TDyMethods, - (PetscEnum)method,(PetscEnum *)&method, &flg); CHKERRQ(ierr); - if (flg && (method != tdy->method)) { + (PetscEnum)method,(PetscEnum *)&method, &flag); CHKERRQ(ierr); + if (flag && (method != tdy->method)) { ierr = TDySetDiscretizationMethod(tdy,method); CHKERRQ(ierr); } + TDyQuadratureType qtype = FULL; ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods", "TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&qtype, - &flg); CHKERRQ(ierr); + &flag); CHKERRQ(ierr); + if (flag && (qtype != tdy->qtype)) { + ierr = TDySetQuadratureType(tdy,qtype); CHKERRQ(ierr); + } ierr = PetscOptionsBool("-tdy_tpf_allow_all_meshes", "Enable to allow non-orthgonal meshes in finite volume TPF method", "",tdy->allow_unsuitable_mesh, &(tdy->allow_unsuitable_mesh),NULL); CHKERRQ(ierr); + TDyMPFAOGmatrixMethod gmatrixmethod = MPFAO_GMATRIX_DEFAULT; ierr = PetscOptionsEnum("-tdy_mpfao_gmatrix_method","MPFA-O gmatrix method", "TDySetMPFAOGmatrixMethod",TDyMPFAOGmatrixMethods,(PetscEnum)gmatrixmethod,(PetscEnum *)&gmatrixmethod, - &flg); CHKERRQ(ierr); - if (flg && (gmatrixmethod != tdy->mpfao_gmatrix_method)) { + &flag); CHKERRQ(ierr); + if (flag && (gmatrixmethod != tdy->mpfao_gmatrix_method)) { ierr = TDySetMPFAOGmatrixMethod(tdy,gmatrixmethod); CHKERRQ(ierr); } + TDyMPFAOBoundaryConditionType bctype = MPFAO_DIRICHLET_BC; ierr = PetscOptionsEnum("-tdy_mpfao_boundary_condition_type","MPFA-O boundary condition type", "TDySetMPFAOBoundaryConditionType",TDyMPFAOBoundaryConditionTypes,(PetscEnum)bctype,(PetscEnum *)&bctype, - &flg); CHKERRQ(ierr); - if (flg && (bctype != tdy->mpfao_bc_type)) { + &flag); CHKERRQ(ierr); + if (flag && (bctype != tdy->mpfao_bc_type)) { ierr = TDySetMPFAOBoundaryConditionType(tdy,bctype); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); - if (flg && (qtype != tdy->qtype)) { ierr = TDySetQuadratureType(tdy,qtype); CHKERRQ(ierr); } ierr = PetscOptionsBool("-tdy_init_with_random_field", "Initialize solution with a random field","", tdy->init_with_random_field, &(tdy->init_with_random_field),NULL); CHKERRQ(ierr); - ierr = PetscOptionsGetString(NULL,NULL,"-tdy_init_file", tdy->init_file, sizeof(tdy->init_file), &tdy->init_from_file); CHKERRQ(ierr); + ierr = PetscOptionsGetString(NULL,NULL,"-tdy_init_file", tdy->init_file, + sizeof(tdy->init_file), + &tdy->init_from_file); CHKERRQ(ierr); ierr = PetscOptionsBool("-tdy_regression_test", "Enable output of a regression file","",tdy->regression_testing, diff --git a/src/tdymaterialproperties.c b/src/tdymaterialproperties.c index c9c7daed..bcaa554e 100644 --- a/src/tdymaterialproperties.c +++ b/src/tdymaterialproperties.c @@ -18,7 +18,7 @@ PetscErrorCode MaterialPropertiesCreate(PetscInt ndim, PetscInt ncells, Material ierr = PetscMalloc(ncells*sizeof(PetscReal),&(*_matprop)->porosity); CHKERRQ(ierr); ierr = PetscMalloc(ncells*sizeof(PetscReal),&(*_matprop)->Cr ); CHKERRQ(ierr); ierr = PetscMalloc(ncells*sizeof(PetscReal),&(*_matprop)->rhosoil ); CHKERRQ(ierr); - + (*_matprop)-> permeability_is_set = 0; (*_matprop)-> porosity_is_set = 0; (*_matprop)-> thermal_conductivity_is_set = 0; @@ -310,7 +310,7 @@ PetscErrorCode TDySetCellPermeability(TDy tdy,PetscInt c,PetscReal *K) { ierr = DMPlexGetHeightStratum(tdy->dm,0,&cStart,&i); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim2); CHKERRQ(ierr); - dim2 *= dim2; + dim2 *= dim2; for(i=0;iK[dim2*(c-cStart)+i] = K[i]; PetscFunctionReturn(0); } @@ -401,7 +401,7 @@ PetscErrorCode TDySoilDensityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *d PetscFunctionBegin; TDyOptions *options = &tdy->options; - *den = options->default_soil_density; + *den = options->soil_density; PetscFunctionReturn(0); } @@ -410,7 +410,7 @@ PetscErrorCode TDySoilSpecificHeatFunctionDefault(TDy tdy, PetscReal *x, PetscRe PetscFunctionBegin; TDyOptions *options = &tdy->options; - *cr = options->default_soil_specific_heat; + *cr = options->soil_specific_heat; PetscFunctionReturn(0); } @@ -427,7 +427,7 @@ PetscErrorCode TDyPermeabilityFunctionDefault(TDy tdy, PetscReal *x, PetscReal * } } for (int i=0; idefault_permeability; + K[i*dim+i] = options->permeability; return 0; PetscFunctionReturn(0); } @@ -445,7 +445,7 @@ PetscErrorCode TDyThermalConductivityFunctionDefault(TDy tdy, PetscReal *x, Pets } } for (int i=0; idefault_thermal_conductivity; + K[i*dim+i] = options->thermal_conductivity; return 0; PetscFunctionReturn(0); } @@ -453,6 +453,6 @@ PetscErrorCode TDyThermalConductivityFunctionDefault(TDy tdy, PetscReal *x, Pets PetscErrorCode TDyPorosityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *por, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; - *por = options->default_porosity; + *por = options->porosity; PetscFunctionReturn(0); } From ff62b36b5893a7895dce1f3623cf95209edd7cb7 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 9 Aug 2021 13:44:05 -0700 Subject: [PATCH 03/10] Renamed TDyXYZFunctionDefault -> TDyConstantXYZFunction (material props) --- include/private/tdymaterialpropertiesimpl.h | 10 +++++----- src/tdydriver.c | 14 +++++++------- src/tdymaterialproperties.c | 12 ++++++------ 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/include/private/tdymaterialpropertiesimpl.h b/include/private/tdymaterialpropertiesimpl.h index a5fe8042..76918784 100644 --- a/include/private/tdymaterialpropertiesimpl.h +++ b/include/private/tdymaterialpropertiesimpl.h @@ -26,11 +26,11 @@ PETSC_INTERN PetscBool TDyIsThermalConductivytSet(TDy); PETSC_INTERN PetscBool TDyIsSoilSpecificHeatSet(TDy); PETSC_INTERN PetscBool TDyIsSoilDensitySet(TDy); -PETSC_INTERN PetscErrorCode TDySoilDensityFunctionDefault(TDy,PetscReal*,PetscReal*,void*); -PETSC_INTERN PetscErrorCode TDySoilSpecificHeatFunctionDefault(TDy,PetscReal*,PetscReal*,void*); -PETSC_INTERN PetscErrorCode TDyPermeabilityFunctionDefault(TDy,PetscReal *,PetscReal *,void*); -PETSC_INTERN PetscErrorCode TDyThermalConductivityFunctionDefault(TDy,PetscReal *,PetscReal *,void*); -PETSC_INTERN PetscErrorCode TDyPorosityFunctionDefault(TDy,PetscReal *,PetscReal *,void*); +PETSC_INTERN PetscErrorCode TDyConstantSoilDensityFunction(TDy,PetscReal*,PetscReal*,void*); +PETSC_INTERN PetscErrorCode TDyConstantSoilSpecificHeatFunction(TDy,PetscReal*,PetscReal*,void*); +PETSC_INTERN PetscErrorCode TDyConstantPermeabilityFunction(TDy,PetscReal *,PetscReal *,void*); +PETSC_INTERN PetscErrorCode TDyConstantThermalConductivityFunction(TDy,PetscReal *,PetscReal *,void*); +PETSC_INTERN PetscErrorCode TDyConstantPorosityFunction(TDy,PetscReal *,PetscReal *,void*); #endif diff --git a/src/tdydriver.c b/src/tdydriver.c index 2de718b9..5f6e7214 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -34,26 +34,26 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { } if (!TDyIsPorositySet(tdy)) { - ierr = TDySetPorosityFunction(tdy,TDyPorosityFunctionDefault,PETSC_NULL); CHKERRQ(ierr); + ierr = TDySetPorosityFunction(tdy,TDyConstantPorosityFunction,PETSC_NULL); CHKERRQ(ierr); } if (!TDyIsPermeabilitySet(tdy)){ - ierr = TDySetPermeabilityFunction(tdy,TDyPermeabilityFunctionDefault,PETSC_NULL); CHKERRQ(ierr); + ierr = TDySetPermeabilityFunction(tdy,TDyConstantPermeabilityFunction,PETSC_NULL); CHKERRQ(ierr); } if (tdy->mode == TH) { if (!TDyIsThermalConductivytSet(tdy)) { - ierr = TDySetThermalConductivityFunction(tdy,TDyThermalConductivityFunctionDefault,PETSC_NULL); CHKERRQ(ierr); + ierr = TDySetThermalConductivityFunction(tdy,TDyConstantThermalConductivityFunction,PETSC_NULL); CHKERRQ(ierr); } if (!TDyIsSoilDensitySet(tdy)) { - //ierr = TDySetSoilDensity(tdy,TDySoilDensityFunctionDefault); CHKERRQ(ierr); - ierr = TDySetSoilDensityFunction(tdy,TDySoilDensityFunctionDefault,PETSC_NULL); CHKERRQ(ierr); + //ierr = TDySetSoilDensity(tdy,TDyConstantSoilDensityFunction); CHKERRQ(ierr); + ierr = TDySetSoilDensityFunction(tdy,TDyConstantSoilDensityFunction,PETSC_NULL); CHKERRQ(ierr); } if (!TDyIsSoilSpecificHeatSet(tdy)) { - ierr = TDySetSoilSpecificHeatFunction(tdy,TDySoilSpecificHeatFunctionDefault,PETSC_NULL); CHKERRQ(ierr); + ierr = TDySetSoilSpecificHeatFunction(tdy,TDyConstantSoilSpecificHeatFunction,PETSC_NULL); CHKERRQ(ierr); } } @@ -88,7 +88,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unrecognized time integration method."); } - // create time integrator + // create time integrator switch(tdy->ti->time_integration_method) { case TDySNES: ierr = SNESCreate(PETSC_COMM_WORLD,&snes); diff --git a/src/tdymaterialproperties.c b/src/tdymaterialproperties.c index bcaa554e..d6a410bc 100644 --- a/src/tdymaterialproperties.c +++ b/src/tdymaterialproperties.c @@ -394,10 +394,10 @@ PetscErrorCode TDyGetPorosityValuesLocal(TDy tdy, PetscInt *ni, PetscScalar y[]) } /* - Default function for material properties + Constant functions for material properties */ -PetscErrorCode TDySoilDensityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *den, void *ctx) { +PetscErrorCode TDyConstantSoilDensityFunction(TDy tdy, PetscReal *x, PetscReal *den, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; @@ -406,7 +406,7 @@ PetscErrorCode TDySoilDensityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *d PetscFunctionReturn(0); } -PetscErrorCode TDySoilSpecificHeatFunctionDefault(TDy tdy, PetscReal *x, PetscReal *cr, void *ctx) { +PetscErrorCode TDyConstantSoilSpecificHeatFunction(TDy tdy, PetscReal *x, PetscReal *cr, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; @@ -414,7 +414,7 @@ PetscErrorCode TDySoilSpecificHeatFunctionDefault(TDy tdy, PetscReal *x, PetscRe PetscFunctionReturn(0); } -PetscErrorCode TDyPermeabilityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *K, void *ctx) { +PetscErrorCode TDyConstantPermeabilityFunction(TDy tdy, PetscReal *x, PetscReal *K, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; PetscInt dim; @@ -432,7 +432,7 @@ PetscErrorCode TDyPermeabilityFunctionDefault(TDy tdy, PetscReal *x, PetscReal * PetscFunctionReturn(0); } -PetscErrorCode TDyThermalConductivityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *K, void *ctx) { +PetscErrorCode TDyConstantThermalConductivityFunction(TDy tdy, PetscReal *x, PetscReal *K, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; PetscErrorCode ierr; @@ -450,7 +450,7 @@ PetscErrorCode TDyThermalConductivityFunctionDefault(TDy tdy, PetscReal *x, Pets PetscFunctionReturn(0); } -PetscErrorCode TDyPorosityFunctionDefault(TDy tdy, PetscReal *x, PetscReal *por, void *ctx) { +PetscErrorCode TDyConstantPorosityFunction(TDy tdy, PetscReal *x, PetscReal *por, void *ctx) { PetscFunctionBegin; TDyOptions *options = &tdy->options; *por = options->porosity; From 82b7210f643f793e18b67cb619ab5b7eb793bdf0 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 9 Aug 2021 17:10:10 -0700 Subject: [PATCH 04/10] Reorganized options into TDyOptions container and added some new ones. --- include/private/tdycoreimpl.h | 21 +--- include/private/tdyoptions.h | 52 +++++--- src/mpfao/3D/tdympfao3D_core.c | 36 +++--- src/mpfao/3D/tdympfao3D_ts.c | 28 ++--- src/mpfao/3D/tdympfao3D_utils.c | 2 +- src/mpfao/tdympfao.c | 47 +++++--- src/tdybdm.c | 2 +- src/tdycore.c | 202 ++++++++++++++++---------------- src/tdydriver.c | 10 +- src/tdyregression.c | 20 ++-- src/tdyrichards.c | 6 +- src/tdyth.c | 2 +- src/tdyti.c | 2 +- src/tdytimers.c | 4 +- src/tdytpf.c | 2 +- 15 files changed, 226 insertions(+), 210 deletions(-) diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index a399b200..2b70ffec 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -43,10 +43,8 @@ struct _p_TDy { TDyTimeIntegrator ti; TDySetupFlags setupflags; TDyIO io; - PetscBool init_with_random_field; - PetscBool init_from_file; - char init_file[PETSC_MAX_PATH_LEN]; + // options that determine the behavior(s) of the dycore TDyOptions options; /* arrays of the size of the Hasse diagram */ @@ -56,9 +54,6 @@ struct _p_TDy { PetscInt ncv,nfv; /* number of {cell|face} vertices */ /* non-linear function of liquid pressure */ - PetscInt rho_type; - PetscInt mu_type; - PetscInt enthalpy_type; PetscReal *rho, *drho_dP, *d2rho_dP2; /* density of water [kg m-3]*/ PetscReal *vis, *dvis_dP, *d2vis_dP2; /* viscosity of water [Pa s] */ PetscReal *h, *dh_dP, *dh_dT; /* enthalpy of water */ @@ -101,12 +96,6 @@ struct _p_TDy { void *soildensityctx; void *soilspecificheatctx; - /* method-specific information*/ - TDyMethod method; - TDyMode mode; - TDyQuadratureType qtype; - PetscBool allow_unsuitable_mesh; - /* Wheeler-Yotov */ PetscInt *vmap; /* [cell,local_vertex] --> global_vertex */ PetscInt *emap; /* [cell,local_vertex,direction] --> global_face */ @@ -122,9 +111,6 @@ struct _p_TDy { PetscInt *faces; /* MPFA-O */ - PetscInt mpfao_gmatrix_method; - PetscInt mpfao_bc_type; - TDyMesh *mesh; PetscReal ****subc_Gmatrix; /* Gmatrix for subcells */ PetscReal ***Trans; @@ -149,12 +135,7 @@ struct _p_TDy { PetscInt *closureSize, **closure, maxClosureSize; - PetscBool output_mesh; - PetscBool regression_testing; TDyRegression *regression; - - /* Timers enabled? */ - PetscBool enable_timers; }; diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index 931306e2..a77b324a 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -3,22 +3,44 @@ #include -typedef struct{ - PetscReal gravity_constant; - - // Material properties - PetscReal porosity; - PetscReal permeability; - PetscReal soil_density; - PetscReal soil_specific_heat; - PetscReal thermal_conductivity; - - // Characteristic curves - PetscReal residual_saturation; - PetscReal gardner_n; - PetscReal vangenuchten_m; - PetscReal vangenuchten_alpha; +typedef struct { + // Model settings + PetscReal gravity_constant; + TDyMode mode; + TDyWaterDensityType rho_type; + PetscInt mu_type; + PetscInt enthalpy_type; + + // Numerics settings + TDyMethod method; + TDyQuadratureType qtype; + PetscInt mpfao_gmatrix_method; + PetscInt mpfao_bc_type; + PetscBool tpf_allow_all_meshes; + + // Material properties + PetscReal porosity; + PetscReal permeability; + PetscReal soil_density; + PetscReal soil_specific_heat; + PetscReal thermal_conductivity; + + // Characteristic curves + PetscReal residual_saturation; + PetscReal gardner_n; + PetscReal vangenuchten_m; + PetscReal vangenuchten_alpha; + + // I/O settings + PetscBool init_with_random_field; + PetscBool init_from_file; + char init_file[PETSC_MAX_PATH_LEN]; + PetscBool output_mesh; + PetscBool regression_testing; + + // Timers enabled? + PetscBool enable_timers; } TDyOptions; #endif diff --git a/src/mpfao/3D/tdympfao3D_core.c b/src/mpfao/3D/tdympfao3D_core.c index 6be4e7b4..6227de07 100644 --- a/src/mpfao/3D/tdympfao3D_core.c +++ b/src/mpfao/3D/tdympfao3D_core.c @@ -69,7 +69,7 @@ PetscErrorCode TDyComputeGMatrixMPFAOFor3DMesh(TDy tdy) { // extract thermal conductivity tensor PetscReal Kappa[3][3]; - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { for (PetscInt ii=0; iiKappa0[icell*dim*dim + ii*dim + jj]; @@ -102,7 +102,7 @@ PetscErrorCode TDyComputeGMatrixMPFAOFor3DMesh(TDy tdy) { ierr = TDyComputeEntryOfGMatrix3D(area, normal, K, nu, subcells->T[subcell_id], dim, &(tdy->subc_Gmatrix[icell][isubcell][ii][jj])); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = TDySubCell_GetIthNuVector(subcells, subcell_id, jj, dim, &nu[0]); CHKERRQ(ierr); ierr = TDyComputeEntryOfGMatrix3D(area, normal, Kappa, @@ -151,7 +151,7 @@ PetscErrorCode TDyComputeGMatrixTPFFor3DMesh(TDy tdy) { // extract thermal conductivity tensor PetscReal Kappa[3][3]; - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { for (ii=0; iiKappa0[icell*dim*dim + ii*dim + jj]; @@ -243,7 +243,7 @@ PetscErrorCode TDyComputeGMatrixTPFFor3DMesh(TDy tdy) { tdy->subc_Gmatrix[icell][isubcell][ii][jj] = area * (dot_prod) * K_aveg/(dist); } - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { if (ii == jj) { ierr = TDySubCell_GetIthNuStarVector(subcells, subcell_id, jj, dim, &nu[0]); CHKERRQ(ierr); @@ -563,7 +563,8 @@ PetscErrorCode ComputeTransmissibilityMatrix_ForNonCornerVertex(TDy tdy, PetscInt npitf_dir_bc_all, npitf_neu_bc_all; - if (tdy->mpfao_bc_type == MPFAO_DIRICHLET_BC || tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC) { + if (tdy->options.mpfao_bc_type == MPFAO_DIRICHLET_BC || + tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC) { nflux_dir_bc_up = nflux_all_bc_up; nflux_dir_bc_dn = nflux_all_bc_dn; npitf_dir_bc_all= npitf_bc_all; @@ -1020,7 +1021,7 @@ PetscErrorCode ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInte // Need to swap entries in (*Trans)[vertices->id][:][:] such that // Step-1. Rows correspond to faces saved in the order of vertices->face_ids[:] // Step-2. Columns are such that boundary cells are followed by the internal cell. - // The boundary cells are in the order of vertices->cell_ids[1:3] + // The boundary cells are in the order of vertices->cell_ids[1:3] // Note: vertices->cell_ids[0] correspond to internal cell // Step-3. Finally move the last column as the first column. So, columns // are in the order of vertices->cell_ids[0:3] @@ -1137,7 +1138,7 @@ PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() + TDY_START_FUNCTION_TIMER() // If a face is shared by two cells that belong to different // regions, zero the rows in the transmissiblity matrix @@ -1155,7 +1156,7 @@ PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy tdy) { PetscInt row[1]; row[0] = iface*num_subfaces + isubface; ierr = MatZeroRows(tdy->Trans_mat,1,row,0.0,0,0); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = MatZeroRows(tdy->Temp_Trans_mat,1,row,0.0,0,0); CHKERRQ(ierr); } } @@ -1180,22 +1181,23 @@ PetscErrorCode TDyComputeTransmissibilityMatrix3DMesh(TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() - + for (ivertex=0; ivertexnum_vertices; ivertex++) { if (!vertices->is_local[ivertex]) continue; if (vertices->num_boundary_faces[ivertex] == 0 || vertices->num_internal_cells[ivertex] > 1) { ierr = ComputeTransmissibilityMatrix_ForNonCornerVertex(tdy, ivertex, cells, 0); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = ComputeTransmissibilityMatrix_ForNonCornerVertex(tdy, ivertex, cells, 1); CHKERRQ(ierr); } } else { // It is assumed that neumann boundary condition is a zero-flux boundary condition. // Thus, compute transmissiblity entries only for dirichlet boundary condition. - if (tdy->mpfao_bc_type == MPFAO_DIRICHLET_BC || tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC) { + if (tdy->options.mpfao_bc_type == MPFAO_DIRICHLET_BC || + tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC) { ierr = ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInternalVertices(tdy, ivertex, cells, 0); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = ComputeTransmissibilityMatrix_ForBoundaryVertex_NotSharedWithInternalVertices(tdy, ivertex, cells, 1); CHKERRQ(ierr); } } @@ -1205,14 +1207,14 @@ PetscErrorCode TDyComputeTransmissibilityMatrix3DMesh(TDy tdy) { ierr = MatAssemblyBegin(tdy->Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(tdy->Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = MatAssemblyBegin(tdy->Temp_Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(tdy->Temp_Trans_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } TDyRegion *region = &mesh->region_connected; if (region->num_cells > 0){ - if (tdy->mpfao_gmatrix_method == MPFAO_GMATRIX_TPF ) { + if (tdy->options.mpfao_gmatrix_method == MPFAO_GMATRIX_TPF ) { ierr = TDyUpdateTransmissibilityMatrix(tdy); CHKERRQ(ierr); } else { PetscPrintf(PETSC_COMM_WORLD,"WARNING -- Connected region option is only supported with MPFA-O TPF\n"); @@ -1437,11 +1439,11 @@ PetscErrorCode ComputeFacePermeabilityTensor(TDy tdy, PetscInt face_id, PetscRea /// g = gravity vector /// GravDis = gravity discretization term that is not dependent on the unknown variable(s) /// such as pressure, temperautre. Thus, this term is precomputed. -/// +/// /// Starnoni, M., Berre, I., Keilegavlen, E., & Nordbotten, J. M. (2019). /// Consistent mpfa discretization for flow in the presence of gravity. Water /// Resources Research, 55, 10105– 10118. https://doi.org/10.1029/2019WR025384 -/// +/// /// @param [in] tdy A TDy struct /// @returns 0 on success, or a non-zero error code on failure PetscErrorCode TDyComputeGravityDiscretizationFor3DMesh(TDy tdy) { @@ -1492,7 +1494,7 @@ PetscErrorCode TDyComputeGravityDiscretizationFor3DMesh(TDy tdy) { // Currently, only zero-flux neumann boundary condition is implemented. // If the boundary condition is neumann, then gravity discretization term is zero - if (tdy->mpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up < 0 || cell_id_dn < 0)) continue; + if (tdy->options.mpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up < 0 || cell_id_dn < 0)) continue; PetscInt cell_id; if (cell_id_up < 0) { diff --git a/src/mpfao/3D/tdympfao3D_ts.c b/src/mpfao/3D/tdympfao3D_ts.c index 9b9d0701..f7e3176c 100644 --- a/src/mpfao/3D/tdympfao3D_ts.c +++ b/src/mpfao/3D/tdympfao3D_ts.c @@ -49,7 +49,7 @@ PetscErrorCode TDyMPFAOIFunction_Vertices_3DMesh(Vec Ul, Vec R, void *ctx) { // Compute = T*P PetscScalar TtimesP[nflux_in + npitf_bc]; - + for (PetscInt irow=0; irowvis_BND[-cell_id_up-1]; ukvr = Kr/vis; - if (tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_up-1] <= tdy->Pref) { + if (tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_up-1] <= tdy->Pref) { set_flow_to_zero = PETSC_TRUE; } } @@ -120,7 +120,7 @@ PetscErrorCode TDyMPFAOIFunction_Vertices_3DMesh(Vec Ul, Vec R, void *ctx) { PetscReal vis = tdy->vis_BND[-cell_id_dn-1]; ukvr = Kr/vis; - if (tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_dn-1] <= tdy->Pref) { + if (tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_dn-1] <= tdy->Pref) { set_flow_to_zero = PETSC_TRUE; } } @@ -211,7 +211,7 @@ PetscErrorCode TDyMPFAOIFunction_3DMesh(TS ts,PetscReal t,Vec U,Vec U_t,Vec R,vo PetscReal S = cc->S[icell]; PetscReal dS_dP = cc->dS_dP[icell]; - PetscReal dmass_dP = + PetscReal dmass_dP = rho * dporosity_dP * S + drho_dP * porosity * S + rho * porosity * dS_dP; @@ -232,7 +232,7 @@ PetscErrorCode TDyMPFAOIFunction_3DMesh(TS ts,PetscReal t,Vec U,Vec U_t,Vec R,vo ierr = VecView(R,viewer); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); #endif - + TDY_STOP_FUNCTION_TIMER() PetscFunctionReturn(0); @@ -281,7 +281,7 @@ PetscErrorCode TDyMPFAOIJacobian_Vertices_3DMesh(Vec Ul, Mat A, void *ctx) { PetscInt nflux_in = vertices->num_faces[ivertex] - vertices->num_boundary_faces[ivertex]; // Compute T*P - PetscScalar TtimesP[nflux_in + npitf_bc]; + PetscScalar TtimesP[nflux_in + npitf_bc]; for (PetscInt irow=0; irow < nflux_in + npitf_bc; irow++) { PetscInt face_id = face_ids[irow]; @@ -308,7 +308,7 @@ PetscErrorCode TDyMPFAOIJacobian_Vertices_3DMesh(Vec Ul, Mat A, void *ctx) { // d(fluxm_ij)/dP_k = // rho_ij + (kr/mu)_{ij,upwind} * T_k // - + for (PetscInt irow=0; irowmpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up<0 || cell_id_dn <0)) continue; + if ( tdy->options.mpfao_bc_type == MPFAO_NEUMANN_BC && (cell_id_up<0 || cell_id_dn <0)) continue; PetscReal dukvr_dPup = 0.0; PetscReal dukvr_dPdn = 0.0; @@ -369,8 +369,8 @@ PetscErrorCode TDyMPFAOIJacobian_Vertices_3DMesh(Vec Ul, Mat A, void *ctx) { ukvr = Kr/vis; dukvr_dPup = 0.0; - - if (tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_up-1] <= tdy->Pref) { + + if (tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_up-1] <= tdy->Pref) { set_jac_to_zero = PETSC_TRUE; } } @@ -396,9 +396,9 @@ PetscErrorCode TDyMPFAOIJacobian_Vertices_3DMesh(Vec Ul, Mat A, void *ctx) { ukvr = Kr/vis; dukvr_dPdn = 0.0; - if (tdy->mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_dn-1] <= tdy->Pref) { - set_jac_to_zero = PETSC_TRUE; - } + if (tdy->options.mpfao_bc_type == MPFAO_SEEPAGE_BC && tdy->P_BND[-cell_id_dn-1] <= tdy->Pref) { + set_jac_to_zero = PETSC_TRUE; + } } } @@ -571,7 +571,7 @@ PetscErrorCode TDyMPFAOIJacobian_3DMesh(TS ts,PetscReal t,Vec U,Vec U_t,PetscRea ierr = MatView(A,viewer); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); #endif - + TDY_STOP_FUNCTION_TIMER() PetscFunctionReturn(0); diff --git a/src/mpfao/3D/tdympfao3D_utils.c b/src/mpfao/3D/tdympfao3D_utils.c index e920aafb..28f5b034 100644 --- a/src/mpfao/3D/tdympfao3D_utils.c +++ b/src/mpfao/3D/tdympfao3D_utils.c @@ -707,7 +707,7 @@ PetscErrorCode TDyMPFAO_SetBoundaryPressure(TDy tdy, Vec Ul) { ierr = VecGetArray(Ul,&u_p); CHKERRQ(ierr); ierr = VecGetArray(tdy->P_vec,&p_vec_ptr); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { for (c=0;cmpfao_gmatrix_method) { + switch (tdy->options.mpfao_gmatrix_method) { case MPFAO_GMATRIX_DEFAULT: ierr = TDyComputeGMatrixMPFAOFor3DMesh(tdy); CHKERRQ(ierr); break; @@ -295,8 +295,8 @@ PetscErrorCode TDyMPFAO_AllocateMemoryForBoundaryValues(TDy tdy) { PetscInt i; PetscReal dden_dP, d2den_dP2, dmu_dP, d2mu_dP2; for (i=0;iPref, tdy->rho_type, &(tdy->rho_BND[i]), &dden_dP, &d2den_dP2); CHKERRQ(ierr); - ierr = ComputeWaterViscosity(tdy->Pref, tdy->mu_type, &(tdy->vis_BND[i]), &dmu_dP, &d2mu_dP2); CHKERRQ(ierr); + ierr = ComputeWaterDensity(tdy->Pref, tdy->options.rho_type, &(tdy->rho_BND[i]), &dden_dP, &d2den_dP2); CHKERRQ(ierr); + ierr = ComputeWaterViscosity(tdy->Pref, tdy->options.mu_type, &(tdy->vis_BND[i]), &dmu_dP, &d2mu_dP2); CHKERRQ(ierr); } TDY_STOP_FUNCTION_TIMER() @@ -321,7 +321,8 @@ PetscErrorCode TDyMPFAO_AllocateMemoryForEnergyBoundaryValues(TDy tdy) { PetscInt i; PetscReal dh_dP, dh_dT; for (i=0;iTref, tdy->Pref,tdy->enthalpy_type, &(tdy->h_BND[i]), &dh_dP, &dh_dT); CHKERRQ(ierr); + ierr = ComputeWaterEnthalpy(tdy->Tref, tdy->Pref,tdy->options.enthalpy_type, + &(tdy->h_BND[i]), &dh_dP, &dh_dT); CHKERRQ(ierr); } TDY_STOP_FUNCTION_TIMER() @@ -410,8 +411,10 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { break; case 3: ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); - if (tdy->mode == TH){ierr = TDyAllocate_RealArray_3D(&tdy->Temp_Trans, - tdy->mesh->num_vertices, tdy->nfv, tdy->nfv); 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 = 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); @@ -430,19 +433,21 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { } ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, - nsubcells, nrow, ncol); CHKERRQ(ierr); - if (tdy->mode == TH) {ierr = TDyAllocate_RealArray_4D(&tdy->Temp_subc_Gmatrix, tdy->mesh->num_cells, - nsubcells, nrow, ncol); CHKERRQ(ierr);} + 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); + } - /* Setup 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->method == MPFA_O_DAE); + use_dae = (tdy->options.method == MPFA_O_DAE); ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); if (!use_dae) { - if (tdy->mode == TH){ + 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); @@ -454,7 +459,9 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); } } else { - if (tdy->mode == TH) {SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_DAE");} + 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); @@ -467,12 +474,14 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { for(p=pStart; pmode == TH){ + if (tdy->options.mode == TH) { ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); } if (use_dae) { - if (tdy->mode == TH) {SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_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); } @@ -506,7 +515,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); ierr = VecZeroEntries(tdy->GravDisVec); - if (tdy->mode == TH){ + 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); @@ -516,7 +525,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} - if (tdy->mode == TH){ + if (tdy->options.mode == TH){ if (tdy->ops->computethermalconductivity) {ierr = SetThermalConductivityFromFunction(tdy); CHKERRQ(ierr);} if (tdy->ops->computesoildensity) { ierr = SetSoilDensityFromFunction(tdy); CHKERRQ(ierr); @@ -535,7 +544,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { if (dim == 3) { ierr = TDyMPFAO_AllocateMemoryForBoundaryValues(tdy); CHKERRQ(ierr); ierr = TDyMPFAO_AllocateMemoryForSourceSinkValues(tdy); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = TDyMPFAO_AllocateMemoryForEnergyBoundaryValues(tdy); CHKERRQ(ierr); ierr = TDyMPFAO_AllocateMemoryForEnergySourceSinkValues(tdy); CHKERRQ(ierr); } @@ -657,7 +666,7 @@ PetscErrorCode TDyMPFAOComputeSystem(TDy tdy,Mat K,Vec F) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() - + ierr = MatZeroEntries(K); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); ierr = MatAssemblyEnd (K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); diff --git a/src/tdybdm.c b/src/tdybdm.c index 2b9ab02a..7b865a7f 100644 --- a/src/tdybdm.c +++ b/src/tdybdm.c @@ -430,7 +430,7 @@ PetscErrorCode TDyBDMComputeSystem(TDy tdy,Mat K,Vec F) { nlocal = dim*nv + 1; /* Get quadrature */ - switch(tdy->qtype) { + switch(tdy->options.qtype) { case FULL: ierr = PetscDTGaussTensorQuadrature(dim ,1,nq1d,-1,+1,& quadrature); CHKERRQ(ierr); ierr = PetscDTGaussTensorQuadrature(dim-1,1,nq1d,-1,+1,&face_quadrature); CHKERRQ(ierr); diff --git a/src/tdycore.c b/src/tdycore.c index 00943cf2..196c9284 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -181,11 +181,22 @@ PetscErrorCode TDyFinalize() { PetscFunctionReturn(0); } -PetscErrorCode TDySetOptionDefaults(TDyOptions *options){ +static PetscErrorCode SetDefaultOptions(TDy tdy) { PetscFunctionBegin; - PetscErrorCode ierr; + TDyOptions *options = &tdy->options; + options->mode = RICHARDS; + options->method = WY; options->gravity_constant = 9.8068; + options->rho_type = WATER_DENSITY_CONSTANT; + options->mu_type = WATER_VISCOSITY_CONSTANT; + options->enthalpy_type = WATER_ENTHALPY_CONSTANT; + + options->mpfao_gmatrix_method = MPFAO_GMATRIX_DEFAULT; + options->mpfao_bc_type = MPFAO_DIRICHLET_BC; + options->tpf_allow_all_meshes = PETSC_FALSE; + + options->qtype = FULL; options->porosity=0.25; options->permeability=1.e-12; @@ -198,10 +209,12 @@ PetscErrorCode TDySetOptionDefaults(TDyOptions *options){ options->vangenuchten_m=0.8; options->vangenuchten_alpha=1.e-4; + options->init_with_random_field = PETSC_FALSE; + options->init_from_file = PETSC_FALSE; + PetscFunctionReturn(0); } - PetscErrorCode TDyCreate(TDy *_tdy) { TDy tdy; PetscErrorCode ierr; @@ -217,7 +230,7 @@ PetscErrorCode TDyCreate(TDy *_tdy) { *_tdy = tdy; tdy->setupflags |= TDyCreated; - ierr = TDySetOptionDefaults(&tdy->options); + SetDefaultOptions(tdy); ierr = TDyIOCreate(&tdy->io); CHKERRQ(ierr); @@ -225,20 +238,11 @@ PetscErrorCode TDyCreate(TDy *_tdy) { tdy->Pref = 101325; tdy->Tref = 25; tdy->gravity[0] = 0; tdy->gravity[1] = 0; tdy->gravity[2] = 0; - tdy->rho_type = WATER_DENSITY_CONSTANT; - tdy->mu_type = WATER_VISCOSITY_CONSTANT; - tdy->enthalpy_type = WATER_ENTHALPY_CONSTANT; - tdy->mpfao_gmatrix_method = MPFAO_GMATRIX_DEFAULT; - tdy->mpfao_bc_type = MPFAO_DIRICHLET_BC; - tdy->allow_unsuitable_mesh = PETSC_FALSE; - tdy->init_with_random_field = PETSC_FALSE; - tdy->init_from_file = PETSC_FALSE; /* initialize method information to null */ tdy->vmap = NULL; tdy->emap = NULL; tdy->Alocal = NULL; tdy->Flocal = NULL; tdy->quad = NULL; tdy->faces = NULL; tdy->LtoG = NULL; tdy->orient = NULL; - tdy->qtype = FULL; tdy->setupflags |= TDyParametersInitialized; PetscFunctionReturn(0); @@ -527,7 +531,13 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetupNumericalMethods()"); } + // Collect options from command line arguments. TDyOptions *options = &tdy->options; + + //------------------------------------------ + // Set options using command line parameters + //------------------------------------------ + PetscValidHeaderSpecific(tdy,TDY_CLASSID,1); ierr = PetscObjectOptionsBegin((PetscObject)tdy); CHKERRQ(ierr); PetscBool flag; @@ -551,92 +561,88 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { // Model options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Model options",""); CHKERRQ(ierr); + ierr = PetscOptionsEnum("-tdy_mode","Flow mode", + "TDySetMode",TDyModes,(PetscEnum)options->mode, + (PetscEnum *)&options->mode, NULL); CHKERRQ(ierr); ierr = PetscOptionsReal("-tdy_gravity", "Magnitude of gravity vector", NULL, options->gravity_constant, &options->gravity_constant, NULL); CHKERRQ(ierr); - TDyWaterDensityType densitytype = WATER_DENSITY_CONSTANT; ierr = PetscOptionsEnum("-tdy_water_density","Water density vertical profile", "TDySetWaterDensityType",TDyWaterDensityTypes, - (PetscEnum)densitytype,(PetscEnum *)&densitytype, - &flag); CHKERRQ(ierr); + (PetscEnum)options->rho_type, + (PetscEnum *)&options->rho_type, NULL); CHKERRQ(ierr); + + // Named boundary condition functions + char func_name[PETSC_MAX_PATH_LEN]; + ierr = PetscOptionsGetString(NULL, NULL, "-tdy_pressure_bc", func_name, + sizeof(func_name), &flag); CHKERRQ(ierr); if (flag) { - ierr = TDySetWaterDensityType(tdy,densitytype); CHKERRQ(ierr); + // TODO: When/where are we supposed to get the context for this function? + ierr = TDySelectBoundaryPressureFn(tdy, func_name, NULL); } - TDyMode mode = RICHARDS; - ierr = PetscOptionsEnum("-tdy_mode","Flow mode", - "TDySetMode",TDyModes,(PetscEnum)mode,(PetscEnum *)&mode, - &flag); CHKERRQ(ierr); - - if (flag && (mode != tdy->mode)) { - ierr = TDySetMode(tdy,mode); CHKERRQ(ierr); + ierr = PetscOptionsGetString(NULL, NULL, "-tdy_velocity_bc", func_name, + sizeof(func_name), &flag); CHKERRQ(ierr); + if (flag) { + // TODO: When/where are we supposed to get the context for this function? + ierr = TDySelectBoundaryVelocityFn(tdy, func_name, NULL); } - ierr = PetscOptionsReal("-tdy_pressure_bc", - "Assigns a named pressure function to the domain boundary", - NULL, options->gravity_constant, &options->gravity_constant, - NULL); CHKERRQ(ierr); + ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); - TDyMethod method = WY; ierr = PetscOptionsEnum("-tdy_method","Discretization method", "TDySetDiscretizationMethod",TDyMethods, - (PetscEnum)method,(PetscEnum *)&method, &flag); CHKERRQ(ierr); - if (flag && (method != tdy->method)) { - ierr = TDySetDiscretizationMethod(tdy,method); CHKERRQ(ierr); - } + (PetscEnum)options->method,(PetscEnum *)&options->method, + NULL); CHKERRQ(ierr); TDyQuadratureType qtype = FULL; ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods", - "TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&qtype, - &flag); CHKERRQ(ierr); - if (flag && (qtype != tdy->qtype)) { - ierr = TDySetQuadratureType(tdy,qtype); CHKERRQ(ierr); - } + "TDySetQuadratureType",TDyQuadratureTypes, + (PetscEnum)qtype,(PetscEnum *)&options->qtype, + NULL); CHKERRQ(ierr); ierr = PetscOptionsBool("-tdy_tpf_allow_all_meshes", "Enable to allow non-orthgonal meshes in finite volume TPF method", - "",tdy->allow_unsuitable_mesh, - &(tdy->allow_unsuitable_mesh),NULL); CHKERRQ(ierr); - TDyMPFAOGmatrixMethod gmatrixmethod = MPFAO_GMATRIX_DEFAULT; + "",options->tpf_allow_all_meshes, + &(options->tpf_allow_all_meshes),NULL); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_mpfao_gmatrix_method","MPFA-O gmatrix method", - "TDySetMPFAOGmatrixMethod",TDyMPFAOGmatrixMethods,(PetscEnum)gmatrixmethod,(PetscEnum *)&gmatrixmethod, - &flag); CHKERRQ(ierr); - if (flag && (gmatrixmethod != tdy->mpfao_gmatrix_method)) { - ierr = TDySetMPFAOGmatrixMethod(tdy,gmatrixmethod); CHKERRQ(ierr); - } + "TDySetMPFAOGmatrixMethod",TDyMPFAOGmatrixMethods, + (PetscEnum)options->mpfao_gmatrix_method, + (PetscEnum *)&options->mpfao_gmatrix_method, + NULL); CHKERRQ(ierr); TDyMPFAOBoundaryConditionType bctype = MPFAO_DIRICHLET_BC; ierr = PetscOptionsEnum("-tdy_mpfao_boundary_condition_type","MPFA-O boundary condition type", - "TDySetMPFAOBoundaryConditionType",TDyMPFAOBoundaryConditionTypes,(PetscEnum)bctype,(PetscEnum *)&bctype, - &flag); CHKERRQ(ierr); - if (flag && (bctype != tdy->mpfao_bc_type)) { + "TDySetMPFAOBoundaryConditionType",TDyMPFAOBoundaryConditionTypes,(PetscEnum)bctype, + (PetscEnum *)&bctype, &flag); CHKERRQ(ierr); + if (flag && (bctype != tdy->options.mpfao_bc_type)) { ierr = TDySetMPFAOBoundaryConditionType(tdy,bctype); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); ierr = PetscOptionsBool("-tdy_init_with_random_field", "Initialize solution with a random field","", - tdy->init_with_random_field, - &(tdy->init_with_random_field),NULL); CHKERRQ(ierr); - - ierr = PetscOptionsGetString(NULL,NULL,"-tdy_init_file", tdy->init_file, - sizeof(tdy->init_file), - &tdy->init_from_file); CHKERRQ(ierr); + options->init_with_random_field, + &(options->init_with_random_field),NULL); CHKERRQ(ierr); + ierr = PetscOptionsGetString(NULL,NULL,"-tdy_init_file", options->init_file, + sizeof(options->init_file), + &options->init_from_file); CHKERRQ(ierr); + + if (options->init_from_file && options->init_with_random_field) { + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + "Only one of -tdy_init_from_file and -tdy_init_with_random_field can be specified"); + } ierr = PetscOptionsBool("-tdy_regression_test", - "Enable output of a regression file","",tdy->regression_testing, - &(tdy->regression_testing),NULL); CHKERRQ(ierr); + "Enable output of a regression file","", + options->regression_testing, + &(options->regression_testing),NULL); CHKERRQ(ierr); ierr = PetscOptionsBool("-tdy_output_mesh", - "Enable output of mesh attributes","",tdy->output_mesh, - &(tdy->output_mesh),NULL); CHKERRQ(ierr); - + "Enable output of mesh attributes","",options->output_mesh, + &(options->output_mesh),NULL); CHKERRQ(ierr); // Wrap up and indicate that options are set. ierr = PetscOptionsEnd(); CHKERRQ(ierr); tdy->setupflags |= TDyOptionsSet; - if (tdy->init_from_file && tdy->init_with_random_field) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Both -tdy_init_from_file and -tdy_init_with_random_field cannot be specified"); - } - ierr = TDyCreateGrid(tdy); CHKERRQ(ierr); PetscFunctionReturn(0); @@ -648,7 +654,7 @@ PetscErrorCode TDySetupDiscretizationScheme(TDy tdy) { PetscValidPointer(tdy,1); PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: ierr = TDyTPFInitialize(tdy); CHKERRQ(ierr); break; @@ -686,13 +692,13 @@ PetscErrorCode TDySetupNumericalMethods(TDy tdy) { TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); - if (tdy->regression_testing) { + if (tdy->options.regression_testing) { /* must come after Sections are set up in TDySetupDiscretizationScheme->XXXInitialize */ ierr = TDyRegressionInitialize(tdy); CHKERRQ(ierr); } - if (tdy->output_mesh) { - if (tdy->method != MPFA_O) { + if (tdy->options.output_mesh) { + if (tdy->options.method != MPFA_O) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, "-tdy_output_mesh only supported for MPFA-O method"); } @@ -706,35 +712,28 @@ PetscErrorCode TDySetupNumericalMethods(TDy tdy) { PetscErrorCode TDySetDiscretizationMethod(TDy tdy,TDyMethod method) { PetscFunctionBegin; - tdy->method = method; + tdy->options.method = method; PetscFunctionReturn(0); } PetscErrorCode TDySetMode(TDy tdy,TDyMode mode) { PetscValidPointer(tdy,1); PetscFunctionBegin; - tdy->mode = mode; + tdy->options.mode = mode; PetscFunctionReturn(0); } PetscErrorCode TDySetQuadratureType(TDy tdy,TDyQuadratureType qtype) { PetscValidPointer(tdy,1); PetscFunctionBegin; - tdy->qtype = qtype; + tdy->options.qtype = qtype; PetscFunctionReturn(0); } -PetscErrorCode TDySetWaterDensityType(TDy tdy,TDyWaterDensityType dentype) { +PetscErrorCode TDySetWaterDensityType(TDy tdy, TDyWaterDensityType dentype) { PetscValidPointer(tdy,1); PetscFunctionBegin; - switch (dentype) { - case WATER_DENSITY_CONSTANT: - tdy->rho_type = WATER_DENSITY_CONSTANT; - break; - case WATER_DENSITY_EXPONENTIAL: - tdy->rho_type = WATER_DENSITY_EXPONENTIAL; - break; - } + tdy->options.rho_type = dentype; PetscFunctionReturn(0); } @@ -742,7 +741,7 @@ PetscErrorCode TDySetMPFAOGmatrixMethod(TDy tdy,TDyMPFAOGmatrixMethod method) { PetscFunctionBegin; PetscValidPointer(tdy,1); - tdy->mpfao_gmatrix_method = method; + tdy->options.mpfao_gmatrix_method = method; PetscFunctionReturn(0); } @@ -751,7 +750,7 @@ PetscErrorCode TDySetMPFAOBoundaryConditionType(TDy tdy,TDyMPFAOBoundaryConditio PetscFunctionBegin; PetscValidPointer(tdy,1); - tdy->mpfao_bc_type = bctype; + tdy->options.mpfao_bc_type = bctype; PetscFunctionReturn(0); } @@ -774,14 +773,14 @@ PetscErrorCode TDySetIFunction(TS ts,TDy tdy) { ierr = PetscSectionGetNumFields(sec, &num_fields); ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"IFunction not implemented for TPF"); break; case MPFA_O: switch (dim) { case 3: - switch (tdy->mode) { + switch (tdy->options.mode) { case RICHARDS: ierr = TSSetIFunction(ts,NULL,TDyMPFAOIFunction_3DMesh,tdy); CHKERRQ(ierr); break; @@ -835,13 +834,13 @@ PetscErrorCode TDySetIJacobian(TS ts,TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)ts,&comm); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"IJacobian not implemented for TPF"); break; case MPFA_O: ierr = TDyCreateJacobian(tdy); CHKERRQ(ierr); - switch (tdy->mode) { + switch (tdy->options.mode) { case RICHARDS: ierr = TSSetIJacobian(ts,tdy->J,tdy->J,TDyMPFAOIJacobian_3DMesh,tdy); CHKERRQ(ierr); break; @@ -883,7 +882,7 @@ PetscErrorCode TDySetSNESFunction(SNES snes,TDy tdy) { ierr = PetscObjectGetComm((PetscObject)snes,&comm); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"SNESFunction not implemented for TPF"); break; @@ -927,7 +926,7 @@ PetscErrorCode TDySetSNESJacobian(SNES snes,TDy tdy) { ierr = PetscObjectGetComm((PetscObject)snes,&comm); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"SNESJacobian not implemented for TPF"); break; @@ -964,7 +963,7 @@ PetscErrorCode TDyComputeSystem(TDy tdy,Mat K,Vec F) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: ierr = TDyTPFComputeSystem(tdy,K,F); CHKERRQ(ierr); break; @@ -1002,7 +1001,7 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&P);CHKERRQ(ierr); ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&temp);CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { for (c=0;cK[i*dim2+j] = matprop->K0[i*dim2+j] * Kr; - ierr = ComputeWaterDensity(P[i], tdy->rho_type, &(tdy->rho[i]), &(tdy->drho_dP[i]), &(tdy->d2rho_dP2[i])); CHKERRQ(ierr); - ierr = ComputeWaterViscosity(P[i], tdy->mu_type, &(tdy->vis[i]), &(tdy->dvis_dP[i]), &(tdy->d2vis_dP2[i])); CHKERRQ(ierr); - if (tdy->mode == TH) { - for(j=0; jKappa[i*dim2+j] = matprop->Kappa0[i*dim2+j]; // update this based on Kersten number, etc. - ierr = ComputeWaterEnthalpy(temp[i], P[i], tdy->enthalpy_type, &(tdy->h[i]), &(tdy->dh_dP[i]), &(tdy->dh_dT[i])); CHKERRQ(ierr); + ierr = ComputeWaterDensity(P[i], tdy->options.rho_type, &(tdy->rho[i]), &(tdy->drho_dP[i]), &(tdy->d2rho_dP2[i])); CHKERRQ(ierr); + ierr = ComputeWaterViscosity(P[i], tdy->options.mu_type, &(tdy->vis[i]), &(tdy->dvis_dP[i]), &(tdy->d2vis_dP2[i])); CHKERRQ(ierr); + if (tdy->options.mode == TH) { + for(j=0; jKappa[i*dim2+j] = matprop->Kappa0[i*dim2+j]; // update this based on Kersten number, etc. + ierr = ComputeWaterEnthalpy(temp[i], P[i], tdy->options.enthalpy_type, &(tdy->h[i]), &(tdy->dh_dP[i]), &(tdy->dh_dT[i])); CHKERRQ(ierr); tdy->u[i] = tdy->h[i] - P[i]/tdy->rho[i]; } } - if ( (tdy->method == MPFA_O || tdy->method == MPFA_O_DAE || tdy->method == MPFA_O_TRANSIENTVAR) && dim == 3) { + if ( (tdy->options.method == MPFA_O || + tdy->options.method == MPFA_O_DAE || + tdy->options.method == MPFA_O_TRANSIENTVAR) && dim == 3) { PetscReal *p_vec_ptr, gz; TDyMesh *mesh = tdy->mesh; TDyCell *cells = &mesh->cells; @@ -1074,7 +1076,7 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { } ierr = VecRestoreArray(tdy->P_vec,&p_vec_ptr); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { PetscReal *t_vec_ptr; ierr = VecGetArray(tdy->Temp_P_vec, &t_vec_ptr); CHKERRQ(ierr); for (c=cStart; cdm),&comm); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: if(normp != NULL) { *normp = TDyTPFPressureNorm(tdy,U); } if(normv != NULL) { *normv = TDyTPFVelocityNorm(tdy,U); } @@ -1400,7 +1402,7 @@ PetscErrorCode TDyOutputRegression(TDy tdy, Vec U) { PetscErrorCode ierr; PetscFunctionBegin; - if (tdy->regression_testing){ + if (tdy->options.regression_testing){ ierr = TDyRegressionOutput(tdy,U); CHKERRQ(ierr); } PetscFunctionReturn(0); @@ -1452,7 +1454,7 @@ PetscErrorCode TDyPreSolveSNESSolver(TDy tdy) { ierr = PetscObjectGetComm((PetscObject)tdy->dm,&comm); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->method) { + switch (tdy->options.method) { case TPF: SETERRQ(comm,PETSC_ERR_SUP,"TDyPreSolveSNESSolver not implemented for TPF"); break; diff --git a/src/tdydriver.c b/src/tdydriver.c index 5f6e7214..64357828 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -21,7 +21,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver currently only supports 3D"); } - switch(tdy->method) { + switch(tdy->options.method) { case TPF: case MPFA_O: break; @@ -41,7 +41,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { ierr = TDySetPermeabilityFunction(tdy,TDyConstantPermeabilityFunction,PETSC_NULL); CHKERRQ(ierr); } - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { if (!TDyIsThermalConductivytSet(tdy)) { ierr = TDySetThermalConductivityFunction(tdy,TDyConstantThermalConductivityFunction,PETSC_NULL); CHKERRQ(ierr); @@ -64,7 +64,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { ierr = TDyCreateJacobian(tdy); CHKERRQ(ierr); // check for unsupported modes - switch (tdy->mode) { + switch (tdy->options.mode) { case RICHARDS: case TH: break; @@ -74,7 +74,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { // check for unsupported time integration methods switch(tdy->ti->time_integration_method) { case TDySNES: - switch (tdy->mode) { + switch (tdy->options.mode) { case RICHARDS: break; case TH: @@ -122,7 +122,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { // ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC); // CHKERRQ(ierr); // mode specific time integrator settings - switch (tdy->mode) { + switch (tdy->options.mode) { case RICHARDS: ierr = SNESLineSearchSetPostCheck(linesearch,TDyRichardsSNESPostCheck, &tdy); CHKERRQ(ierr); diff --git a/src/tdyregression.c b/src/tdyregression.c index d01a183b..2a02beee 100644 --- a/src/tdyregression.c +++ b/src/tdyregression.c @@ -24,7 +24,7 @@ PetscErrorCode TDyRegressionInitialize(TDy tdy) { ndof_per_cell = 1; - if (tdy->mode == TH) ndof_per_cell = 2; + if (tdy->options.mode == TH) ndof_per_cell = 2; regression = (TDyRegression *) malloc(sizeof(TDyRegression)); @@ -61,7 +61,7 @@ PetscErrorCode TDyRegressionInitialize(TDy tdy) { ierr = PetscMalloc(ndof_per_cell*size*regression->num_cells_per_process*sizeof(PetscInt),&(regression->cells_per_process_natural_ids)); CHKERRQ(ierr); } - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { increment = floor(vecsize_local/regression->num_cells_per_process/2); } else { @@ -170,14 +170,14 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { PetscReal *temp_p, *pres_p, *u_p; ierr = VecGetLocalSize(U,&vec_local_size); CHKERRQ(ierr); ierr = VecCreate(PetscObjectComm((PetscObject)dm),&U_pres); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { vec_local_size = vec_local_size/2; } ierr = VecSetSizes(U_pres,vec_local_size,PETSC_DECIDE); CHKERRQ(ierr); ierr = VecSetFromOptions(U_pres); CHKERRQ(ierr); ierr = VecGetArray(U_pres,&pres_p); CHKERRQ(ierr); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = VecCreate(PetscObjectComm((PetscObject)dm),&U_temp); CHKERRQ(ierr); ierr = VecSetSizes(U_temp,vec_local_size,PETSC_DECIDE); CHKERRQ(ierr); ierr = VecSetFromOptions(U_temp); CHKERRQ(ierr); @@ -187,7 +187,7 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { ierr = VecGetArray(U,&u_p); CHKERRQ(ierr); - if (tdy->mode == TH && num_fields == 2) { + if (tdy->options.mode == TH && num_fields == 2) { for (c=0;cmode == TH) { + if (tdy->options.mode == TH) { ierr = VecRestoreArray(U_temp,&temp_p); CHKERRQ(ierr); } @@ -211,7 +211,7 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { ierr = VecGetSize(U_pres,&pres_vec_size); CHKERRQ(ierr); mean_pres_val = mean_pres_val/pres_vec_size; - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { ierr = VecMax(U_temp,NULL,&max_temp_val); CHKERRQ(ierr); ierr = VecMin(U_temp,NULL,&min_temp_val); CHKERRQ(ierr); ierr = VecSum(U_temp,&mean_temp_val); CHKERRQ(ierr); @@ -237,7 +237,7 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { fprintf(fp," Min: %21.13e\n",min_pres_val); fprintf(fp," Mean: %21.13e\n",mean_pres_val); - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { for (i=0; icells_per_process_natural_ids[2*i]/2,vec_ptr[2*i]); } @@ -247,7 +247,7 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { } } - if (tdy->mode == TH) { + if (tdy->options.mode == TH) { fprintf(fp,"-- GENERIC: Temperature --\n"); fprintf(fp," Max: %21.13e\n",max_temp_val); fprintf(fp," Min: %21.13e\n",min_temp_val); @@ -263,7 +263,7 @@ PetscErrorCode TDyRegressionOutput(TDy tdy, Vec U) { } ierr = VecDestroy(&U_pres); CHKERRQ(ierr); - if (tdy->mode == TH) ierr = VecDestroy(&U_temp); CHKERRQ(ierr); + if (tdy->options.mode == TH) ierr = VecDestroy(&U_temp); CHKERRQ(ierr); TDY_STOP_FUNCTION_TIMER() PetscFunctionReturn(0); diff --git a/src/tdyrichards.c b/src/tdyrichards.c index eb52b4a8..89814a9e 100644 --- a/src/tdyrichards.c +++ b/src/tdyrichards.c @@ -8,15 +8,15 @@ PetscErrorCode TDyRichardsInitialize(TDy tdy) { PetscPrintf(PETSC_COMM_WORLD,"Running Richards mode.\n"); - if (tdy->init_with_random_field) { + if (tdy->options.init_with_random_field) { PetscRandom rand; ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand); CHKERRQ(ierr); ierr = PetscRandomSetInterval(rand,1.e4,1.e6); CHKERRQ(ierr); ierr = VecSetRandom(tdy->solution,rand); CHKERRQ(ierr); ierr = PetscRandomDestroy(&rand); CHKERRQ(ierr); - } else if (tdy->init_from_file) { + } else if (tdy->options.init_from_file) { PetscViewer viewer; - ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,tdy->init_file,FILE_MODE_READ,&viewer); CHKERRQ(ierr); + ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,tdy->options.init_file,FILE_MODE_READ,&viewer); CHKERRQ(ierr); ierr = VecLoad(tdy->solution,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); } else { diff --git a/src/tdyth.c b/src/tdyth.c index 18d5d0b6..21eac4fa 100644 --- a/src/tdyth.c +++ b/src/tdyth.c @@ -11,7 +11,7 @@ PetscErrorCode TDyTHInitialize(TDy tdy) { PetscPrintf(PETSC_COMM_WORLD,"Running TH mode.\n"); - if (tdy->init_with_random_field) { + if (tdy->options.init_with_random_field) { PetscRandom rand; ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand); CHKERRQ(ierr); ierr = VecGetLocalSize(tdy->solution,&local_size); CHKERRQ(ierr); diff --git a/src/tdyti.c b/src/tdyti.c index 40ab1185..0576ff86 100644 --- a/src/tdyti.c +++ b/src/tdyti.c @@ -117,7 +117,7 @@ PetscErrorCode TDyTimeIntegratorRunToTime(TDy tdy,PetscReal sync_time) { ierr = TDyTimeIntegratorSetTargetTime(ti,sync_time); CHKERRQ(ierr); ierr = TDySetDtimeForSNESSolver(tdy,ti->dt); CHKERRQ(ierr); if (!rank){ - switch (tdy->mode){ + switch (tdy->options.mode){ case RICHARDS: printf("===== RICHARDS MODE ==============================\n"); break; diff --git a/src/tdytimers.c b/src/tdytimers.c index cc62451a..b50c5c8f 100644 --- a/src/tdytimers.c +++ b/src/tdytimers.c @@ -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->method; - md->mode = tdy->mode; + md->method = tdy->options.method; + md->mode = tdy->options.mode; if (tdy->mesh != NULL) { md->num_cells = TDyMeshGetNumberOfLocalCells(tdy->mesh); } else { diff --git a/src/tdytpf.c b/src/tdytpf.c index ec54b96a..a2e83d2c 100644 --- a/src/tdytpf.c +++ b/src/tdytpf.c @@ -52,7 +52,7 @@ PetscErrorCode TDyTPFComputeSystem(TDy tdy,Mat K,Vec F) { DM dm = tdy->dm; PetscFunctionBegin; TDY_START_FUNCTION_TIMER() - if(!tdy->allow_unsuitable_mesh) { + if(!tdy->options.tpf_allow_all_meshes) { ierr = TDyTPFCheckMeshSuitability(tdy); CHKERRQ(ierr); } ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); From 39f8873cc743e7cb8f2b36dd31b97622ac5d8273 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 08:38:31 -0700 Subject: [PATCH 05/10] Touched up boundary condition args. --- include/private/tdyconditionsimpl.h | 10 ++++++++++ include/private/tdyoptions.h | 9 +++++++-- src/tdyconditions.c | 18 ++++++++++++++++++ src/tdycore.c | 25 ++++++++++++++++++++++--- 4 files changed, 57 insertions(+), 5 deletions(-) create mode 100644 include/private/tdyconditionsimpl.h diff --git a/include/private/tdyconditionsimpl.h b/include/private/tdyconditionsimpl.h new file mode 100644 index 00000000..44b3aca3 --- /dev/null +++ b/include/private/tdyconditionsimpl.h @@ -0,0 +1,10 @@ +#if !defined(TDYCONDITIONS_H) +#define TDYCONDITIONS_H + +#include + +PETSC_INTERN PetscErrorCode TDyConstantBoundaryPressureFn(TDy,PetscReal*,PetscReal*,void*); +PETSC_INTERN PetscErrorCode TDyConstantBoundaryVelocityFn(TDy,PetscReal*,PetscReal*,void*); + +#endif + diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index a77b324a..1d1fb308 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -19,19 +19,24 @@ typedef struct { PetscInt mpfao_bc_type; PetscBool tpf_allow_all_meshes; - // Material properties + // Constant material properties PetscReal porosity; PetscReal permeability; PetscReal soil_density; PetscReal soil_specific_heat; PetscReal thermal_conductivity; - // Characteristic curves + // Characteristic curve parameters PetscReal residual_saturation; PetscReal gardner_n; PetscReal vangenuchten_m; PetscReal vangenuchten_alpha; + // Constant boundary values + PetscReal boundary_pressure; + PetscReal boundary_temperature; + PetscReal boundary_velocity; // (normal component) + // I/O settings PetscBool init_with_random_field; PetscBool init_from_file; diff --git a/src/tdyconditions.c b/src/tdyconditions.c index 95c43009..951259f3 100644 --- a/src/tdyconditions.c +++ b/src/tdyconditions.c @@ -145,3 +145,21 @@ PetscErrorCode TDySetEnergySourceSinkValuesLocal(TDy tdy, PetscInt ni, const Pet PetscFunctionReturn(0); } +PetscErrorCode TDyConstantBoundaryPressureFn(TDy tdy, PetscReal *x, PetscReal *p, void *ctx) { + PetscFunctionBegin; + TDyOptions *options = &tdy->options; + + *p = options->boundary_pressure; + + PetscFunctionReturn(0); +} + +PetscErrorCode TDyConstantBoundaryVelocityFn(TDy tdy, PetscReal *x, PetscReal *v, void *ctx) { + PetscFunctionBegin; + TDyOptions *options = &tdy->options; + + *v = options->boundary_velocity; + + PetscFunctionReturn(0); +} + diff --git a/src/tdycore.c b/src/tdycore.c index 196c9284..8699f5de 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -209,6 +210,10 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { options->vangenuchten_m=0.8; options->vangenuchten_alpha=1.e-4; + options->boundary_pressure = 0.0; + options->boundary_temperature = 273.0; + options->boundary_velocity = 0.0; + options->init_with_random_field = PETSC_FALSE; options->init_from_file = PETSC_FALSE; @@ -572,19 +577,33 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { (PetscEnum)options->rho_type, (PetscEnum *)&options->rho_type, NULL); CHKERRQ(ierr); - // Named boundary condition functions + // Boundary conditions. char func_name[PETSC_MAX_PATH_LEN]; - ierr = PetscOptionsGetString(NULL, NULL, "-tdy_pressure_bc", func_name, + ierr = PetscOptionsGetString(NULL, NULL, "-tdy_pressure_bc_func", func_name, sizeof(func_name), &flag); CHKERRQ(ierr); if (flag) { // TODO: When/where are we supposed to get the context for this function? ierr = TDySelectBoundaryPressureFn(tdy, func_name, NULL); + } else { + ierr = PetscOptionsReal("-tdy_pressure_bc_value", "Constant boundary pressure", NULL, + options->boundary_pressure, &options->boundary_pressure, + NULL); CHKERRQ(ierr); + // Even if not given, we can set the boundary pressure to its default. + ierr = TDySetBoundaryPressureFn(tdy, TDyConstantBoundaryPressureFn, + PETSC_NULL); CHKERRQ(ierr); } - ierr = PetscOptionsGetString(NULL, NULL, "-tdy_velocity_bc", func_name, + ierr = PetscOptionsGetString(NULL, NULL, "-tdy_velocity_bc_func", func_name, sizeof(func_name), &flag); CHKERRQ(ierr); if (flag) { // TODO: When/where are we supposed to get the context for this function? ierr = TDySelectBoundaryVelocityFn(tdy, func_name, NULL); + } else { + ierr = PetscOptionsReal("-tdy_velocity_bc_value", "Constant normal boundary velocity", + NULL, options->boundary_pressure, &options->boundary_pressure, + NULL); CHKERRQ(ierr); + // Even if not given, we can set the boundary velocity to its default. + ierr = TDySetBoundaryVelocityFn(tdy, TDyConstantBoundaryVelocityFn, + PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); From 8a6171d9643e8ca1812d555ca1399894a1d8977e Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 09:19:45 -0700 Subject: [PATCH 06/10] Adding mesh generation/reading options. Not finished yet. --- include/private/tdyoptions.h | 9 +++++++- src/tdycore.c | 45 +++++++++++++++++++++++++++++++++--- 2 files changed, 50 insertions(+), 4 deletions(-) diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index 1d1fb308..48465502 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -37,10 +37,17 @@ typedef struct { PetscReal boundary_temperature; PetscReal boundary_velocity; // (normal component) - // I/O settings + // Initial conditions PetscBool init_with_random_field; PetscBool init_from_file; char init_file[PETSC_MAX_PATH_LEN]; + + // Mesh-related options + PetscBool generate_mesh; + PetscBool read_mesh; + char mesh_file[PETSC_MAX_PATH_LEN]; + + // I/O settings PetscBool output_mesh; PetscBool regression_testing; diff --git a/src/tdycore.c b/src/tdycore.c index 8699f5de..e0db1600 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -243,6 +243,7 @@ PetscErrorCode TDyCreate(TDy *_tdy) { tdy->Pref = 101325; tdy->Tref = 25; tdy->gravity[0] = 0; tdy->gravity[1] = 0; tdy->gravity[2] = 0; + tdy->dm = NULL; /* initialize method information to null */ tdy->vmap = NULL; tdy->emap = NULL; tdy->Alocal = NULL; tdy->Flocal = NULL; @@ -255,6 +256,15 @@ PetscErrorCode TDyCreate(TDy *_tdy) { PetscErrorCode TDySetDM(TDy tdy, DM dm) { PetscFunctionBegin; + 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 + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + "Can't call TDySetDM: A DM has already been read from a given file."); + } + } if (!dm) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A DM must be created prior to TDySetDM()"); } @@ -351,8 +361,14 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { if (!tdy->dm) { DM dm; - ierr = TDyCreateDM(&dm); CHKERRQ(ierr); - ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); + if (tdy->options.generate_mesh) { + // TODO: Let's jettison our own options in favor of PETSc's DM options. + ierr = TDyCreateDM(&dm); CHKERRQ(ierr); + ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); + } else { // if (tdy->options.read_mesh) + ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, tdy->options.mesh_file, + PETSC_TRUE, &dm); CHKERRQ(ierr); + } tdy->dm = dm; } @@ -605,7 +621,6 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = TDySetBoundaryVelocityFn(tdy, TDyConstantBoundaryVelocityFn, PETSC_NULL); CHKERRQ(ierr); } - ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options @@ -650,6 +665,30 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { "Only one of -tdy_init_from_file and -tdy_init_with_random_field can be specified"); } + // Mesh-related options + 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(PETSC_COMM_WORLD,PETSC_ERR_USER, + "Only one of -tdy_generate_mesh and -tdy_read_mesh can be specified"); + } + if ((tdy->dm != NULL) && (options->generate_mesh || options->read_mesh)) { + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + "TDySetDM was called before TDySetFromOptions: can't generate or " + "read a mesh"); + } + if ((tdy->dm == NULL) && !(options->generate_mesh || options->read_mesh)) { + SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, + "No mesh is available for TDycore: please use TDySetDM, " + "-tdy_generate_mesh, or -tdy_read_mesh to specify a mesh"); + } + + // Other options ierr = PetscOptionsBool("-tdy_regression_test", "Enable output of a regression file","", options->regression_testing, From 1de8f72ceec920c8a4adf30123d1f65a0dd2f124 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 09:53:23 -0700 Subject: [PATCH 07/10] Fixed a failing test. --- src/tdycore.c | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/tdycore.c b/src/tdycore.c index e0db1600..e6f9f7d0 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -603,10 +603,12 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { } else { ierr = PetscOptionsReal("-tdy_pressure_bc_value", "Constant boundary pressure", NULL, options->boundary_pressure, &options->boundary_pressure, - NULL); CHKERRQ(ierr); - // Even if not given, we can set the boundary pressure to its default. - ierr = TDySetBoundaryPressureFn(tdy, TDyConstantBoundaryPressureFn, - PETSC_NULL); CHKERRQ(ierr); + &flag); CHKERRQ(ierr); + if (flag) { + ierr = TDySetBoundaryPressureFn(tdy, TDyConstantBoundaryPressureFn, + PETSC_NULL); CHKERRQ(ierr); + } else { // TODO: what goes here?? + } } ierr = PetscOptionsGetString(NULL, NULL, "-tdy_velocity_bc_func", func_name, sizeof(func_name), &flag); CHKERRQ(ierr); @@ -616,10 +618,12 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { } else { ierr = PetscOptionsReal("-tdy_velocity_bc_value", "Constant normal boundary velocity", NULL, options->boundary_pressure, &options->boundary_pressure, - NULL); CHKERRQ(ierr); - // Even if not given, we can set the boundary velocity to its default. - ierr = TDySetBoundaryVelocityFn(tdy, TDyConstantBoundaryVelocityFn, - PETSC_NULL); CHKERRQ(ierr); + &flag); CHKERRQ(ierr); + if (flag) { + ierr = TDySetBoundaryVelocityFn(tdy, TDyConstantBoundaryVelocityFn, + PETSC_NULL); CHKERRQ(ierr); + } else { // TODO: what goes here?? + } } ierr = PetscOptionsEnd(); CHKERRQ(ierr); From 8f4433afc80c48941353aee9c7b1c6f0671a545c Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 14:29:42 -0700 Subject: [PATCH 08/10] We now use PETSc's mesh args for TDySetFromOptions. --- demo/mpfao/mpfao.cfg | 2 +- src/tdyconditions.c | 4 ++-- src/tdycore.c | 11 +++++++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/demo/mpfao/mpfao.cfg b/demo/mpfao/mpfao.cfg index 755fe2cd..2a796977 100644 --- a/demo/mpfao/mpfao.cfg +++ b/demo/mpfao/mpfao.cfg @@ -27,7 +27,7 @@ input_arguments=-problem 3 -tdy_regression_test -tdy_regression_test_num_cells_p input_arguments=-problem 4 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob4 [mpfao-64xy-3z-wedge] -input_arguments=-problem 4 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-64xy-3z-wedge -mesh_filename ../../share/meshes/64xy_3z_wedge.exo +input_arguments=-problem 4 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-64xy-3z-wedge -tdy_read_mesh ../../share/meshes/64xy_3z_wedge.exo [mpfao-3d-prob] input_arguments=-N 4 -dim 3 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-3d-prob diff --git a/src/tdyconditions.c b/src/tdyconditions.c index 951259f3..66a4539c 100644 --- a/src/tdyconditions.c +++ b/src/tdyconditions.c @@ -39,12 +39,12 @@ PetscErrorCode TDyGetFunction(const char* name, Function* f) { *f = kh_val(funcs_, iter); } else { ierr = -1; - SETERRQ(MPI_COMM_WORLD, ierr, "Function not found!"); + SETERRQ(PETSC_COMM_WORLD, ierr, "Function not found!"); return ierr; } } else { ierr = -1; - SETERRQ(MPI_COMM_WORLD, ierr, "No functions have been registered!"); + SETERRQ(PETSC_COMM_WORLD, ierr, "No functions have been registered!"); return ierr; } PetscFunctionReturn(0); diff --git a/src/tdycore.c b/src/tdycore.c index e6f9f7d0..24eb5af6 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -355,18 +355,21 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { PetscScalar *coords; PetscErrorCode ierr; PetscFunctionBegin; + MPI_Comm comm = PETSC_COMM_WORLD; if ((tdy->setupflags & TDyOptionsSet) == 0) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Options must be set prior to TDyCreateGrid()"); + SETERRQ(comm,PETSC_ERR_USER,"Options must be set prior to TDyCreateGrid()"); } if (!tdy->dm) { DM dm; if (tdy->options.generate_mesh) { - // TODO: Let's jettison our own options in favor of PETSc's DM options. - ierr = TDyCreateDM(&dm); CHKERRQ(ierr); + // 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(PETSC_COMM_WORLD, tdy->options.mesh_file, + ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, PETSC_TRUE, &dm); CHKERRQ(ierr); } tdy->dm = dm; From 25c73ca91a198d7e84a38874726a31ef9544a8f0 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 15:09:37 -0700 Subject: [PATCH 09/10] Eliminated TDyCreateDM. The only demo using this function was mpfao, which has been modified not to need it. Note the PETSc -dm_plex_* functions used as demo arguments. --- demo/mpfao/mpfao.c | 15 ++-- demo/mpfao/mpfao.cfg | 16 ++--- include/private/tdydmimpl.h | 1 - src/interface/ftn/tdydmf.c | 12 ---- src/tdydm.c | 134 ------------------------------------ 5 files changed, 17 insertions(+), 161 deletions(-) diff --git a/demo/mpfao/mpfao.c b/demo/mpfao/mpfao.c index 8c1570c2..57a80f0d 100644 --- a/demo/mpfao/mpfao.c +++ b/demo/mpfao/mpfao.c @@ -208,15 +208,18 @@ int main(int argc, char **argv) { ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); /* Create and distribute the mesh */ + MPI_Comm comm = PETSC_COMM_WORLD; DM dm; - ierr = TDyCreateDM(&dm); CHKERRQ(ierr); + ierr = DMCreate(comm, &dm); CHKERRQ(ierr); + ierr = DMSetType(dm, DMPLEX); CHKERRQ(ierr); + ierr = DMSetFromOptions(dm); CHKERRQ(ierr); if (perturb) {ierr = PerturbInteriorVertices(dm,perturbation); CHKERRQ(ierr);} ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); ierr = TDySetDM(tdy,dm); CHKERRQ(ierr); ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt ("-problem","Problem number","",problem,&problem,NULL); CHKERRQ(ierr); @@ -286,18 +289,18 @@ int main(int argc, char **argv) { #if defined(DEBUG) PetscViewer viewer; - ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"K.mat",&viewer); CHKERRQ(ierr); + ierr = PetscViewerASCIIOpen(comm,"K.mat",&viewer); CHKERRQ(ierr); ierr = MatView(K,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); - ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"trans.mat",&viewer); CHKERRQ(ierr); + ierr = PetscViewerASCIIOpen(comm,"trans.mat",&viewer); CHKERRQ(ierr); ierr = MatView(tdy->Trans_mat,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); #endif // Solve system KSP ksp; - ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr); + ierr = KSPCreate(comm,&ksp); CHKERRQ(ierr); ierr = KSPSetOperators(ksp,K,K); CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); ierr = KSPSetUp(ksp); CHKERRQ(ierr); @@ -305,7 +308,7 @@ int main(int argc, char **argv) { PetscReal normp, normv; ierr = TDyComputeErrorNorms(tdy,U,&normp,&normv); - ierr = PetscPrintf(PETSC_COMM_WORLD,"%e %e\n",normp,normv);CHKERRQ(ierr); + ierr = PetscPrintf(comm,"%e %e\n",normp,normv);CHKERRQ(ierr); ierr = TDyOutputRegression(tdy,U); diff --git a/demo/mpfao/mpfao.cfg b/demo/mpfao/mpfao.cfg index 2a796977..17b9737b 100644 --- a/demo/mpfao/mpfao.cfg +++ b/demo/mpfao/mpfao.cfg @@ -15,27 +15,27 @@ standard_parallel= pressure = 1.0e-12 relative [mpfao-prob1] -input_arguments=-problem 1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob1 -tdy_output_mesh +input_arguments=-problem 1 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob1 -tdy_output_mesh [mpfao-prob2] -input_arguments=-problem 2 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob2 +input_arguments=-problem 2 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob2 [mpfao-prob3] -input_arguments=-problem 3 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3 +input_arguments=-problem 3 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3 [mpfao-prob4] -input_arguments=-problem 4 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob4 +input_arguments=-problem 4 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob4 [mpfao-64xy-3z-wedge] -input_arguments=-problem 4 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-64xy-3z-wedge -tdy_read_mesh ../../share/meshes/64xy_3z_wedge.exo +input_arguments=-problem 4 -dm_plex_ѕimplex 0 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-64xy-3z-wedge -dm_plex_filename ../../share/meshes/64xy_3z_wedge.exo [mpfao-3d-prob] -input_arguments=-N 4 -dim 3 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-3d-prob +input_arguments=-dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 4,4,4 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-3d-prob [mpfao-prob3-np3] np=3 -input_arguments=-problem 3 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3-np3 +input_arguments=-problem 3 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3-np3 [mpfao-prob3-perturb-np3] np=3 -input_arguments=-problem 3 -perturb -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3-perturb-np3 +input_arguments=-problem 3 -dm_plex_simplex 0 -dm_plex_box_faces 8,8 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -perturb -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename mpfao-prob3-perturb-np3 diff --git a/include/private/tdydmimpl.h b/include/private/tdydmimpl.h index f7cb77ed..c39862b2 100644 --- a/include/private/tdydmimpl.h +++ b/include/private/tdydmimpl.h @@ -3,7 +3,6 @@ #include -PETSC_EXTERN PetscErrorCode TDyCreateDM(DM*); PETSC_EXTERN PetscErrorCode TDyDistributeDM(DM*); #endif diff --git a/src/interface/ftn/tdydmf.c b/src/interface/ftn/tdydmf.c index 6ceef371..9745bb2b 100644 --- a/src/interface/ftn/tdydmf.c +++ b/src/interface/ftn/tdydmf.c @@ -9,23 +9,11 @@ #include "tdycore.h" #ifdef PETSC_HAVE_FORTRAN_CAPS -#define tdycreatedm_ TDYCREATEDM #define tdydistributedm_ TDYDISTRIBUTEDM #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) && !defined(FORTRANDOUBLEUNDERSCORE) -#define tdycreatedm_ tdycreatedm #define tdydistributedm_ tdydistributedm #endif -#if defined(__cplusplus) -extern "C" { -#endif -PETSC_EXTERN void tdycreatedm_(DM *dm,int *__ierr){ -*__ierr = TDyCreateDM(dm); -} -#if defined(__cplusplus) -} -#endif - #if defined(__cplusplus) extern "C" { #endif diff --git a/src/tdydm.c b/src/tdydm.c index afeb8c4e..5c347240 100644 --- a/src/tdydm.c +++ b/src/tdydm.c @@ -1,140 +1,6 @@ #include #include -PetscErrorCode TDyCreateDM(DM *dm) { - PetscErrorCode ierr; - PetscInt dim = 2; - PetscReal lower_bound_x = 0., lower_bound_y = lower_bound_x, - lower_bound_z = lower_bound_x; - PetscReal upper_bound_x = 1., upper_bound_y = upper_bound_x, - upper_bound_z = upper_bound_x; - PetscInt Nx = -999, Ny = -999, Nz = -999; - PetscBool found; - char mesh_filename[PETSC_MAX_PATH_LEN]; - char string[512]; - char word[32]; - - mesh_filename[0] = '\0'; - - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"DM Options",""); - CHKERRQ(ierr); - ierr = PetscOptionsInt("-dim","Problem dimension","",dim,&dim,NULL); - CHKERRQ(ierr); - ierr = PetscOptionsInt("-N","Number of elements in 1D","",Nx,&Nx,&found); - CHKERRQ(ierr); - if (found) {Ny = Nx; Nz = Nx;} - ierr = PetscOptionsInt("-Nx","Number of elements in X","",Nx,&Nx,NULL); - CHKERRQ(ierr); - ierr = PetscOptionsInt("-Ny","Number of elements in Y","",Ny,&Ny,NULL); - CHKERRQ(ierr); - ierr = PetscOptionsInt("-Nz","Number of elements in Z","",Nz,&Nz,NULL); - CHKERRQ(ierr); - ierr = PetscOptionsReal("-lower_bound","Lower bound","",lower_bound_x, - &lower_bound_x,&found); CHKERRQ(ierr); - if (found) {lower_bound_y = lower_bound_x; lower_bound_z = lower_bound_x;} - ierr = PetscOptionsReal("-lower_bound_x","Lower bound in X","",lower_bound_x, - &lower_bound_x,NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-lower_bound_y","Lower bound in Y","",lower_bound_y, - &lower_bound_y,NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-lower_bound_z","Lower bound in Z","",lower_bound_z, - &lower_bound_z,NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-upper_bound","Upper bound","",upper_bound_x, - &upper_bound_x,&found); CHKERRQ(ierr); - if (found) {upper_bound_y = upper_bound_x; upper_bound_z = upper_bound_x;} - ierr = PetscOptionsReal("-upper_bound_x","Upper bound in X","",upper_bound_x, - &upper_bound_x,NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-upper_bound_y","Upper bound in Y","",upper_bound_y, - &upper_bound_y,NULL); CHKERRQ(ierr); - ierr = PetscOptionsReal("-upper_bound_z","Upper bound in Z","",upper_bound_z, - &upper_bound_z,NULL); CHKERRQ(ierr); - ierr = PetscOptionsString("-mesh_filename", "The mesh file", "", - mesh_filename, mesh_filename, PETSC_MAX_PATH_LEN, - NULL); CHKERRQ(ierr); - ierr = PetscOptionsEnd(); CHKERRQ(ierr); - - size_t len; - ierr = PetscStrlen(mesh_filename, &len); CHKERRQ(ierr); - if (!len){ - if (dim == 1) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "ERROR: Only two or three dimensions currently supported."); - int i = 0; - if (Nx > 0) i++; - if (Ny > 0) i++; - if (Nz > 0) i++; - if (i > 1) { - strcpy(string,"ERROR: Number of grid cells must be defined in only one dimension for a 1D problem:"); - if (Nx > 0) sprintf(word," %d",Nx); else sprintf(word," ?"); - strcat(string,word); - if (Ny > 0) sprintf(word," %d",Ny); else sprintf(word," ?"); - strcat(string,word); - if (Nz > 0) sprintf(word," %d\n",Nz); else sprintf(word," ?"); - strcat(string,word); - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,string); - } - if (Nx > 0) {Ny = 1; Nz = 1;} - else if (Ny > 0) {Nx = 1; Nz = 1;} - else if (Nz > 0) {Nx = 1; Ny = 1;} - else {Nx = 8; Ny = 1; Nz = 1;} - } else if (dim == 2) { - if (!((Nx > 0 && Ny > 0) || - (Nx > 0 && Nz > 0) || - (Ny > 0 && Nz > 0) || - (Nx > 0 && Ny < 0 && Nz < 0) || - (Nx < 0 && Ny < 0 && Nz < 0))) { - strcpy(string,"ERROR: Number of grid cells must be defined in one (-N #) or two dimensions for a 2D problem:"); - if (Nx > 0) sprintf(word," %d",Nx); else sprintf(word," ?"); - strcat(string,word); - if (Ny > 0) sprintf(word," %d",Ny); else sprintf(word," ?"); - strcat(string,word); - if (Nz > 0) sprintf(word," %d\n",Nz); else sprintf(word," ?"); - strcat(string,word); - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,string); - } - if (Nx > 0 && Ny > 0) Nz = 1; - else if (Nx > 0 && Nz > 0) Ny = 1; - else if (Ny > 0 && Nz > 0) Nx = 1; - else {Nx = 8; Ny = Nx; Nz = 1;} - } - else if (dim == 3) { - if ((Nx < 0 && (Ny > 0 || Nz > 0)) || - (Nx > 0 && ((Ny > 0 && Nz < 0) || (Ny < 0 && Nz > 0)))) { - strcpy(string,"ERROR: Number of grid cells must be defined in one (-N #) or three dimensions for a 3D problem:"); - if (Nx > 0) sprintf(word," %d",Nx); else sprintf(word," ?"); - strcat(string,word); - if (Ny > 0) sprintf(word," %d",Ny); else sprintf(word," ?"); - strcat(string,word); - if (Nz > 0) sprintf(word," %d\n",Nz); else sprintf(word," ?"); - strcat(string,word); - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,string); - } - if (Nx < 0) Nx = 8; - if (Ny < 0) Ny = Nx; - if (Nz < 0) Nz = Nx; - } - const PetscInt faces[3] = {Nx,Ny,Nz}; - const PetscReal lower[3] = {lower_bound_x,lower_bound_y,lower_bound_z}; - const PetscReal upper[3] = {upper_bound_x,upper_bound_y,upper_bound_z}; - - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, - lower, upper, NULL, PETSC_TRUE, dm); - CHKERRQ(ierr); - PetscInt rank; - ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)*dm),&rank); - if (!rank) { - int num_cell_in_set; - num_cell_in_set = Nx * Ny * Nz; - for (int c=0; c < num_cell_in_set; ++c) - ierr = DMSetLabelValue(*dm, "Cell Sets", c, 1); CHKERRQ(ierr); - } - } else { - ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh_filename, PETSC_TRUE, - dm); CHKERRQ(ierr); - } - - PetscFunctionReturn(0); -} - PetscErrorCode TDyDistributeDM(DM *dm) { DM dmDist; PetscErrorCode ierr; From 8c5c5336d28d0a74118c2e532d757307ba0d6e49 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 10 Aug 2021 15:51:51 -0700 Subject: [PATCH 10/10] Updated richards and th regression tests. --- demo/richards/richards.cfg | 8 ++++---- demo/richards/richards_driver.c | 6 +++--- demo/th/th.cfg | 8 ++++---- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/demo/richards/richards.cfg b/demo/richards/richards.cfg index 10fd49ee..adcef86a 100644 --- a/demo/richards/richards.cfg +++ b/demo/richards/richards.cfg @@ -11,17 +11,17 @@ standard_parallel= pressure = 1.0e-12 relative [richards-driver-snes-prob1] -input_arguments=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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 [richards-driver-ts-prob1] -input_arguments=-dim 3 -Nx 3 -Ny 3 -Nz 3 -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=-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 [richards-driver-snes-prob1-np4] np=4 -input_arguments=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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 [richards-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-dim 3 -Nx 3 -Ny 3 -Nz 3 -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=-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 diff --git a/demo/richards/richards_driver.c b/demo/richards/richards_driver.c index 9458da7f..1db1a184 100644 --- a/demo/richards/richards_driver.c +++ b/demo/richards/richards_driver.c @@ -9,7 +9,7 @@ int main(int argc, char **argv) { PetscBool print_intermediate = PETSC_FALSE; PetscMPIInt rank, size; TDy tdy = PETSC_NULL; - TDyIOFormat format = HDF5Format; + TDyIOFormat format = HDF5Format; ierr = TDyInit(argc, argv); CHKERRQ(ierr); ierr = TDyCreate(&tdy); CHKERRQ(ierr); @@ -19,7 +19,7 @@ int main(int argc, char **argv) { 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 = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code", "Code passed on successful completion","", @@ -44,7 +44,7 @@ int main(int argc, char **argv) { 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); diff --git a/demo/th/th.cfg b/demo/th/th.cfg index 46a9a165..de037b6a 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=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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 [th-driver-ts-prob1] -input_arguments=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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 #[th-driver-snes-prob1-np4] #np=4 -#input_arguments=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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 [th-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-dim 3 -Nx 2 -Ny 2 -Nz 2 -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=-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