From 53bb3c89e66b06dd424fa72ce2330dbf77daa109 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 27 Jul 2018 17:16:52 +0200
Subject: [PATCH] [2pncni][test] avoid code duplication by reusing main and
problem file
---
.../2pnc/implicit/CMakeLists.txt | 4 +-
.../2pnc/implicit/fuelcellniproblem.hh | 341 ------------------
.../2pnc/implicit/fuelcellproblem.hh | 21 +-
.../2pnc/implicit/test_2pncni_fv.cc | 232 ------------
4 files changed, 21 insertions(+), 577 deletions(-)
delete mode 100644 test/porousmediumflow/2pnc/implicit/fuelcellniproblem.hh
delete mode 100644 test/porousmediumflow/2pnc/implicit/test_2pncni_fv.cc
diff --git a/test/porousmediumflow/2pnc/implicit/CMakeLists.txt b/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
index 74aba18232..2fd12c6134 100644
--- a/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2pnc/implicit/CMakeLists.txt
@@ -38,8 +38,8 @@ dune_add_test(NAME test_cc2pnc_fickslaw
--command "${CMAKE_CURRENT_BINARY_DIR}/test_cc2pnc_fickslaw test_2pnc_diffusion.input -Problem.Name test_cc2pnc_fickslaw")
dune_add_test(NAME test_2pncni_box
- SOURCES test_2pncni_fv.cc
- COMPILE_DEFINITIONS TYPETAG=FuelCellNIBoxTypeTag
+ SOURCES test_2pnc_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=FuelCellNIBoxTypeTag NONISOTHERMAL=1
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/fuelcell2pncboxni-reference.vtu
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcellniproblem.hh b/test/porousmediumflow/2pnc/implicit/fuelcellniproblem.hh
deleted file mode 100644
index fb2f38b604..0000000000
--- a/test/porousmediumflow/2pnc/implicit/fuelcellniproblem.hh
+++ /dev/null
@@ -1,341 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup TwoPNCTests
- * \brief Definition of a problem for water management in PEM fuel cells.
- */
-#ifndef DUMUX_FUELCELL_NI_PROBLEM_HH
-#define DUMUX_FUELCELL_NI_PROBLEM_HH
-
-#include <dune/grid/yaspgrid.hh>
-
-#include <dumux/discretization/elementsolution.hh>
-#include <dumux/discretization/cellcentered/tpfa/properties.hh>
-#include <dumux/discretization/box/properties.hh>
-#include <dumux/porousmediumflow/2pnc/model.hh>
-#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/material/fluidsystems/h2on2o2.hh>
-#include <dumux/material/chemistry/electrochemistry/electrochemistryni.hh>
-
-#include "fuelcellspatialparams.hh"
-
-namespace Dumux {
-
-/*!
- * \ingroup TwoPNCTests
- * \brief Definition of a problem for water management in PEM fuel cells.
- */
-template <class TypeTag>
-class FuelCellNIProblem;
-
-namespace Properties {
-NEW_TYPE_TAG(FuelCellNITypeTag, INHERITS_FROM(TwoPNCNI));
-NEW_TYPE_TAG(FuelCellNIBoxTypeTag, INHERITS_FROM(BoxModel, FuelCellNITypeTag));
-NEW_TYPE_TAG(FuelCellNICCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, FuelCellNITypeTag));
-
-// Set the grid type
-SET_TYPE_PROP(FuelCellNITypeTag, Grid, Dune::YaspGrid<2>);
-
-// Set the problem property
-SET_TYPE_PROP(FuelCellNITypeTag, Problem, FuelCellNIProblem<TypeTag>);
-
-// Set the spatial parameters
-SET_TYPE_PROP(FuelCellNITypeTag, SpatialParams, FuelCellSpatialParams<TypeTag>);
-
-// Set the primary variable combination for the 2pnc model
-SET_PROP(FuelCellNITypeTag, Formulation)
-{ static constexpr auto value = TwoPFormulation::p1s0; };
-
-// Set fluid configuration
-SET_PROP(FuelCellNITypeTag, FluidSystem)
-{
-private:
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
- using type = FluidSystems::H2ON2O2<Scalar>;
-};
-} // end namespace Properties
-
-/*!
- * \ingroup TwoPNCModel
- * \ingroup ImplicitTestProblems
- * \brief Problem or water management in PEM fuel cells.
- *
- * To run the simulation execute the following line in shell:
- * <tt>./test_box2pnc</tt>
- */
-template <class TypeTag>
-class FuelCellNIProblem : public PorousMediumFlowProblem<TypeTag>
-{
- using ParentType = PorousMediumFlowProblem<TypeTag>;
-
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Element = typename GridView::template Codim<0>::Entity;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- // Select the electrochemistry method
- using ElectroChemistry = typename Dumux::ElectroChemistryNI<Scalar, Indices, FluidSystem, FVGridGeometry, ElectroChemistryModel::Ochs>;
-
- static constexpr int dim = GridView::dimension;
- static constexpr int dimWorld = GridView::dimensionworld;
- static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
- using GlobalPosition = typename SubControlVolume::GlobalPosition;
-
- enum { dofCodim = isBox ? dim : 0 };
-public:
- /*!
- * \brief The constructor
- *
- * \param timeManager The time manager
- * \param gridView The grid view
- */
- FuelCellNIProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
- : ParentType(fvGridGeometry)
- {
- nTemperature_ = getParam<int>("Problem.NTemperature");
- nPressure_ = getParam<int>("Problem.NPressure");
- pressureLow_ = getParam<Scalar>("Problem.PressureLow");
- pressureHigh_ = getParam<Scalar>("Problem.PressureHigh");
- temperatureLow_ = getParam<Scalar>("Problem.TemperatureLow");
- temperatureHigh_ = getParam<Scalar>("Problem.TemperatureHigh");
- temperature_ = getParam<Scalar>("Problem.InitialTemperature");
-
- name_ = getParam<std::string>("Problem.Name");
-
- pO2Inlet_ = getParam<Scalar>("ElectroChemistry.pO2Inlet");
-
- FluidSystem::init(/*Tmin=*/temperatureLow_,
- /*Tmax=*/temperatureHigh_,
- /*nT=*/nTemperature_,
- /*pmin=*/pressureLow_,
- /*pmax=*/pressureHigh_,
- /*np=*/nPressure_);
- }
-
- /*!
- * \name Problem parameters
- */
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const std::string& name() const
- { return name_; }
-
- /*!
- * \brief Returns the temperature within the domain.
- *
- * This problem assumes a temperature of 10 degrees Celsius.
- */
- Scalar temperature() const
- { return temperature_; }
-
- //! \copydoc Dumux::FVProblem::source()
- NumEqVector source(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv) const
- {
- NumEqVector values(0.0);
- const auto& globalPos = scv.dofPosition();
-
- //reaction sources from electro chemistry
- if(inReactionLayer_(globalPos))
- {
- const auto& volVars = elemVolVars[scv];
- auto currentDensity = ElectroChemistry::calculateCurrentDensity(volVars);
- ElectroChemistry::reactionSource(values, currentDensity);
- }
-
- return values;
- }
-
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment
- *
- * \param globalPos The global position
- */
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
- {
- BoundaryTypes bcTypes;
- bcTypes.setAllNeumann();
-
- if (onUpperBoundary_(globalPos)){
- bcTypes.setAllDirichlet();
- }
- return bcTypes;
- }
-
- /*!
- * \brief Evaluates the boundary conditions for a Dirichlet
- * boundary segment
- *
- * \param globalPos The global position
- */
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
- {
- auto priVars = initial_(globalPos);
-
- if(onUpperBoundary_(globalPos))
- {
- Scalar pn = 1.0e5;
- priVars[Indices::pressureIdx] = pn;
- priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
- priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
- priVars[Indices::temperatureIdx] = 293.15;
- }
-
- return priVars;
- }
-
- /*!
- * \name Volume terms
- */
-
-
- /*!
- * \brief Evaluates the initial values for a control volume
- *
- * \param globalPos The global position
- */
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
- { return initial_(globalPos); }
-
-
- /*!
- * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
- */
- template<class VTKWriter>
- void addVtkFields(VTKWriter& vtk)
- {
- const auto& gridView = this->fvGridGeometry().gridView();
- currentDensity_.resize(gridView.size(dofCodim));
- reactionSourceH2O_.resize(gridView.size(dofCodim));
- reactionSourceO2_.resize(gridView.size(dofCodim));
- Kxx_.resize(gridView.size(dofCodim));
- Kyy_.resize(gridView.size(dofCodim));
-
- vtk.addField(currentDensity_, "currentDensity [A/cm^2]");
- vtk.addField(reactionSourceH2O_, "reactionSourceH2O [mol/(sm^2)]");
- vtk.addField(reactionSourceO2_, "reactionSourceO2 [mol/(sm^2)]");
- vtk.addField(Kxx_, "Kxx");
- vtk.addField(Kyy_, "Kyy");
- }
-
- void updateVtkFields(const SolutionVector& curSol)
- {
- for (const auto& element : elements(this->fvGridGeometry().gridView()))
- {
- auto elemSol = elementSolution(element, curSol, this->fvGridGeometry());
-
- auto fvGeometry = localView(this->fvGridGeometry());
- fvGeometry.bindElement(element);
-
- for (auto&& scv : scvs(fvGeometry))
- {
- VolumeVariables volVars;
- volVars.update(elemSol, *this, element, scv);
- const auto& globalPos = scv.dofPosition();
- const auto dofIdxGlobal = scv.dofIndex();
-
- if(inReactionLayer_(globalPos))
- {
- //reactionSource Output
- PrimaryVariables source;
- auto i = ElectroChemistry::calculateCurrentDensity(volVars);
- ElectroChemistry::reactionSource(source, i);
-
- reactionSourceH2O_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::H2OIdx];
- reactionSourceO2_[dofIdxGlobal] = source[Indices::conti0EqIdx + FluidSystem::O2Idx];
-
- //Current Output in A/cm^2
- currentDensity_[dofIdxGlobal] = i/10000;
- }
- else
- {
- reactionSourceH2O_[dofIdxGlobal] = 0.0;
- reactionSourceO2_[dofIdxGlobal] = 0.0;
- currentDensity_[dofIdxGlobal] = 0.0;
- }
- Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
- Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
- }
- }
- }
-
-private:
-
- PrimaryVariables initial_(const GlobalPosition &globalPos) const
- {
- PrimaryVariables priVars(0.0);
- priVars.setState(Indices::bothPhases);
-
- Scalar pn = 1.0e5;
- priVars[Indices::pressureIdx] = pn;
- priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
- priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
- priVars[Indices::temperatureIdx] = 293.15;
-
- return priVars;
- }
-
- bool onUpperBoundary_(const GlobalPosition &globalPos) const
- { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
-
- bool inReactionLayer_(const GlobalPosition& globalPos) const
- { return globalPos[1] < 0.1*(this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]) + eps_; }
-
- Scalar temperature_;
- static constexpr Scalar eps_ = 1e-6;
- int nTemperature_;
- int nPressure_;
- std::string name_ ;
- Scalar pressureLow_, pressureHigh_;
- Scalar temperatureLow_, temperatureHigh_;
- Scalar pO2Inlet_;
- std::vector<double> currentDensity_;
- std::vector<double> reactionSourceH2O_;
- std::vector<double> reactionSourceO2_;
- std::vector<double> Kxx_;
- std::vector<double> Kyy_;
-};
-
-} //end namespace Dumux
-
-#endif
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
index b82263d893..9d72928587 100644
--- a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
+++ b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
@@ -32,8 +32,11 @@
#include <dumux/porousmediumflow/2pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2o2.hh>
+#ifdef NONISOTHERMAL
+#include <dumux/material/chemistry/electrochemistry/electrochemistryni.hh>
+#else
#include <dumux/material/chemistry/electrochemistry/electrochemistry.hh>
-
+#endif
#include "fuelcellspatialparams.hh"
namespace Dumux {
@@ -46,9 +49,14 @@ template <class TypeTag>
class FuelCellProblem;
namespace Properties {
+#ifdef NONISOTHERMAL
+NEW_TYPE_TAG(FuelCellTypeTag, INHERITS_FROM(TwoPNCNI));
+NEW_TYPE_TAG(FuelCellNIBoxTypeTag, INHERITS_FROM(BoxModel, FuelCellTypeTag));
+#else
NEW_TYPE_TAG(FuelCellTypeTag, INHERITS_FROM(TwoPNC));
NEW_TYPE_TAG(FuelCellBoxTypeTag, INHERITS_FROM(BoxModel, FuelCellTypeTag));
NEW_TYPE_TAG(FuelCellCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, FuelCellTypeTag));
+#endif
// Set the grid type
SET_TYPE_PROP(FuelCellTypeTag, Grid, Dune::YaspGrid<2>);
@@ -101,8 +109,11 @@ class FuelCellProblem : public PorousMediumFlowProblem<TypeTag>
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
// Select the electrochemistry method
+#ifdef NONISOTHERMAL
+ using ElectroChemistry = typename Dumux::ElectroChemistryNI<Scalar, Indices, FluidSystem, FVGridGeometry, ElectroChemistryModel::Ochs>;
+#else
using ElectroChemistry = typename Dumux::ElectroChemistry<Scalar, Indices, FluidSystem, FVGridGeometry, ElectroChemistryModel::Ochs>;
-
+#endif
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
@@ -218,6 +229,9 @@ public:
priVars[Indices::pressureIdx] = pn;
priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
+#ifdef NONISOTHERMAL
+ priVars[Indices::temperatureIdx] = 293.15;
+#endif
}
return priVars;
@@ -309,6 +323,9 @@ private:
priVars[Indices::pressureIdx] = pn;
priVars[Indices::switchIdx] = 0.3;//Sw for bothPhases
priVars[Indices::switchIdx+1] = pO2Inlet_/4.315e9; //moleFraction xlO2 for bothPhases
+#ifdef NONISOTHERMAL
+ priVars[Indices::temperatureIdx] = 293.15;
+#endif
return priVars;
}
diff --git a/test/porousmediumflow/2pnc/implicit/test_2pncni_fv.cc b/test/porousmediumflow/2pnc/implicit/test_2pncni_fv.cc
deleted file mode 100644
index 58ceab7d60..0000000000
--- a/test/porousmediumflow/2pnc/implicit/test_2pncni_fv.cc
+++ /dev/null
@@ -1,232 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Test for the 2pnc cc model used for water management in PEM fuel cells.
- */
-#include <config.h>
-
-#include <ctime>
-#include <iostream>
-
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/timer.hh>
-#include <dune/grid/io/file/dgfparser/dgfexception.hh>
-#include <dune/grid/io/file/vtk.hh>
-#include <dune/istl/io.hh>
-
-#include "fuelcellproblem.hh"
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/common/valgrind.hh>
-#include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/defaultusagemessage.hh>
-
-#include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
-
-#include <dumux/assembly/fvassembler.hh>
-#include <dumux/assembly/diffmethod.hh>
-
-#include <dumux/discretization/methods.hh>
-
-#include <dumux/io/vtkoutputmodule.hh>
-#include <dumux/io/grid/gridmanager.hh>
-
-#include "fuelcellniproblem.hh"
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- * reading in parameters.
- *
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
- if (errorMsg.size() > 0) {
- std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
- "\t-ParameterFile Parameter file (Input file) \n";
-
- std::cout << errorMessageOut
- << "\n";
- }
-}
-
-int main(int argc, char** argv) try
-{
- using namespace Dumux;
-
- // define the type tag for this problem
- using TypeTag = TTAG(TYPETAG);
-
- // initialize MPI, finalize is done automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
-
- // print dumux start message
- if (mpiHelper.rank() == 0)
- DumuxMessage::print(/*firstCall=*/true);
-
- // parse command line arguments and input file
- Parameters::init(argc, argv, usage);
-
- // try to create a grid (from the given grid file or the input file)
- GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
- gridManager.init();
-
- ////////////////////////////////////////////////////////////
- // run instationary non-linear problem on this grid
- ////////////////////////////////////////////////////////////
-
- // we compute on the leaf grid view
- const auto& leafGridView = gridManager.grid().leafGridView();
-
- // create the finite volume grid geometry
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
- auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
- fvGridGeometry->update();
-
- // the problem (initial and boundary conditions)
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- auto problem = std::make_shared<Problem>(fvGridGeometry);
-
- // the solution vector
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- SolutionVector x(fvGridGeometry->numDofs());
- problem->applyInitialSolution(x);
- auto xOld = x;
-
- // the grid variables
- using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
- auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
- gridVariables->init(x, xOld);
-
- // get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
- auto dt = getParam<Scalar>("TimeLoop.DtInitial");
-
- // check if we are about to restart a previously interrupted simulation
- Scalar restartTime = 0;
- if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
- restartTime = getParam<Scalar>("TimeLoop.Restart");
-
- // initialize the vtk output module
- using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
- VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
- using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
- vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
- VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
- problem->addVtkFields(vtkWriter); //!< Add problem specific output fields
- vtkWriter.write(0.0);
-
- // instantiate time loop
- auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
-
- // the assembler with time loop for instationary problem
- using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
-
- // the linear solver
- using LinearSolver = AMGBackend<TypeTag>;
- auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
-
- // the non-linear solver
- using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
- typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch)>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
-
- // time loop
- timeLoop->start(); do
- {
- // set previous solution for storage evaluations
- assembler->setPreviousSolution(xOld);
-
- // solve the non-linear system with time step control
- nonLinearSolver.solve(x, *timeLoop);
-
- // make the new solution the old solution
- xOld = x;
- gridVariables->advanceTimeStep();
-
- // advance to the time loop to the next step
- timeLoop->advanceTimeStep();
-
- // update the output fields with the current solution before write
- problem->updateVtkFields(x);
-
- // write vtk output
- vtkWriter.write(timeLoop->time());
-
- // report statistics of this time step
- timeLoop->reportTimeStep();
-
- // set new dt as suggested by the newton solver
- timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
-
- } while (!timeLoop->finished());
-
- timeLoop->finalize(leafGridView.comm());
-
- ////////////////////////////////////////////////////////////
- // finalize, print dumux message to say goodbye
- ////////////////////////////////////////////////////////////
-
- // print dumux end message
- if (mpiHelper.rank() == 0)
- {
- Parameters::print();
- DumuxMessage::print(/*firstCall=*/false);
- }
-
- return 0;
-} // end main
-catch (Dumux::ParameterException &e)
-{
- std::cerr << std::endl << e << " ---> Abort!" << std::endl;
- return 1;
-}
-catch (Dune::DGFException & e)
-{
- std::cerr << "DGF exception thrown (" << e <<
- "). Most likely, the DGF file name is wrong "
- "or the DGF file is corrupted, "
- "e.g. missing hash at end of file or wrong number (dimensions) of entries."
- << " ---> Abort!" << std::endl;
- return 2;
-}
-catch (Dune::Exception &e)
-{
- std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
- return 3;
-}
-catch (...)
-{
- std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
- return 4;
-}
--
GitLab