diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
index 05f2a24c02a866646d0b49ba1e9b9c5bc1690427..d4b1559ecc2a74ce75272e3474a9892e536e7ced 100644
--- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
@@ -103,17 +103,12 @@ public:
static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1,
"2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!");
- // we know that these values are constant throughout the simulation
const auto poreVolume = scv.volume()*curVolVars.porosity();
- static const std::array<Scalar, numPhases> rhos = { curVolVars.density(0), curVolVars.density(1) };
+ const auto dt = this->timeLoop().timeStepSize();
- for (int pIdx = 0; pIdx < ModelTraits::numFluidPhases(); ++pIdx)
- {
- // partial derivative of wetting phase storage term w.r.t. p_w
- partialDerivatives[conti0EqIdx+pIdx][pressureIdx] += 0.0;
- // partial derivative of wetting phase storage term w.r.t. S_n
- partialDerivatives[conti0EqIdx+pIdx][saturationIdx] -= poreVolume*rhos[pIdx]/this->timeLoop().timeStepSize();
- }
+ // partial derivative of the phase storage terms w.r.t. S_n
+ partialDerivatives[conti0EqIdx+0][saturationIdx] -= poreVolume*curVolVars.density(0)/dt;
+ partialDerivatives[conti0EqIdx+1][saturationIdx] += poreVolume*curVolVars.density(1)/dt;
}
/*!
@@ -216,12 +211,12 @@ public:
// derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
const auto insideSw = insideVolVars.saturation(0);
const auto outsideSw = outsideVolVars.saturation(0);
- const auto dKrw_dSn_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
- const auto dKrw_dSn_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
- const auto dKrn_dSn_inside = MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
- const auto dKrn_dSn_outside = MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
- const auto dpc_dSn_inside = MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
- const auto dpc_dSn_outside = MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
+ const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
+ const auto dKrw_dSn_outside = -1.0*MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
+ const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
+ const auto dKrn_dSn_outside = -1.0*MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
+ const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+ const auto dpc_dSn_outside = -1.0*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
const auto tij = elemFluxVarsCache[scvf].advectionTij();
@@ -238,20 +233,20 @@ public:
dI_dJ[conti0EqIdx+0][pressureIdx] -= tij_up_w;
// partial derivative of the wetting phase flux w.r.t. S_n
- dI_dI[conti0EqIdx+0][saturationIdx] -= rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
- dI_dJ[conti0EqIdx+0][saturationIdx] -= rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
+ dI_dI[conti0EqIdx+0][saturationIdx] += rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
+ dI_dJ[conti0EqIdx+0][saturationIdx] += rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
// partial derivative of the non-wetting phase flux w.r.t. p_w
dI_dI[conti0EqIdx+1][pressureIdx] += tij_up_n;
dI_dJ[conti0EqIdx+1][pressureIdx] -= tij_up_n;
// partial derivative of the non-wetting phase flux w.r.t. S_n (relative permeability derivative contribution)
- dI_dI[conti0EqIdx+1][saturationIdx] -= rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
- dI_dJ[conti0EqIdx+1][saturationIdx] -= rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
+ dI_dI[conti0EqIdx+1][saturationIdx] += rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
+ dI_dJ[conti0EqIdx+1][saturationIdx] += rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
// partial derivative of the non-wetting phase flux w.r.t. S_n (capillary pressure derivative contribution)
- dI_dI[conti0EqIdx+1][saturationIdx] -= tij_up_n*dpc_dSn_inside;
- dI_dJ[conti0EqIdx+1][saturationIdx] += tij_up_n*dpc_dSn_outside;
+ dI_dI[conti0EqIdx+1][saturationIdx] += tij_up_n*dpc_dSn_inside;
+ dI_dJ[conti0EqIdx+1][saturationIdx] -= tij_up_n*dpc_dSn_outside;
}
/*!
@@ -369,39 +364,39 @@ public:
{
// partial derivative of the wetting phase flux w.r.t. S_n
const auto insideSw = insideVolVars.saturation(0);
- const auto dKrw_dSn_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
+ const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
- dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
- dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
+ dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
+ dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
// partial derivative of the non-wetting phase flux w.r.t. S_n (k_rn contribution)
- const auto dKrn_dSn_inside = MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
+ const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
- dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
- dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxN_dSnJ_krn;
+ dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn;
+ dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
// partial derivative of the non-wetting phase flux w.r.t. S_n (p_c contribution)
- const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
- dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
- dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
+ const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
+ dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
+ dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
}
else if (localJ == outsideScvIdx)
{
// see comments for (globalJ == insideScvIdx)
const auto outsideSw = outsideVolVars.saturation(0);
- const auto dKrw_dSn_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
+ const auto dKrw_dSn_outside = -1.0*MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
- dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
- dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
+ dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ;
+ dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ;
- const auto dKrn_dSn_outside = MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
+ const auto dKrn_dSn_outside = -1.0*MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
- dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
- dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn;
+ dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn;
+ dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn;
- const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
- dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
- dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
+ const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
+ dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc;
+ dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc;
}
else
{
diff --git a/test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt b/test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt
index 42f790e2e498448e0a0d08bb3903affb5e878dc2..88209079d5c9ce0a90409942d8b7068451cd1ae3 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt
+++ b/test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt
@@ -13,6 +13,16 @@ dumux_add_test(NAME test_2p_incompressible_tpfa
${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa params.input -Problem.Name test_2p_incompressible_tpfa")
+dumux_add_test(NAME test_2p_incompressible_tpfa_analytic
+ SOURCES main.cc
+ LABELS porousmediumflow 2p
+ COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleTpfa DIFFMETHOD=DiffMethod::analytic
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2p_incompressible_cc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa_analytic-00007.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa_analytic params.input -Problem.Name test_2p_incompressible_tpfa_analytic -Newton.EnablePartialReassembly false")
+
# using tpfa
dumux_add_test(NAME test_2p_incompressible_tpfa_restart
TARGET test_2p_incompressible_tpfa
@@ -37,6 +47,16 @@ dumux_add_test(NAME test_2p_incompressible_box
${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box params.input -Problem.Name test_2p_incompressible_box")
+dumux_add_test(NAME test_2p_incompressible_box_analytic
+ LABELS porousmediumflow 2p
+ SOURCES main.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleBox DIFFMETHOD=DiffMethod::analytic
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2p_incompressible_box-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box_analytic-00007.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box params.input -Problem.Name test_2p_incompressible_box_analytic -Newton.EnablePartialReassembly false")
+
# using box with interface solver
dumux_add_test(NAME test_2p_incompressible_box_ifsolver
LABELS porousmediumflow 2p
diff --git a/test/porousmediumflow/2p/implicit/incompressible/main.cc b/test/porousmediumflow/2p/implicit/incompressible/main.cc
index 676ac8e40e558f71b06f8ede10540557ec167ab2..847924989c5e5150b428d20033eb3f1184de70a4 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/main.cc
+++ b/test/porousmediumflow/2p/implicit/incompressible/main.cc
@@ -52,6 +52,10 @@
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/loadsolution.hh>
+#ifndef DIFFMETHOD
+#define DIFFMETHOD DiffMethod::numeric
+#endif
+
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -172,7 +176,7 @@ int main(int argc, char** argv) try
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
- using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ using Assembler = FVAssembler<TypeTag, DIFFMETHOD>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver