Commit 3f347747 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/richards-analytical-jacobian' into 'master'

Feature/richards analytical jacobian

See merge request !1697
parents a43a81c4 356e90e2
......@@ -264,24 +264,43 @@ public:
*/
static Scalar dswe_dpc(const Params &params, Scalar pc)
{
// TODO: This is most probably still buggy!!
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
// Genuchten's law
Scalar sw;
if (pc < 0)
sw = 1.5; // make sure we regularize below
else
sw = VanGenuchten::sw(params, pc);
const Scalar swThLow = params.pcLowSw();
const Scalar swThHigh = params.pcHighSw();
if (pc <= 0)
{
// for swThHigh = 1.0 the slope gets infinity
// swThHigh > 1.0 are not sensible threshold values
// setting swThHigh = 1.0 is a way to disable regularization
// return the maximum representable number with the given precision
if (swThHigh > 1.0 - std::numeric_limits<Scalar>::epsilon())
return std::numeric_limits<Scalar>::max();
Scalar yTh = VanGenuchten::pc(params, swThHigh);
Scalar m1 = (0.0 - yTh)/(1.0 - swThHigh)*2;
return 1.0/m1;
}
const auto sw = VanGenuchten::sw(params, pc);
// derivative of the regularization
if (sw < params.pcLowSw()) {
// same as in dpc_dswe() but inverted
if (sw <= swThLow)
return 1/mLow_(params);
}
if (sw > params.pcHighSw()) {
// same as in dpc_dswe() but inverted
return 1/mHigh_(params);
if (sw > swThHigh)
{
const Scalar yTh = VanGenuchten::pc(params, swThHigh);
const Scalar m1 = (0.0 - yTh)/(1.0 - swThHigh)*2;
// invert spline between threshold swe and 1.0
const Scalar mTh = VanGenuchten::dpc_dswe(params, swThHigh);
const Spline<Scalar> sp(swThHigh, 1.0, // x0, x1
yTh, 0, // m0, m1
mTh, m1); // m0, m1
const auto swReg = sp.intersectInterval(swThHigh, 1.0,
0, 0, 0, pc);
// derivative of the inverse of the function is one over derivative of the function
return 1.0/sp.evalDerivative(swReg);
}
return VanGenuchten::dswe_dpc(params, pc);
......@@ -337,15 +356,14 @@ public:
// phase from the parameters
const Scalar swThHigh = params.krwHighSw();
// derivative of the regualarization
if (swe < 0) {
// the slope is zero
// derivative of the regularization
// the slope is zero below sw=0.0 and above sw=1.0
if (swe < 0)
return 0.0;
}
else if (swe > 1 - std::numeric_limits<Scalar>::epsilon()) {
// the slope is zero
else if (swe > 1 - std::numeric_limits<Scalar>::epsilon())
return 0.0;
}
// for high sw we need the slope of the interpolation spline
else if (swe > swThHigh) {
using Spline = Dumux::Spline<Scalar>;
Spline sp(swThHigh, 1.0, // x1, x2
......
......@@ -27,6 +27,7 @@
#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
#include <dumux/common/typetraits/typetraits.hh>
namespace Dumux {
......@@ -68,6 +69,15 @@ class RichardsLocalResidual : public GetPropType<TypeTag, Properties::BaseLocalR
static constexpr bool enableWaterDiffusionInAir
= getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>();
//! An element solution that does not compile if the [] operator is used
struct InvalidElemSol
{
template<class Index>
double operator[] (const Index i) const
{ static_assert(AlwaysFalse<Index>::value, "Solution-dependent material parameters not supported with analytical differentiation"); return 0.0; }
};
public:
using ParentType::ParentType;
......@@ -150,6 +160,324 @@ public:
return flux;
}
/*!
* \brief Adds the storage derivative
*
* \param partialDerivatives The partial derivatives
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curVolVars The current volume variables
* \param scv The sub control volume
*/
template<class PartialDerivativeMatrix>
void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curVolVars,
const SubControlVolume& scv) const
{
static_assert(!enableWaterDiffusionInAir,
"richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
static_assert(!FluidSystem::isCompressible(0),
"richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
const auto poreVolume = scv.volume()*curVolVars.porosity();
static const auto rho = curVolVars.density(0);
// material law for pc-sw
const auto& matParams = problem.spatialParams().materialLawParams(element, scv, InvalidElemSol{});
// partial derivative of storage term w.r.t. p_w
// d(Sw*rho*phi*V/dt)/dpw = rho*phi*V/dt*dsw/dpw = rho*phi*V/dt*dsw/dpc*dpc/dpw = -rho*phi*V/dt*dsw/dpc
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*MaterialLaw::dsw_dpc(matParams, curVolVars.capillaryPressure());
}
/*!
* \brief Adds source derivatives for wetting and non-wetting phase.
*
* \param partialDerivatives The partial derivatives
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curVolVars The current volume variables
* \param scv The sub control volume
*
* \todo Maybe forward to problem for the user to implement the source derivatives?
*/
template<class PartialDerivativeMatrix>
void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curVolVars,
const SubControlVolume& scv) const
{ /* TODO maybe forward to problem for the user to implement the source derivatives?*/ }
/*!
* \brief Adds flux derivatives for wetting and non-wetting phase for cell-centered FVM using TPFA
*
* Compute derivatives for the wetting and the non-wetting phase flux with respect to \f$p_w\f$
* and \f$S_n\f$.
*
* \param derivativeMatrices The partial derivatives
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curElemVolVars The current element volume variables
* \param elemFluxVarsCache The element flux variables cache
* \param scvf The sub control volume face
*/
template<class PartialDerivativeMatrices, class T = TypeTag>
std::enable_if_t<GetPropType<T, Properties::FVGridGeometry>::discMethod == DiscretizationMethod::cctpfa, void>
addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
static_assert(!enableWaterDiffusionInAir,
"richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
static_assert(!FluidSystem::isCompressible(0),
"richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
static_assert(!FluidSystem::viscosityIsConstant(0),
"richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
// get references to the two participating vol vars & parameters
const auto insideScvIdx = scvf.insideScvIdx();
const auto outsideScvIdx = scvf.outsideScvIdx();
const auto outsideElement = fvGeometry.fvGridGeometry().element(outsideScvIdx);
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& insideVolVars = curElemVolVars[insideScvIdx];
const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement, outsideScv, InvalidElemSol{});
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
static const auto rho = insideVolVars.density(0);
static const auto mu = insideVolVars.viscosity(0);
static const auto rho_mu = rho/mu;
// upwind term
// evaluate the current wetting phase Darcy flux and resulting upwind weights
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
const auto outsideWeight = 1.0 - insideWeight;
const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
// material law derivatives
const auto insideSw = insideVolVars.saturation(0);
const auto outsideSw = outsideVolVars.saturation(0);
const auto insidePc = insideVolVars.capillaryPressure();
const auto outsidePc = outsideVolVars.capillaryPressure();
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
const auto dkrw_dsw_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
const auto dsw_dpw_outside = -MaterialLaw::dsw_dpc(outsideMaterialParams, outsidePc);
// the transmissibility
const auto tij = elemFluxVarsCache[scvf].advectionTij();
// get references to the two participating derivative matrices
auto& dI_dI = derivativeMatrices[insideScvIdx];
auto& dI_dJ = derivativeMatrices[outsideScvIdx];
// partial derivative of the wetting phase flux w.r.t. pw
dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
}
/*!
* \brief Adds flux derivatives for box method
*
* \param A The Jacobian Matrix
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curElemVolVars The current element volume variables
* \param elemFluxVarsCache The element flux variables cache
* \param scvf The sub control volume face
*/
template<class JacobianMatrix, class T = TypeTag>
std::enable_if_t<GetPropType<T, Properties::FVGridGeometry>::discMethod == DiscretizationMethod::box, void>
addFluxDerivatives(JacobianMatrix& A,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
static_assert(!enableWaterDiffusionInAir,
"richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
static_assert(!FluidSystem::isCompressible(0),
"richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
static_assert(!FluidSystem::viscosityIsConstant(0),
"richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
// get references to the two participating vol vars & parameters
const auto insideScvIdx = scvf.insideScvIdx();
const auto outsideScvIdx = scvf.outsideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& insideVolVars = curElemVolVars[insideScvIdx];
const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(element, outsideScv, InvalidElemSol{});
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
static const auto rho = insideVolVars.density(0);
static const auto mu = insideVolVars.viscosity(0);
static const auto rho_mu = rho/mu;
// upwind term
// evaluate the current wetting phase Darcy flux and resulting upwind weights
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
const auto outsideWeight = 1.0 - insideWeight;
const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
// material law derivatives
const auto insideSw = insideVolVars.saturation(0);
const auto outsideSw = outsideVolVars.saturation(0);
const auto insidePc = insideVolVars.capillaryPressure();
const auto outsidePc = outsideVolVars.capillaryPressure();
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
const auto dkrw_dsw_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
const auto dsw_dpw_outside = -MaterialLaw::dsw_dpc(outsideMaterialParams, outsidePc);
// so far it was the same as for tpfa
// the transmissibilities (flux derivatives with respect to all pw-dofs on the element)
const auto ti = AdvectionType::calculateTransmissibilities(
problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]
);
// get the rows of the jacobian matrix for the inside/outside scv
auto& dI_dJ_inside = A[insideScv.dofIndex()];
auto& dI_dJ_outside = A[outsideScv.dofIndex()];
// add the partial derivatives w.r.t all scvs in the element
for (const auto& scvJ : scvs(fvGeometry))
{
// the transmissibily associated with the scvJ
const auto& tj = ti[scvJ.indexInElement()];
const auto globalJ = scvJ.dofIndex();
// partial derivative of the wetting phase flux w.r.t. p_w
dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm;
dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm;
}
const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
// additional constribution of the upwind term only for inside and outside dof
A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside;
A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside;
A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside;
A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside;
}
/*!
* \brief Adds cell-centered Dirichlet flux derivatives for wetting and non-wetting phase
*
* Compute derivatives for the wetting and the non-wetting phase flux with respect to \f$p_w\f$
* and \f$S_n\f$.
*
* \param derivativeMatrices The matrices containing the derivatives
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curElemVolVars The current element volume variables
* \param elemFluxVarsCache The element flux variables cache
* \param scvf The sub control volume face
*/
template<class PartialDerivativeMatrices>
void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
static_assert(!enableWaterDiffusionInAir,
"richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
static_assert(!FluidSystem::isCompressible(0),
"richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
static_assert(!FluidSystem::viscosityIsConstant(0),
"richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
// get references to the two participating vol vars & parameters
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = curElemVolVars[insideScvIdx];
const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
static const auto rho = insideVolVars.density(0);
static const auto mu = insideVolVars.viscosity(0);
static const auto rho_mu = rho/mu;
// upwind term
// evaluate the current wetting phase Darcy flux and resulting upwind weights
using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>;
static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
const auto outsideWeight = 1.0 - insideWeight;
const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
// material law derivatives
const auto insideSw = insideVolVars.saturation(0);
const auto insidePc = insideVolVars.capillaryPressure();
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
// the transmissibility
const auto tij = elemFluxVarsCache[scvf].advectionTij();
// partial derivative of the wetting phase flux w.r.t. pw
derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
}
/*!
* \brief Adds Robin flux derivatives for wetting and non-wetting phase
*
* \param derivativeMatrices The matrices containing the derivatives
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume element geometry
* \param curElemVolVars The current element volume variables
* \param elemFluxVarsCache The element flux variables cache
* \param scvf The sub control volume face
*/
template<class PartialDerivativeMatrices>
void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{ /* TODO maybe forward to problem for the user to implement the Robin derivatives?*/ }
private:
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
......
......@@ -2,10 +2,10 @@ add_input_file_links(FILES params.input)
# isothermal tests
add_executable(test_richards_lens_tpfa EXCLUDE_FROM_ALL main.cc)
target_compile_definitions(test_richards_lens_tpfa PUBLIC TYPETAG=RichardsLensCC)
target_compile_definitions(test_richards_lens_tpfa PUBLIC TYPETAG=RichardsLensCC PUBLIC DIFFMETHOD=DiffMethod::numeric)
add_executable(test_richards_lens_box EXCLUDE_FROM_ALL main.cc)
target_compile_definitions(test_richards_lens_box PUBLIC TYPETAG=RichardsLensBox)
target_compile_definitions(test_richards_lens_box PUBLIC TYPETAG=RichardsLensBox PUBLIC DIFFMETHOD=DiffMethod::numeric)
dune_add_test(TARGET test_richards_lens_box
LABELS porousmediumflow richards
......@@ -15,6 +15,17 @@ dune_add_test(TARGET test_richards_lens_box
${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_box-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_box params.input -Problem.Name test_richards_lens_box")
dune_add_test(NAME test_richards_lens_box_analyticdiff
SOURCES main.cc
LABELS porousmediumflow richards
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic
COMPILE_DEFINITIONS TYPETAG=RichardsLensBox
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_richards_lens_box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_box_analyticdiff-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_box_analyticdiff params.input -Problem.Name test_richards_lens_box_analyticdiff")
dune_add_test(TARGET test_richards_lens_tpfa
LABELS porousmediumflow richards
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
......@@ -23,6 +34,17 @@ dune_add_test(TARGET test_richards_lens_tpfa
${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_tpfa-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_tpfa params.input -Problem.Name test_richards_lens_tpfa")
dune_add_test(NAME test_richards_lens_tpfa_analyticdiff
SOURCES main.cc
LABELS porousmediumflow richards
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic
COMPILE_DEFINITIONS TYPETAG=RichardsLensCC
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_richards_lens_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_tpfa_analyticdiff-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_richards_lens_tpfa_analyticdiff params.input -Problem.Name test_richards_lens_tpfa_analyticdiff")
dune_add_test(NAME test_richards_lens_box_parallel_yasp
TARGET test_richards_lens_box
LABELS porousmediumflow richards parallel
......
......@@ -50,6 +50,10 @@
#include "problem.hh"
#ifndef DIFFMETHOD
#define DIFFMETHOD DiffMethod::numeric
#endif
////////////////////////
// the main function
////////////////////////
......@@ -136,7 +140,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, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
......@@ -176,6 +180,9 @@ int main(int argc, char** argv) try
timeLoop->finalize(leafGridView.comm());
// write Newton report
nonLinearSolver.report();
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment