From 62ff7c8194369cb17934f3112de4ecd0b0f79387 Mon Sep 17 00:00:00 2001
From: Melanie Lipp <melanie.lipp@iws.uni-stuttgart.de>
Date: Thu, 16 Sep 2021 15:45:57 +0200
Subject: [PATCH 1/3] [appl][angeli] Fresh copy from dumux-adaptivestaggered.
---
appl/freeflow/navierstokes/angeli/main.cc | 189 ++++++-----
.../freeflow/navierstokes/angeli/params.input | 23 +-
appl/freeflow/navierstokes/angeli/problem.hh | 311 +++++-------------
3 files changed, 189 insertions(+), 334 deletions(-)
diff --git a/appl/freeflow/navierstokes/angeli/main.cc b/appl/freeflow/navierstokes/angeli/main.cc
index 7714809..cfc9c13 100644
--- a/appl/freeflow/navierstokes/angeli/main.cc
+++ b/appl/freeflow/navierstokes/angeli/main.cc
@@ -22,11 +22,7 @@
* \brief Test for the instationary staggered grid Navier-Stokes model with analytical solution (Angeli et al., 2017)
*/
-bool doFirstOrderLocalTruncErrorGlobalRefinement = true;
-
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
-#pragma GCC diagnostic ignored "-Wdeprecated-copy"
+bool doFirstOrderLocalTruncErrorGlobalRefinement = false;
#include <config.h>
@@ -48,16 +44,18 @@ bool doFirstOrderLocalTruncErrorGlobalRefinement = true;
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
-#include <dumux/assembly/cvdstaggeredrefinedfvassembler.hh>
+#include <dumux/assembly/staggeredrefinedfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/mystaggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
-#include <dumux/io/readdofbasedresult.hh>
-#include <dumux/freeflow/navierstokes/staggered/cvdresidualCalc.hh>
+#include <dumux/freeflow/navierstokes/staggered/residualCalc.hh>
+#include <dumux/freeflow/navierstokes/errors.hh>
+
+#include <dumux/io/outputFacility.hh>
#include "problem.hh"
@@ -110,6 +108,11 @@ int main(int argc, char** argv) try
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
+ bool adapt = getParam<bool>("Adaptivity.Adapt", false);
+ bool doFirstOrderLocalTruncError = getParam<bool>("Numerics.DoFirstOrderLocalTruncError", false);
+ //when adapt is true refinement is local, not only global
+ doFirstOrderLocalTruncErrorGlobalRefinement = !adapt && doFirstOrderLocalTruncError;
+
using HelpingGridManager = Dumux::GridManager<GetPropType<TypeTag, Properties::HelpingGrid>>;
HelpingGridManager helpingGridManager;
helpingGridManager.init("Helper");
@@ -129,35 +132,41 @@ int main(int argc, char** argv) try
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
- bool adapt = getParam<bool>("Adaptivity.Adapt", false);
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
if (adapt)
{
- gridManager.grid().preAdapt();
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- for (const auto& element : elements(leafGridView))
- {
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- GlobalPosition pos = element.geometry().center();
+ Scalar leftX = getParam<Scalar>("Adaptivity.LeftX");
+ Scalar rightX = getParam<Scalar>("Adaptivity.RightX");
+ Scalar lowerY = getParam<Scalar>("Adaptivity.LowerY");
+ Scalar upperY = getParam<Scalar>("Adaptivity.UpperY");
- auto x = pos[0];
- auto y = pos[1];
+ gridManager.grid().preAdapt();
- if((x > 0.4 && x < 0.6) && (y > 0.4 && y < 0.6))
+ for (const auto& element : elements(leafGridView))
{
- gridManager.grid().mark( 1, element);
+ GlobalPosition pos = element.geometry().center();
+
+ auto x = pos[0];
+ auto y = pos[1];
+
+ if((leftX < x) && (x < rightX) && (y > lowerY) && (y < upperY))
+ {
+ gridManager.grid().mark( 1, element);
+ }
}
- }
- gridManager.grid().adapt();
- gridManager.grid().postAdapt();
+ gridManager.grid().adapt();
+ gridManager.grid().postAdapt();
}
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
- gridGeometry->update();
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -203,20 +212,15 @@ int main(int argc, char** argv) try
vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
-
- if (getParam<bool>("Problem.ReadDofBasedPressure", false))
- {
- using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
- CellCenterSolutionVector IBAMRpres;
- IBAMRpres.resize(gridGeometry->numCellCenterDofs());
- readDofBasedResult(IBAMRpres, gridGeometry, "pressure");
- vtkWriter.addField(IBAMRpres, "IBAMRpres");
- }
-
vtkWriter.write(0.0);
+ using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
+ CellCenterSolutionVector numericalCCResidualWithAnalytical;
+ numericalCCResidualWithAnalytical.resize(numDofsCellCenter);
+ vtkWriter.addField(numericalCCResidualWithAnalytical, "continuityResidual");
+
// the assembler with time loop for instationary problem
- using Assembler = CVDStaggeredRefinedFVAssembler<TypeTag, DiffMethod::numeric>;
+ using Assembler = StaggeredRefinedFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
@@ -227,83 +231,91 @@ int main(int argc, char** argv) try
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
- if (timeLoop->time() != 0)
- DUNE_THROW(Dune::InvalidStateException, "Was expecting to start with time zero.");
-
- //here, the temporal starting value is for t=0, the iteration starting value is the analytical solution for t=deltaT
- problem->createAnalyticalSolution(timeLoop->timeStepSize());
- SolutionVector analyticalFirstTimeStep;
- analyticalFirstTimeStep[GridGeometry::faceIdx()] = problem->getAnalyticalFaceSolutionVector();
- analyticalFirstTimeStep[GridGeometry::cellCenterIdx()] = problem->getAnalyticalPressureSolution();
- gridVariables->update(analyticalFirstTimeStep);
- nonLinearSolver.assembleLinearSystem(analyticalFirstTimeStep);
- const auto faceResidualWithAnalyticalSol = (*assembler).residual()[Dune::index_constant<0>()];
- const auto ccResidualWithAnalyticalSol = (*assembler).residual()[Dune::index_constant<1>()];
-// std::string filenamevars0 = "analyticalsol0";
-// std::ofstream s;
-// s.open(filenamevars0 + ".csv");
-// writeVectorToCSV(s, analyticalFirstTimeStep[Dune::index_constant<0>()]);
-// s.close();
-// std::string filenamevars1 = "analyticalsol1";
-// s.open(filenamevars1 + ".csv");
-// writeVectorToCSV(s, analyticalFirstTimeStep[Dune::index_constant<1>()]);
-// s.close();
- gridVariables->update(x);
+ using HostGrid = Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming >;
+ using HostGridManager = GridManager<HostGrid>;
+
+ std::array<std::string, 4> paramGroups = {"finexCVs", "fineyCVs", "coarsexCVs", "coarseyCVs"}; //xfine...
+ std::array<HostGridManager, 4> cvGridManagers = {};
+ std::array<std::map<GlobalPosition, std::vector<unsigned int>, ContainerCmpClass<GlobalPosition>>,4> cvCenterScvfIndicesMaps;
+ CVOutputFacility<TypeTag> cvOutputFacility(problem);
+ cvOutputFacility.iniCVGridManagers(cvGridManagers, cvCenterScvfIndicesMaps, paramGroups);
+ std::vector<std::pair<unsigned int, Scalar>> pvdFileInfos;
+
+ NavierStokesRefinedErrors errors(problem, x, 0.0, 0);
+ NavierStokesRefinedErrorCSVWriter errorCSVWriter(problem);
+ unsigned int timeIndex = 0;
// time loop
- const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
- int timeIndex = 0;
timeLoop->start(); do
{
+ Scalar actualCurrentTime = timeLoop->time()+timeLoop->timeStepSize();
+
// 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();
+ //output
+ if (getParam<bool>("Problem.PrintErrors", true))
+ {
+ errors.update(x, actualCurrentTime, timeIndex);
+ errorCSVWriter.printErrors(errors);
+ }
+
+ using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
+ FaceSolutionVector numericalFaceResidualWithAnalytical;
+
+ std::optional<std::array<std::vector<Scalar>, 4>> optionalPerfectInterpolationFaceResiduals;
+ std::array<std::vector<Scalar>, 4> perfectInterpolationFaceResiduals;
+ optionalPerfectInterpolationFaceResiduals.emplace(perfectInterpolationFaceResiduals);
+
+ ResidualCalc<TypeTag> residualCalc(problem);
+ residualCalc.calcResidualsInstationary(timeLoop, gridGeometry, cvGridManagers, cvCenterScvfIndicesMaps, numericalCCResidualWithAnalytical, numericalFaceResidualWithAnalytical, optionalPerfectInterpolationFaceResiduals);
- if(printL2Error)
+ bool readDofBasedValues = getParam<bool>("Problem.ReadDofBasedValues", false) && scalarCmp(getParam<Scalar>("Problem.TimeDofBasedValues"), actualCurrentTime);
+ if (readDofBasedValues)
{
- using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-
- const auto l2error = L2Error::calculateL2Error(*problem, x);
- const int numCellCenterDofs = (*gridGeometry).numCellCenterDofs();
- const int numFaceDofs = (*gridGeometry).numFaceDofs();
- std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
- << std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
- << std::scientific
- << "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
- << ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
- << ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
- << std::endl;
+ CellCenterSolutionVector IBAMRpres;
+ IBAMRpres.resize(gridGeometry->numCellCenterDofs());
+ ReadDofBasedResult<TypeTag> reader(problem);
+ reader.readDofBasedResult(IBAMRpres, "pressure");
+ vtkWriter.addField(IBAMRpres, "IBAMRpres");
}
- // advance to the time loop to the next step
- timeLoop->advanceTimeStep();
- problem->updateTime(timeLoop->time());
+ problem->createAnalyticalSolution(actualCurrentTime);
+ vtkWriter.write(actualCurrentTime);
- problem->createAnalyticalSolution(timeLoop->time());
+ if (readDofBasedValues)
+ vtkWriter.removeField();
- // write vtk output
+ cvOutputFacility.onCVsOutputResidualsAndPrimVars(
+ gridGeometry,
+ cvGridManagers,
+ cvCenterScvfIndicesMaps,
+ paramGroups,
+ x,
+ timeIndex,/*0 for stationary*/
+ numericalFaceResidualWithAnalytical,
+ optionalPerfectInterpolationFaceResiduals,
+ readDofBasedValues);
- vtkWriter.write(timeLoop->time());
+ pvdFileInfos.push_back(std::make_pair(timeIndex, actualCurrentTime));
- // report statistics of this time step
+ // prepare next time step
+ xOld = x;
+ gridVariables->advanceTimeStep();
+ timeLoop->advanceTimeStep();
+ problem->updateTime(timeLoop->time());
timeLoop->reportTimeStep();
-
- // set new dt as suggested by newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
problem->updateTimeStepSize(timeLoop->timeStepSize());
++timeIndex;
} while (!timeLoop->finished());
+ cvOutputFacility.printPvdFiles(paramGroups, pvdFileInfos);
+
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
@@ -343,4 +355,3 @@ catch (...)
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
-#pragma GCC diagnostic pop
diff --git a/appl/freeflow/navierstokes/angeli/params.input b/appl/freeflow/navierstokes/angeli/params.input
index f7e7dce..f2099f4 100644
--- a/appl/freeflow/navierstokes/angeli/params.input
+++ b/appl/freeflow/navierstokes/angeli/params.input
@@ -28,10 +28,13 @@ File = grid.dgf
[Problem]
Name = test_angeli # name passed to the output routines
EnableGravity = false
-PrintL2Error = false
+PrintErrors = true
EnableInertiaTerms = true
-AveragedAnalyticalSolution = false
+ReadDofBasedPressure = false
ReadDofBasedVelocity = false
+UseScvfLineAveraging = true
+UseContiSourceAveraging = true
+QuadratureOrder = 6 #among others for continuity-residual to become closer to zero
[Component]
LiquidDensity = 1
@@ -43,19 +46,21 @@ NumericDifference.BaseEpsilon = 1e-8
[Newton]
MaxSteps = 20
-TargetSteps = 20
MaxRelativeShift = 1e-5
-ExitAfterFirstIteration = false
-OutputMatrices = false
[Vtk]
WriteFaceData = false
+[Numerics]
+DoFirstOrderLocalTruncError = false
+
[Adaptivity]
-Adapt = false
-Conservation = false
-LinearInterpolation = false
-DebugDof = -1
+Conservation = true
+Adapt = true
+LeftX = .4
+RightX = .6
+LowerY = .4
+UpperY = .6
[VelocityReadIn]
NumFileEntries = 260
diff --git a/appl/freeflow/navierstokes/angeli/problem.hh b/appl/freeflow/navierstokes/angeli/problem.hh
index 2ee3874..77037c1 100644
--- a/appl/freeflow/navierstokes/angeli/problem.hh
+++ b/appl/freeflow/navierstokes/angeli/problem.hh
@@ -21,44 +21,39 @@
* \ingroup NavierStokesTests
* \brief Test for the instationary staggered grid Navier-Stokes model with analytical solution (Angeli et al., 2017)
*/
-#ifndef DUMUX_CVD_ANGELI_TEST_PROBLEM_HH
-#define DUMUX_CVD_ANGELI_TEST_PROBLEM_HH
+#ifndef DUMUX_ANGELI_TEST_PROBLEM_HH
+#define DUMUX_ANGELI_TEST_PROBLEM_HH
#include <dune/alugrid/grid.hh>
#include <dune/grid/yaspgrid.hh>
-#include <dumux/assembly/cvdstaggeredrefinedlocalresidual.hh>
+#include <dumux/assembly/staggeredrefinedlocalresidual.hh>
+
+#include <dumux/common/numeqvector.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/model.hh>
-#include <dumux/freeflow/navierstokes/cvdproblem.hh>
-#include <dumux/freeflow/navierstokes/cvdlocalresidual.hh>
-#include <dumux/freeflow/navierstokes/cvdfluxvariableswrapper.hh>
+#include <dumux/freeflow/navierstokes/myproblem.hh>
+#include <dumux/freeflow/navierstokes/mylocalresidual.hh>
+#include <dumux/freeflow/navierstokes/myfluxvariables.hh>
#include <dumux/freeflow/navierstokes/myvolumevariables.hh>
#include <dumux/freeflow/navierstokes/myiofields.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/staggered/myfacesolution.hh>
+#include <dumux/discretization/staggered/mygridvariables.hh>
+#include <dumux/discretization/staggered/myfvgridgeometry.hh>
#include <dumux/discretization/staggered/mygridfluxvariablescache.hh>
+#include <dumux/discretization/staggered/freeflow/fvgridgeometrytraits1.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/discretization/staggered/freeflow/myfacevariables.hh>
-#include <dumux/discretization/staggered/freeflow/cvdlocalfacevariables.hh>
#include <dumux/discretization/staggered/freeflow/elementvolumevariables1.hh>
#include <dumux/discretization/staggered/freeflow/myvelocityoutput.hh>
#include <dumux/discretization/staggered/freeflow/mygridvolumevariables.hh>
-#include <dumux/discretization/staggered/cvdelementfacesolution.hh>
-#include <dumux/discretization/staggered/cvdgridvariables.hh>
-#include <dumux/discretization/staggered/cvdlocalgridfacevariables.hh>
-#include <dumux/discretization/staggered/cvdfvgridgeometry.hh>
-#include <dumux/discretization/staggered/freeflow/cvdfvgridgeometrytraits.hh>
-
-
-#include <appl/freeflow/navierstokes/l2error.hh>
-
namespace Dumux
{
template <class TypeTag>
@@ -120,15 +115,6 @@ public:
using type = MyStaggeredFaceVariables<ModelTraits, FacePrimaryVariables, GridView::dimension, upwindSchemeOrder>;
};
-template<class TypeTag>
-struct LocalFaceVariables<TypeTag, TTag::AngeliTestTypeTag>
-{
-private:
- using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
-public:
- using type = CVDStaggeredLocalFaceVariables<FacePrimaryVariables>;
-};
-
//! The velocity output
template<class TypeTag>
struct VelocityOutput<TypeTag, TTag::AngeliTestTypeTag>
@@ -192,28 +178,6 @@ public:
using type = Dumux::MyStaggeredFaceSolution<FaceSolutionVector>;
};
-
-template<class TypeTag>
-struct ElementFaceSolution<TypeTag, TTag::AngeliTestTypeTag>
-{
-private:
- using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
-public:
- using type = Dumux::CVDElementFaceSolution<FaceSolutionVector>;
-};
-
-//! Set the GridLocalFaceVariables
-template<class TypeTag>
-struct GridLocalFaceVariables<TypeTag, TTag::AngeliTestTypeTag>
-{
- private:
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using LocalFaceVariables = GetPropType<TypeTag, Properties::LocalFaceVariables>;
- static constexpr auto enableCache = getPropValue<TypeTag, Properties::EnableGridFaceVariablesCache>();
- public:
- using type = CVDStaggeredLocalGridFaceVariables<Problem, LocalFaceVariables, enableCache>;
-};
-
//! Set the grid variables (volume, flux and face variables)
template<class TypeTag>
struct GridVariables<TypeTag, TTag::AngeliTestTypeTag>
@@ -223,9 +187,8 @@ private:
using GVV = GetPropType<TypeTag, Properties::GridVolumeVariables>;
using GFVC = GetPropType<TypeTag, Properties::GridFluxVariablesCache>;
using GFV = GetPropType<TypeTag, Properties::GridFaceVariables>;
- using GLFV = GetPropType<TypeTag, Properties::GridLocalFaceVariables>;
public:
- using type = CVDStaggeredGridVariables<GG, GVV, GFVC, GFV, GLFV>;
+ using type = MyStaggeredGridVariables<GG, GVV, GFVC, GFV>;
};
//! The default fv grid geometry
@@ -236,9 +199,9 @@ private:
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
- using Traits = CVDStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
+ using Traits = MyStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
public:
- using type = CVDStaggeredGridGeometry<GridView, enableCache, Traits>;
+ using type = MyStaggeredGridGeometry<GridView, enableCache, Traits>;
};
//! Set the volume variables property
@@ -299,9 +262,9 @@ private:
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming >::LeafGridView;
- using Traits = CVDStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
+ using Traits = MyStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
public:
- using type = CVDStaggeredGridGeometry<GridView, enableCache, Traits>;
+ using type = MyStaggeredGridGeometry<GridView, enableCache, Traits>;
};
}
@@ -319,9 +282,8 @@ class AngeliTestProblem : public MyNavierStokesProblem<TypeTag>
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
@@ -332,8 +294,12 @@ class AngeliTestProblem : public MyNavierStokesProblem<TypeTag>
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GetPropType<TypeTag, Properties::GridGeometry>::SubControlVolumeFace;
public:
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
AngeliTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry), eps_(1e-6)
{
@@ -433,127 +399,62 @@ public:
/*!
- * \brief Return dirichlet boundary values at a given position
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume face (velocities)
*
- * \param globalPos The global position
- */
- PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
- {
- // use the values of the analytical solution
- return this->instationaryAnalyticalSolution(globalPos, time_+timeStepSize_);
- }
-
- /*!
- * \brief Return the analytical solution of the problem at a given position
+ * \param element The finite element
+ * \param scvf the sub control volume face
*
- * \param globalPos The global position
+ * Concerning the usage of averaged velocity, see the explanation of the initial function. The average of the nondirection velocity along an scvf is required in the context of dirichlet boundary conditions for the functions getParallelVelocityFromBoundary_ and getParallelVelocityFromOtherBoundary_ within myfluxvariables.hh. There dirichlet is called for ghost faces built from the normalFace but the value then is taken in the the direction scvf.directionIndex, which is along that normalFace.
*/
- PrimaryVariables instationaryAnalyticalSolutionAtPos(const GlobalPosition& globalPos, const Scalar time) const
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
- const Scalar x = globalPos[0];
- const Scalar y = globalPos[1];
- const Scalar t = time;
+ Scalar time = time_+timeStepSize_;
- PrimaryVariables values = {};
+ PrimaryVariables priVars(0.0);
- for (unsigned int j = 0; j < 2; ++j)
- {
- for (unsigned int i = 0; i < values.size(); ++i)
- {
- values[i] += xFactorAnalyticalSolutionAtPos(x,t)[j][i]*yFactorAnalyticalSolutionAtPos(y,t)[j][i];
- }
- }
+ unsigned int dirIdx = scvf.directionIndex();
+ priVars[Indices::velocity(dirIdx)] = this->instationaryAnalyticalVelocitySolutionAtPos(scvf, time);
+ unsigned int nonDirIdx = (dirIdx == 0)?1:0;
+ priVars[Indices::velocity(nonDirIdx)] = this->instationaryAnalyticalNonDirVelocitySolutionAtPos(scvf, time);
- return values;
+ return priVars;
}
/*!
- * \brief Returns the analytical solution of the problem at a given position.
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume face (pressure)
*
- * \param globalPos The global position
+ * \param element The finite element
+ * \param scv the sub control volume
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
{
- return this->instationaryAnalyticalSolution(globalPos, time()+timeStepSize_);
+ Scalar time = time_+timeStepSize_;
+
+ PrimaryVariables priVars(0.0);
+ priVars[Indices::pressureIdx] = this->instationaryAnalyticalPressureSolutionAtPos(scv, time);
+ return priVars;
}
/*!
- * \brief Returns the analytical solution of the problem at a given position.
+ * \brief Return the analytical solution of the problem at a given position
*
* \param globalPos The global position
*/
- PrimaryVariables analyticalSolutionAtPos(const GlobalPosition& globalPos) const
- {
- return this->instationaryAnalyticalSolutionAtPos(globalPos, time()+timeStepSize_);
- }
-
-
- std::array<PrimaryVariables, 2> xFactorAnalyticalSolutionAtPos(Scalar x, Scalar t) const
- {
- PrimaryVariables values1;
- values1[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * 4.0 * M_PI * M_PI * std::cos(2.0 * M_PI * x)*rho_;
- values1[Indices::velocityXIdx] = f_(t) * f2_(x) ;
- values1[Indices::velocityYIdx] = - f_(t) * df2_(x);
-
- PrimaryVariables values2 = {};
- values2[Indices::pressureIdx] = 1.;
-
- std::array<PrimaryVariables, 2> retPair;
- retPair[0] = values1;
- retPair[1] = values2;
-
- return retPair;
- }
-
- std::array<PrimaryVariables, 2> yFactorAnalyticalSolutionAtPos(Scalar y, Scalar t) const
- {
- PrimaryVariables values1;
- values1[Indices::pressureIdx] = 1.;
- values1[Indices::velocityXIdx] = df3_(y);
- values1[Indices::velocityYIdx] = f3_(y);
-
- PrimaryVariables values2 = {};
- values2[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * M_PI * M_PI * std::cos(4.0 * M_PI * y)*rho_;
-
- std::array<PrimaryVariables, 2> retPair;
- retPair[0] = values1;
- retPair[1] = values2;
-
- return retPair;
- }
-
- std::array<PrimaryVariables, 2> xFactorAnalyticalSolutionAntiderivativeAtPos(Scalar x, Scalar t) const
- {
- PrimaryVariables values1;
- values1[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * 2.0 * M_PI * std::sin(2.0 * M_PI * x)*rho_;
- values1[Indices::velocityXIdx] = f_(t) * intf2_(x) ;
- values1[Indices::velocityYIdx] = - f_(t) * f2_(x);
-
- PrimaryVariables values2 = {};
- values2[Indices::pressureIdx] = x;
-
- std::array<PrimaryVariables, 2> retPair;
- retPair[0] = values1;
- retPair[1] = values2;
-
- return retPair;
- }
-
- std::array<PrimaryVariables, 2> yFactorAnalyticalSolutionAntiderivativeAtPos(Scalar y, Scalar t) const
+ PrimaryVariables instationaryAnalyticalSolutionAtPos(const GlobalPosition& globalPos, const Scalar time) const
{
- PrimaryVariables values1;
- values1[Indices::pressureIdx] = y;
- values1[Indices::velocityXIdx] = f3_(y);
- values1[Indices::velocityYIdx] = intf3_(y);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ const Scalar t = time;
- PrimaryVariables values2 = {};
- values2[Indices::pressureIdx] = - 1./16 * f_(t) * f_(t) * M_PI * std::sin(4.0 * M_PI * y)*rho_;
+ PrimaryVariables values = {};
- std::array<PrimaryVariables, 2> retPair;
- retPair[0] = values1;
- retPair[1] = values2;
+ values[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * M_PI * M_PI *rho_ * (4.0 * std::cos(2.0 * M_PI * x) - std::cos(4.0 * M_PI * y));
+ values[Indices::velocityXIdx] = f_(t) * f2_(x)*df3_(y);
+ values[Indices::velocityYIdx] = - f_(t) * df2_(x)* f3_(y);
- return retPair;
+ return values;
}
// \}
@@ -563,16 +464,32 @@ public:
*/
// \{
- /*!
- * \brief Evaluate the initial value for a control volume.
+ /*!
+ * \brief Evaluates the initial value for a sub control volume face (velocities)
*
- * \param globalPos The global position
+ * Simply assigning the value of the analytical solution at the face center
+ * gives a discrete solution that is not divergence-free. For small initial
+ * time steps, this has a negative impact on the pressure solution
+ * after the first time step. The flag UseVelocityAveraging triggers the
+ * function averagedVelocity_ which uses a higher order quadrature formula to
+ * bring the discrete solution sufficiently close to being divergence-free.
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ PrimaryVariables initial(const SubControlVolumeFace& scvf) const
{
- return this->instationaryAnalyticalSolution(globalPos, 0);
+ PrimaryVariables priVars(0.0);
+ priVars[Indices::velocity(scvf.directionIndex())] = this->instationaryAnalyticalVelocitySolutionAtPos(scvf, 0.);
+ return priVars;
}
+ /*!
+ * \brief Evaluates the initial value for a control volume (pressure)
+ */
+ PrimaryVariables initial(const SubControlVolume& scv) const
+ {
+ PrimaryVariables priVars(0.0);
+ priVars[Indices::pressureIdx] = this->instationaryAnalyticalPressureSolutionAtPos(scv, 0.);
+ return priVars;
+ }
/*!
* \brief Updates the time
*/
@@ -589,79 +506,11 @@ public:
timeStepSize_ = timeStepSize;
}
- /*!
- * \brief Returns the analytical solution for the pressure
- */
- auto& getAnalyticalPressureSolution() const
- {
- return analyticalPressure_;
- }
-
- /*!
- * \brief Returns the analytical solution for the velocity
- */
- auto& getAnalyticalVelocitySolution() const
- {
- return analyticalVelocity_;
- }
-
- auto& getAnalyticalFaceSolutionVector() const
- {
- return analyticalFaceSolutionVector_;
- }
-
- /*!
- * \brief Returns the analytical solution for the velocity at the faces
- */
- auto& getAnalyticalVelocitySolutionOnFace() const
- {
- return analyticalVelocityOnFace_;
- }
-
Scalar time() const
{
return time_;
}
- /*!
- * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
- */
- void createAnalyticalSolution(Scalar t)
- {
- analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
- analyticalFaceSolutionVector_.resize(this->gridGeometry().numFaceDofs());
-
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- auto fvGeometry = localView(this->gridGeometry());
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- auto ccDofIdx = scv.dofIndex();
- auto ccDofPosition = scv.dofPosition();
- auto analyticalSolutionAtCc = this->instationaryAnalyticalSolution(ccDofPosition, t);
-
- // velocities on faces
- for (auto&& scvf : scvfs(fvGeometry))
- {
- const auto faceDofIdx = scvf.dofIndex();
- const auto faceDofPosition = scvf.center();
- const auto dirIdx = scvf.directionIndex();
- const auto analyticalSolutionAtFace = this->instationaryAnalyticalSolution(faceDofPosition, t);
- analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
- analyticalFaceSolutionVector_[faceDofIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
- }
-
- analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
-
- for(int dirIdx = 0; dirIdx < ModelTraits::dim(); ++dirIdx)
- analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
- }
- }
- }
-
private:
Scalar f_(Scalar t) const
{ return std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * t); }
@@ -669,18 +518,12 @@ private:
Scalar df_(Scalar t) const
{ return - 5.0 * kinematicViscosity_ * M_PI * M_PI * f_(t); }
- Scalar intf2_ (Scalar x) const
- { return std::sin(M_PI * x)/M_PI; }
-
Scalar f2_ (Scalar x) const
{ return std::cos(M_PI * x); }
Scalar df2_ (Scalar x) const
{ return -M_PI * std::sin(M_PI * x); }
- Scalar intf3_ ( Scalar x) const
- { return std::sin(2.0 * M_PI * x)/(2.0 * M_PI); }
-
Scalar f3_ ( Scalar x) const
{ return std::cos(2.0 * M_PI * x); }
@@ -697,10 +540,6 @@ private:
Scalar kinematicViscosity_;
Scalar rho_;
- CellCenterSolutionVector analyticalPressure_;
- std::vector<VelocityVector> analyticalVelocity_;
- std::vector<VelocityVector> analyticalVelocityOnFace_;
- FaceSolutionVector analyticalFaceSolutionVector_;
Scalar time_ = 0;
Scalar timeStepSize_ = 0;
--
GitLab
From 27b72b25f4887b5fd14dea34323deaa63db8315f Mon Sep 17 00:00:00 2001
From: Melanie Lipp <melanie.lipp@iws.uni-stuttgart.de>
Date: Thu, 16 Sep 2021 16:26:34 +0200
Subject: [PATCH 2/3] [appl][angeli] Adapt to the needs of the cvd refinement.
---
appl/freeflow/navierstokes/angeli/main.cc | 30 +++++----
.../freeflow/navierstokes/angeli/params.input | 7 +-
appl/freeflow/navierstokes/angeli/problem.hh | 64 +++++++++++++++----
3 files changed, 74 insertions(+), 27 deletions(-)
diff --git a/appl/freeflow/navierstokes/angeli/main.cc b/appl/freeflow/navierstokes/angeli/main.cc
index cfc9c13..443fabf 100644
--- a/appl/freeflow/navierstokes/angeli/main.cc
+++ b/appl/freeflow/navierstokes/angeli/main.cc
@@ -24,6 +24,10 @@
bool doFirstOrderLocalTruncErrorGlobalRefinement = false;
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
+#pragma GCC diagnostic ignored "-Wdeprecated-copy"
+
#include <config.h>
#include <ctime>
@@ -37,6 +41,7 @@ bool doFirstOrderLocalTruncErrorGlobalRefinement = false;
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/coursevelocityproperties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
@@ -44,18 +49,19 @@ bool doFirstOrderLocalTruncErrorGlobalRefinement = false;
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
-#include <dumux/assembly/staggeredrefinedfvassembler.hh>
+#include <dumux/assembly/cvdstaggeredrefinedfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/mystaggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
+#include <dumux/io/grid/cvdcontrolvolumegrids.hh>
-#include <dumux/freeflow/navierstokes/staggered/residualCalc.hh>
+#include <dumux/freeflow/navierstokes/staggered/cvdresidualCalc.hh>
#include <dumux/freeflow/navierstokes/errors.hh>
-#include <dumux/io/outputFacility.hh>
+#include <dumux/io/cvdoutputFacility.hh>
#include "problem.hh"
@@ -220,7 +226,7 @@ int main(int argc, char** argv) try
vtkWriter.addField(numericalCCResidualWithAnalytical, "continuityResidual");
// the assembler with time loop for instationary problem
- using Assembler = StaggeredRefinedFVAssembler<TypeTag, DiffMethod::numeric>;
+ using Assembler = CVDStaggeredRefinedFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
@@ -234,14 +240,15 @@ int main(int argc, char** argv) try
using HostGrid = Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming >;
using HostGridManager = GridManager<HostGrid>;
- std::array<std::string, 4> paramGroups = {"finexCVs", "fineyCVs", "coarsexCVs", "coarseyCVs"}; //xfine...
- std::array<HostGridManager, 4> cvGridManagers = {};
- std::array<std::map<GlobalPosition, std::vector<unsigned int>, ContainerCmpClass<GlobalPosition>>,4> cvCenterScvfIndicesMaps;
- CVOutputFacility<TypeTag> cvOutputFacility(problem);
+ std::array<std::string, 2> paramGroups = {"coarsexCVs", "coarseyCVs"};
+ std::array<HostGridManager, 2> cvGridManagers = {};
+ std::array<std::map<GlobalPosition, std::vector<unsigned int>, ContainerCmpClass<GlobalPosition>>,2> cvCenterScvfIndicesMaps;
+ CVDCVOutputFacility<TypeTag> cvOutputFacility(problem);
cvOutputFacility.iniCVGridManagers(cvGridManagers, cvCenterScvfIndicesMaps, paramGroups);
std::vector<std::pair<unsigned int, Scalar>> pvdFileInfos;
NavierStokesRefinedErrors errors(problem, x, 0.0, 0);
+ errors.setCoarseVelocityDiscretization(true);
NavierStokesRefinedErrorCSVWriter errorCSVWriter(problem);
unsigned int timeIndex = 0;
@@ -266,12 +273,10 @@ int main(int argc, char** argv) try
using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
FaceSolutionVector numericalFaceResidualWithAnalytical;
- std::optional<std::array<std::vector<Scalar>, 4>> optionalPerfectInterpolationFaceResiduals;
- std::array<std::vector<Scalar>, 4> perfectInterpolationFaceResiduals;
- optionalPerfectInterpolationFaceResiduals.emplace(perfectInterpolationFaceResiduals);
+ std::optional<std::array<std::vector<Scalar>, 2>> optionalPerfectInterpolationFaceResiduals;
ResidualCalc<TypeTag> residualCalc(problem);
- residualCalc.calcResidualsInstationary(timeLoop, gridGeometry, cvGridManagers, cvCenterScvfIndicesMaps, numericalCCResidualWithAnalytical, numericalFaceResidualWithAnalytical, optionalPerfectInterpolationFaceResiduals);
+ residualCalc.calcResidualsInstationary(timeLoop, gridGeometry, numericalCCResidualWithAnalytical, numericalFaceResidualWithAnalytical);
bool readDofBasedValues = getParam<bool>("Problem.ReadDofBasedValues", false) && scalarCmp(getParam<Scalar>("Problem.TimeDofBasedValues"), actualCurrentTime);
if (readDofBasedValues)
@@ -355,3 +360,4 @@ catch (...)
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
+#pragma GCC diagnostic pop
diff --git a/appl/freeflow/navierstokes/angeli/params.input b/appl/freeflow/navierstokes/angeli/params.input
index f2099f4..b0c3753 100644
--- a/appl/freeflow/navierstokes/angeli/params.input
+++ b/appl/freeflow/navierstokes/angeli/params.input
@@ -46,7 +46,10 @@ NumericDifference.BaseEpsilon = 1e-8
[Newton]
MaxSteps = 20
+TargetSteps = 20
MaxRelativeShift = 1e-5
+ExitAfterFirstIteration = false
+OutputMatrices = false
[Vtk]
WriteFaceData = false
@@ -55,12 +58,14 @@ WriteFaceData = false
DoFirstOrderLocalTruncError = false
[Adaptivity]
-Conservation = true
+Conservation = false
Adapt = true
LeftX = .4
RightX = .6
LowerY = .4
UpperY = .6
+LinearInterpolation = false
+DebugDof = -1
[VelocityReadIn]
NumFileEntries = 260
diff --git a/appl/freeflow/navierstokes/angeli/problem.hh b/appl/freeflow/navierstokes/angeli/problem.hh
index 77037c1..5f05e7c 100644
--- a/appl/freeflow/navierstokes/angeli/problem.hh
+++ b/appl/freeflow/navierstokes/angeli/problem.hh
@@ -21,13 +21,13 @@
* \ingroup NavierStokesTests
* \brief Test for the instationary staggered grid Navier-Stokes model with analytical solution (Angeli et al., 2017)
*/
-#ifndef DUMUX_ANGELI_TEST_PROBLEM_HH
-#define DUMUX_ANGELI_TEST_PROBLEM_HH
+#ifndef DUMUX_CVD_ANGELI_TEST_PROBLEM_HH
+#define DUMUX_CVD_ANGELI_TEST_PROBLEM_HH
#include <dune/alugrid/grid.hh>
#include <dune/grid/yaspgrid.hh>
-#include <dumux/assembly/staggeredrefinedlocalresidual.hh>
+#include <dumux/assembly/cvdstaggeredrefinedlocalresidual.hh>
#include <dumux/common/numeqvector.hh>
@@ -36,24 +36,28 @@
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/model.hh>
-#include <dumux/freeflow/navierstokes/myproblem.hh>
-#include <dumux/freeflow/navierstokes/mylocalresidual.hh>
-#include <dumux/freeflow/navierstokes/myfluxvariables.hh>
+#include <dumux/freeflow/navierstokes/cvdproblem.hh>
+#include <dumux/freeflow/navierstokes/cvdlocalresidual.hh>
+#include <dumux/freeflow/navierstokes/cvdfluxvariableswrapper.hh>
#include <dumux/freeflow/navierstokes/myvolumevariables.hh>
#include <dumux/freeflow/navierstokes/myiofields.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/staggered/myfacesolution.hh>
-#include <dumux/discretization/staggered/mygridvariables.hh>
-#include <dumux/discretization/staggered/myfvgridgeometry.hh>
#include <dumux/discretization/staggered/mygridfluxvariablescache.hh>
-#include <dumux/discretization/staggered/freeflow/fvgridgeometrytraits1.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/discretization/staggered/freeflow/myfacevariables.hh>
+#include <dumux/discretization/staggered/freeflow/cvdlocalfacevariables.hh>
#include <dumux/discretization/staggered/freeflow/elementvolumevariables1.hh>
#include <dumux/discretization/staggered/freeflow/myvelocityoutput.hh>
#include <dumux/discretization/staggered/freeflow/mygridvolumevariables.hh>
+#include <dumux/discretization/staggered/cvdelementfacesolution.hh>
+#include <dumux/discretization/staggered/cvdgridvariables.hh>
+#include <dumux/discretization/staggered/cvdlocalgridfacevariables.hh>
+#include <dumux/discretization/staggered/cvdfvgridgeometry.hh>
+#include <dumux/discretization/staggered/freeflow/cvdfvgridgeometrytraits.hh>
+
namespace Dumux
{
template <class TypeTag>
@@ -115,6 +119,15 @@ public:
using type = MyStaggeredFaceVariables<ModelTraits, FacePrimaryVariables, GridView::dimension, upwindSchemeOrder>;
};
+template<class TypeTag>
+struct LocalFaceVariables<TypeTag, TTag::AngeliTestTypeTag>
+{
+private:
+ using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
+public:
+ using type = CVDStaggeredLocalFaceVariables<FacePrimaryVariables>;
+};
+
//! The velocity output
template<class TypeTag>
struct VelocityOutput<TypeTag, TTag::AngeliTestTypeTag>
@@ -178,6 +191,28 @@ public:
using type = Dumux::MyStaggeredFaceSolution<FaceSolutionVector>;
};
+
+template<class TypeTag>
+struct ElementFaceSolution<TypeTag, TTag::AngeliTestTypeTag>
+{
+private:
+ using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
+public:
+ using type = Dumux::CVDElementFaceSolution<FaceSolutionVector>;
+};
+
+//! Set the GridLocalFaceVariables
+template<class TypeTag>
+struct GridLocalFaceVariables<TypeTag, TTag::AngeliTestTypeTag>
+{
+ private:
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using LocalFaceVariables = GetPropType<TypeTag, Properties::LocalFaceVariables>;
+ static constexpr auto enableCache = getPropValue<TypeTag, Properties::EnableGridFaceVariablesCache>();
+ public:
+ using type = CVDStaggeredLocalGridFaceVariables<Problem, LocalFaceVariables, enableCache>;
+};
+
//! Set the grid variables (volume, flux and face variables)
template<class TypeTag>
struct GridVariables<TypeTag, TTag::AngeliTestTypeTag>
@@ -187,8 +222,9 @@ private:
using GVV = GetPropType<TypeTag, Properties::GridVolumeVariables>;
using GFVC = GetPropType<TypeTag, Properties::GridFluxVariablesCache>;
using GFV = GetPropType<TypeTag, Properties::GridFaceVariables>;
+ using GLFV = GetPropType<TypeTag, Properties::GridLocalFaceVariables>;
public:
- using type = MyStaggeredGridVariables<GG, GVV, GFVC, GFV>;
+ using type = CVDStaggeredGridVariables<GG, GVV, GFVC, GFV, GLFV>;
};
//! The default fv grid geometry
@@ -199,9 +235,9 @@ private:
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
- using Traits = MyStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
+ using Traits = CVDStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
public:
- using type = MyStaggeredGridGeometry<GridView, enableCache, Traits>;
+ using type = CVDStaggeredGridGeometry<GridView, enableCache, Traits>;
};
//! Set the volume variables property
@@ -262,9 +298,9 @@ private:
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming >::LeafGridView;
- using Traits = MyStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
+ using Traits = CVDStaggeredFreeFlowDefaultGridGeometryTraits<GridView, upwindSchemeOrder>;
public:
- using type = MyStaggeredGridGeometry<GridView, enableCache, Traits>;
+ using type = CVDStaggeredGridGeometry<GridView, enableCache, Traits>;
};
}
--
GitLab
From c728ef3b032939dc289d45ef5fcbd24e87366052 Mon Sep 17 00:00:00 2001
From: Melanie Lipp <melanie.lipp@iws.uni-stuttgart.de>
Date: Thu, 16 Sep 2021 18:07:40 +0200
Subject: [PATCH 3/3] [appl][angeli] Correct sign. This is a bug that was
introduced by commit 32e0daac. With this fix I also finally retrieve the L2
errors that I gave to Caroline.
---
appl/freeflow/navierstokes/angeli/problem.hh | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/appl/freeflow/navierstokes/angeli/problem.hh b/appl/freeflow/navierstokes/angeli/problem.hh
index 5f05e7c..066486a 100644
--- a/appl/freeflow/navierstokes/angeli/problem.hh
+++ b/appl/freeflow/navierstokes/angeli/problem.hh
@@ -486,7 +486,7 @@ public:
PrimaryVariables values = {};
- values[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * M_PI * M_PI *rho_ * (4.0 * std::cos(2.0 * M_PI * x) - std::cos(4.0 * M_PI * y));
+ values[Indices::pressureIdx] = - 0.25 * f_(t) * f_(t) * M_PI * M_PI *rho_ * (4.0 * std::cos(2.0 * M_PI * x) + std::cos(4.0 * M_PI * y));
values[Indices::velocityXIdx] = f_(t) * f2_(x)*df3_(y);
values[Indices::velocityYIdx] = - f_(t) * df2_(x)* f3_(y);
--
GitLab