diff --git a/dumux/io/plotmateriallaw.hh b/dumux/io/plotmateriallaw.hh index 7271fe152f3db3d022293acd032a40a831b8e090..ef87efa51631c7c0988f0f2a498df4351b0125c4 100644 --- a/dumux/io/plotmateriallaw.hh +++ b/dumux/io/plotmateriallaw.hh @@ -24,6 +24,8 @@ #ifndef DUMUX_PLOT_FLUID_MATRIX_LAW_HH #define DUMUX_PLOT_FLUID_MATRIX_LAW_HH +#include <dune/common/deprecated.hh> + #include <dumux/common/basicproperties.hh> #include <dumux/io/gnuplotinterface.hh> @@ -73,19 +75,24 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sw(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> sw; + std::vector<Scalar> pc; Scalar satInterval = upperSat - lowerSat; Scalar pcMin = 0.0; Scalar pcMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar swTemp, pcTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - pc[i] = MaterialLaw::pc(params, sw[i]); - pcMin = std::min(pcMin, pc[i]); - pcMax = std::max(pcMax, pc[i]); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + pcTemp = MaterialLaw::pc(params, swTemp); + if (checkValues_(swTemp, pcTemp)) + { + sw.push_back(swTemp); + pc.push_back(pcTemp); + pcMin = std::min(pcMin, pcTemp); + pcMax = std::max(pcMax, pcTemp); + } } gnuplotpcsw_.setXRange(lowerSat, upperSat); @@ -109,25 +116,31 @@ public: Scalar upperpc = 5000.0, std::string plotName = "") { - std::vector<Scalar> sat(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> sw; + std::vector<Scalar> pc; Scalar pcInterval = upperpc - lowerpc; Scalar swMin = 1e100; Scalar swMax = -1e100; + Scalar pcTemp, swTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - pc[i] = lowerpc + pcInterval * Scalar(i) / Scalar(numIntervals_); - sat[i] = MaterialLaw::sw(params, pc[i]); - swMin = std::min(swMin, sat[i]); - swMax = std::max(swMax, sat[i]); + pcTemp = lowerpc + pcInterval * Scalar(i) / Scalar(numIntervals_); + swTemp = MaterialLaw::sw(params, pcTemp); + if (checkValues_(pcTemp, swTemp)) + { + pc.push_back(pcTemp); + sw.push_back(swTemp); + swMin = std::min(swMin, swTemp); + swMax = std::max(swMax, swTemp); + } } gnuplotswpc_.setXRange(lowerpc, upperpc); gnuplotswpc_.setYRange(swMin, swMax); gnuplotswpc_.setXlabel("capillary pressure [Pa]"); gnuplotswpc_.setYlabel("wetting phase saturation [-]"); - gnuplotswpc_.addDataSetToPlot(pc, sat, plotName + "_Sw-pc"); + gnuplotswpc_.addDataSetToPlot(pc, sw, plotName + "_Sw-pc"); gnuplotswpc_.plot("sw-pc"); } @@ -144,26 +157,31 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sat(numIntervals_ + 1); - std::vector<Scalar> dpcdsw(numIntervals_ + 1); + std::vector<Scalar> sw; + std::vector<Scalar> dpcdsw; Scalar satInterval = upperSat - lowerSat; Scalar dpcdswMin = 1e100; Scalar dpcdswMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar swTemp, dpcdswTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sat[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - dpcdsw[i] = MaterialLaw::dpc_dsw(params, sat[i]); - dpcdswMin = std::min(dpcdswMin, dpcdsw[i]); - dpcdswMax = std::max(dpcdswMax, dpcdsw[i]); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + dpcdswTemp = MaterialLaw::dpc_dsw(params, swTemp); + if (checkValues_(swTemp, dpcdsw)) + { + sw.push_back(swTemp); + dpcdsw.push_back(dpcdswTemp); + dpcdswMin = std::min(dpcdswMin, dpcdswTemp); + dpcdswMax = std::max(dpcdswMax, dpcdswTemp); + } } gnuplotdpcdsw_.setXRange(lowerSat, upperSat); gnuplotdpcdsw_.setYRange(dpcdswMin, dpcdswMax); gnuplotdpcdsw_.setXlabel("wetting phase saturation [-]"); gnuplotdpcdsw_.setYlabel("gradient of the pc-Sw curve [Pa]"); - gnuplotdpcdsw_.addDataSetToPlot(sat, dpcdsw, plotName + "_dpcdSw-Sw"); + gnuplotdpcdsw_.addDataSetToPlot(sw, dpcdsw, plotName + "_dpcdSw-Sw"); gnuplotdpcdsw_.plot("dpcdsw"); } @@ -180,18 +198,24 @@ public: Scalar upperpc = 5000.0, std::string plotName = "") { - std::vector<Scalar> dswdpc(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> pc; + std::vector<Scalar> dswdpc; Scalar pcInterval = upperpc - lowerpc; Scalar dswdpcMin = 1e100; Scalar dswdpcMax = -1e100; + Scalar dswdpcTemp, pcTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - pc[i] = lowerpc + pcInterval * Scalar(i) / Scalar(numIntervals_); - dswdpc[i] = MaterialLaw::dsw_dpc(params, pc[i]); - dswdpcMin = std::min(dswdpcMin, dswdpc[i]); - dswdpcMax = std::max(dswdpcMax, dswdpc[i]); + pcTemp = lowerpc + pcInterval * Scalar(i) / Scalar(numIntervals_); + dswdpcTemp = MaterialLaw::dsw_dpc(params, pcTemp); + if (checkValues_(pcTemp, dswdpcTemp)) + { + pc.push_back(pcTemp); + dswdpc.push_back(dswdpcTemp); + dswdpcMin = std::min(dswdpcMin, dswdpcTemp); + dswdpcMax = std::max(dswdpcMax, dswdpcTemp); + } } gnuplotdswdpc_.setXRange(lowerpc, upperpc); @@ -215,21 +239,27 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sw(numIntervals_ + 1); - std::vector<Scalar> krw(numIntervals_ + 1); - std::vector<Scalar> krn(numIntervals_ + 1); + std::vector<Scalar> sw; + std::vector<Scalar> krw; + std::vector<Scalar> krn; Scalar satInterval = upperSat - lowerSat; Scalar krMin = 1e100; Scalar krMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName + "_kr"); + Scalar swTemp, krwTemp, krnTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - krw[i] = MaterialLaw::krw(params, sw[i]); - krn[i] = MaterialLaw::krn(params, sw[i]); - krMin = std::min(krMin, std::min(krw[i], krn[i])); - krMax = std::max(krMax, std::max(krw[i], krn[i])); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + krwTemp = MaterialLaw::krw(params, swTemp); + krnTemp = MaterialLaw::krn(params, swTemp); + if (checkValues_(swTemp, krwTemp) && checkValues_(swTemp, krnTemp)) + { + sw.push_back(swTemp); + krw.push_back(krwTemp); + krn.push_back(krnTemp); + krMin = std::min({krMin, krwTemp, krnTemp}); + krMax = std::max({krMax, krwTemp, krnTemp}); + } } gnuplotkr_.setXRange(lowerSat, upperSat); @@ -254,21 +284,27 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sw(numIntervals_+1); - std::vector<Scalar> dkrw_dsw(numIntervals_+1); - std::vector<Scalar> dkrn_dsw(numIntervals_+1); + std::vector<Scalar> sw; + std::vector<Scalar> dkrw_dsw; + std::vector<Scalar> dkrn_dsw; Scalar satInterval = upperSat - lowerSat; Scalar dkrdswMin = 1e100; Scalar dkrdswMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName + "_dkr_dsw"); + Scalar swTemp, dkrwdswTemp, dkrndswTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - dkrw_dsw[i] = MaterialLaw::dkrw_dsw(params, sw[i]); - dkrn_dsw[i] = MaterialLaw::dkrn_dsw(params, sw[i]); - dkrdswMin = std::min(dkrdswMin, std::min(dkrw_dsw[i], dkrn_dsw[i])); - dkrdswMax = std::max(dkrdswMax, std::max(dkrw_dsw[i], dkrn_dsw[i])); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + dkrwdswTemp = MaterialLaw::dkrw_dsw(params, swTemp); + dkrndswTemp = MaterialLaw::dkrn_dsw(params, swTemp); + if (checkValues_(swTemp, dkrwdswTemp) && checkValues_(swTemp, dkrndswTemp)) + { + sw.push_back(swTemp); + dkrw_dsw.push_back(dkrwdswTemp); + dkrn_dsw.push_back(dkrndswTemp); + dkrdswMin = std::min({dkrdswMin, dkrwdswTemp, dkrndswTemp}); + dkrdswMax = std::max({dkrdswMax, dkrwdswTemp, dkrndswTemp}); + } } gnuplotkrdsw_.setXRange(lowerSat, upperSat); @@ -289,6 +325,7 @@ public: * \param upperSat Maximum x-value * \param plotName Name of the plotted curve */ + DUNE_DEPRECATED_MSG("checkEffectiveSaturation is deprecated.") void checkEffectiveSaturation(const MaterialLawParams ¶ms, Scalar lowerSat, Scalar upperSat, @@ -301,6 +338,18 @@ public: } private: + /*! + * \brief Check the values for occurrences of nan and inf + * + * \param value1 A data point value + * \param value2 An other data point value + */ + bool checkValues_(Scalar value1, Scalar value2) + { + return !std::isnan(value1) && !std::isinf(value1) + && !std::isnan(value2) && !std::isinf(value2); + } + int numIntervals_; GnuplotInterface<Scalar> gnuplotpcsw_; GnuplotInterface<Scalar> gnuplotswpc_; diff --git a/dumux/io/plotmateriallaw3p.hh b/dumux/io/plotmateriallaw3p.hh index 1cb520bb757777a1a870c6392fbcaa64e0084f1d..7da7ab693d89fc0b3906a4294ab8a1c1839cd644 100644 --- a/dumux/io/plotmateriallaw3p.hh +++ b/dumux/io/plotmateriallaw3p.hh @@ -24,6 +24,8 @@ #ifndef DUMUX_PLOT_FLUID_MATRIX_LAW_HH #define DUMUX_PLOT_FLUID_MATRIX_LAW_HH +#include <dune/common/deprecated.hh> + #include <dumux/common/basicproperties.hh> #include <dumux/io/gnuplotinterface.hh> @@ -91,19 +93,24 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sw(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> sw; + std::vector<Scalar> pc; Scalar satInterval = upperSat - lowerSat; Scalar pcMin = 0.0; Scalar pcMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar swTemp, pcTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - pc[i] = MaterialLaw::pcgw(params, sw[i]); - pcMin = std::min(pcMin, pc[i]); - pcMax = std::max(pcMax, pc[i]); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + pcTemp = MaterialLaw::pcgw(params, swTemp); + if (checkValues_(swTemp, pcTemp)) + { + sw.push_back(swTemp); + pc.push_back(pcTemp); + pcMin = std::min(pcMin, pcTemp); + pcMax = std::max(pcMax, pcTemp); + } } gnuplotpc_.setXRange(lowerSat, upperSat); @@ -127,19 +134,24 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> sw(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> sw; + std::vector<Scalar> pc; Scalar satInterval = upperSat - lowerSat; Scalar pcMin = 0.0; Scalar pcMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar swTemp, pcTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - pc[i] = MaterialLaw::pcnw(params, sw[i]); - pcMin = std::min(pcMin, pc[i]); - pcMax = std::max(pcMax, pc[i]); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + pcTemp = MaterialLaw::pcnw(params, swTemp); + if (checkValues_(swTemp, pcTemp)) + { + sw.push_back(swTemp); + pc.push_back(pcTemp); + pcMin = std::min(pcMin, pcTemp); + pcMax = std::max(pcMax, pcTemp); + } } gnuplotpc_.setXRange(lowerSat, upperSat); @@ -163,19 +175,24 @@ public: Scalar upperSat = 1.0, std::string plotName = "") { - std::vector<Scalar> st(numIntervals_+1); - std::vector<Scalar> pc(numIntervals_+1); + std::vector<Scalar> st; + std::vector<Scalar> pc; Scalar satInterval = upperSat - lowerSat; Scalar pcMin = 0.0; Scalar pcMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar stTemp, pcTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - st[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - pc[i] = MaterialLaw::pcgn(params, st[i]); - pcMin = std::min(pcMin, pc[i]); - pcMax = std::max(pcMax, pc[i]); + stTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + pcTemp = MaterialLaw::pcgn(params, stTemp); + if (checkValues_(stTemp, pcTemp)) + { + st.push_back(stTemp); + pc.push_back(pcTemp); + pcMin = std::min(pcMin, pcTemp); + pcMax = std::max(pcMax, pcTemp); + } } gnuplotpc_.setXRange(lowerSat, upperSat); @@ -207,16 +224,25 @@ public: Scalar satInterval = upperSat - lowerSat; Scalar krMin = 1e100; Scalar krMax = -1e100; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName + "_kr"); + Scalar swTemp, krwTemp, krnTemp, krgTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sw[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - krw[i] = MaterialLaw::krw(params, sw[i], 0.0); - krn[i] = MaterialLaw::krn(params, sw[i], 1-sw[i]); - krg[i] = MaterialLaw::krg(params, sw[i], 0.0); - krMin = std::min(krMin, std::min({krw[i], krn[i], krg[i]})); - krMax = std::max(krMax, std::max({krw[i], krn[i], krg[i]})); + swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + krwTemp = MaterialLaw::krw(params, swTemp, 0.0); + krnTemp = MaterialLaw::krn(params, swTemp, 1.0 - swTemp); + krgTemp = MaterialLaw::krg(params, swTemp, 0.0); + if (checkValues_(swTemp, krwTemp) + && checkValues_(swTemp, krnTemp) + && checkValues_(swTemp, krgTemp)) + { + sw.push_back(swTemp); + krw.push_back(krwTemp); + krn.push_back(krnTemp); + krg.push_back(krgTemp); + krMin = std::min({krMin, krwTemp, krnTemp, krgTemp}); + krMax = std::max({krMax, krwTemp, krnTemp, krgTemp}); + } } gnuplotkr_.setXRange(lowerSat, upperSat); @@ -247,14 +273,19 @@ public: Scalar satInterval = upperSat - lowerSat; Scalar alphaMin = -2; Scalar alphaMax = 2; - checkEffectiveSaturation(params, lowerSat, upperSat, plotName); + Scalar snTemp, alphaTemp = 0.0; for (int i = 0; i <= numIntervals_; i++) { - sn[i] = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); - alpha[i] = MaterialLaw::pcAlpha(params, sn[i]); - alphaMin = std::min(alphaMin, alpha[i]); - alphaMax = std::max(alphaMax, alpha[i]); + snTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_); + alphaTemp = MaterialLaw::pcAlpha(params, snTemp); + if (checkValues_(snTemp, alphaTemp)) + { + sn.push_back(snTemp); + alpha.push_back(alphaTemp); + alphaMin = std::min(alphaMin, alphaTemp); + alphaMax = std::max(alphaMax, alphaTemp); + } } gnuplotpcAlpha_.setXRange(lowerSat, upperSat); @@ -274,6 +305,7 @@ public: * \param upperSat Maximum x-value * \param plotName Name of the plotted curve */ + DUNE_DEPRECATED_MSG("checkEffectiveSaturation is deprecated.") void checkEffectiveSaturation(const MaterialLawParams ¶ms, Scalar lowerSat, Scalar upperSat, @@ -286,6 +318,18 @@ public: } private: + /*! + * \brief Check the values for occurrences of nan and inf + * + * \param value1 A data point value + * \param value2 An other data point value + */ + bool checkValues_(Scalar value1, Scalar value2) + { + return !std::isnan(value1) && !std::isinf(value1) + && !std::isnan(value2) && !std::isinf(value2); + } + int numIntervals_; GnuplotInterface<Scalar> gnuplotpc_; GnuplotInterface<Scalar> gnuplotpcAlpha_; diff --git a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh index 7a73d40545438b5a058dddb3f8fcf506ef051079..4dc9f7878080c4f4d9208f7e69c811b93b8d6176 100644 --- a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh +++ b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh @@ -143,6 +143,8 @@ public: FluidSystem::init(282.15, 284.15, 3, 8e4, 3e5, 200); name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name); + + this->spatialParams().plotMaterialLaw(); } /*! diff --git a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh index ea8b9fcc01d14286f9cc02088065c42da44d187b..c7bc3e18a863d3b6d0f2337fa4e992a66d89b201 100644 --- a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh +++ b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh @@ -30,6 +30,7 @@ #include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh> #include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3pparams.hh> #include <dumux/material/fluidmatrixinteractions/3p/efftoabslaw.hh> +#include <dumux/io/plotmateriallaw3p.hh> namespace Dumux { @@ -117,11 +118,26 @@ public: // parameters for adsorption materialParams_.setKdNAPL(0.); materialParams_.setRhoBulk(1500.); + + plotFluidMatrixInteractions_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Output, + PlotFluidMatrixInteractions); } ~InfiltrationThreePSpatialParams() {} + /*! + * \brief This is called from the problem and creates a gnuplot output + * of e.g the pc-Sw curve + */ + void plotMaterialLaw() + { + PlotMaterialLaw<TypeTag> plotMaterialLaw(plotFluidMatrixInteractions_); + + plotMaterialLaw.plotpc(materialParams_); + plotMaterialLaw.plotkr(materialParams_); + } + /*! * \brief Intrinsic permability * @@ -188,6 +204,8 @@ private: Scalar porosity_; MaterialLawParams materialParams_; + + bool plotFluidMatrixInteractions_; }; } diff --git a/test/porousmediumflow/3p/implicit/test_box3p.input b/test/porousmediumflow/3p/implicit/test_box3p.input index edd8a86f233adbb0cbaa2229763fe06de8a7f24e..6445dff4ae034abcd0ce6f15d05dc78e0c9e4349 100644 --- a/test/porousmediumflow/3p/implicit/test_box3p.input +++ b/test/porousmediumflow/3p/implicit/test_box3p.input @@ -15,3 +15,6 @@ porosity = 0.40 vanGenuchtenAlpha = 0.0005 vanGenuchtenN = 4.0 +[Output] +PlotFluidMatrixInteractions = false + diff --git a/test/porousmediumflow/3p/implicit/test_cc3p.input b/test/porousmediumflow/3p/implicit/test_cc3p.input index 08652caa0f1966db45b69feded08c298bc2ec9e7..764e1a3f0c4ab0036e27c25bc2f9177e9a33673e 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3p.input +++ b/test/porousmediumflow/3p/implicit/test_cc3p.input @@ -15,3 +15,6 @@ porosity = 0.40 vanGenuchtenAlpha = 0.0005 vanGenuchtenN = 4.0 +[Output] +PlotFluidMatrixInteractions = false +