Commit 09218942 authored by Katharina Heck's avatar Katharina Heck
Browse files

[cleanup] delete unnecessary input

cleanup
parent c993aecc
...@@ -20,6 +20,6 @@ dumux_add_test(NAME test_ff_navierstokes_angeli_averaged ...@@ -20,6 +20,6 @@ dumux_add_test(NAME test_ff_navierstokes_angeli_averaged
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli params.input --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli params.input
-TimeLoop.TEnd 1e-6 -TimeLoop.DtInitial 1e-6 -TimeLoop.TEnd 1e-6 -TimeLoop.DtInitial 1e-6
-Problem.Name test_ff_navierstokes_angeli_averaged -Problem.Name test_ff_navierstokes_angeli_averaged
-Problem.UseQuad true") -Problem.InterpolateExactVelocity true")
dune_symlink_to_source_files(FILES "params.input") dune_symlink_to_source_files(FILES "params.input")
...@@ -32,58 +32,25 @@ ...@@ -32,58 +32,25 @@
#include <dune/common/parallel/mpihelper.hh> #include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh> #include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/common/dumuxmessage.hh> #include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh> #include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh> #include <dumux/common/properties.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/linear/seqsolverbackend.hh> #include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh> #include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh> #include <test/freeflow/navierstokes/errors.hh>
#include "properties.hh" #include "properties.hh"
template<class MomentumProblem, class MassProblem, class SolutionVector>
void printL2Error(const MomentumProblem& momentumProblem,
const MassProblem& massProblem,
const SolutionVector& x)
{
const auto velocityL2error = calculateL2Error(momentumProblem, x[Dune::index_constant<0>()]);
const auto pressureL2error = calculateL2Error(massProblem, x[Dune::index_constant<1>()]);
static constexpr auto pressureIdx = MassProblem::Indices::pressureIdx;
static constexpr auto velocityXIdx = MomentumProblem::Indices::velocityXIdx;
static constexpr auto velocityYIdx = MomentumProblem::Indices::velocityYIdx;
const auto numCellCenterDofs = massProblem.gridGeometry().numDofs();
const auto numFaceDofs = momentumProblem.gridGeometry().numDofs();
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) = " << pressureL2error.absolute[pressureIdx] << " / " << pressureL2error.relative[pressureIdx]
<< " , L2(vx) = " << velocityL2error.absolute[velocityXIdx] << " / " << velocityL2error.relative[velocityXIdx]
<< " , L2(vy) = " << velocityL2error.absolute[velocityYIdx] << " / " << velocityL2error.relative[velocityYIdx]
<< std::endl;
// write the norm into a log file
std::ofstream logFile;
logFile.open(momentumProblem.name() + ".log", std::ios::app);
logFile << "[ConvergenceTest] L2(p) = " << pressureL2error.absolute[pressureIdx] << " L2(vx) = "
<< velocityL2error.absolute[velocityXIdx] << " L2(vy) = " << velocityL2error.absolute[velocityYIdx] << std::endl;
logFile.close();
}
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
using namespace Dumux; using namespace Dumux;
...@@ -146,12 +113,10 @@ int main(int argc, char** argv) ...@@ -146,12 +113,10 @@ int main(int argc, char** argv)
timeLoop->setMaxTimeStepSize(maxDt); timeLoop->setMaxTimeStepSize(maxDt);
// the solution vector // the solution vector
constexpr auto momentumIdx = Dune::index_constant<0>(); constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = Dune::index_constant<1>(); constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector; using SolutionVector = typename Traits::SolutionVector;
SolutionVector x; SolutionVector x;
x[momentumIdx].resize(momentumGridGeometry->numDofs());
x[massIdx].resize(massGridGeometry->numDofs());
momentumProblem->applyInitialSolution(x[momentumIdx]); momentumProblem->applyInitialSolution(x[momentumIdx]);
massProblem->applyInitialSolution(x[massIdx]); massProblem->applyInitialSolution(x[massIdx]);
auto xOld = x; auto xOld = x;
...@@ -186,15 +151,6 @@ int main(int argc, char** argv) ...@@ -186,15 +151,6 @@ int main(int argc, char** argv)
std::make_tuple(momentumGridVariables, massGridVariables), std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager, timeLoop, xOld); couplingManager, timeLoop, xOld);
assembler->assembleResidual(x);
std::cout << "Initial momentum defect " << std::endl;
for (const auto s : assembler->residual()[momentumIdx])
std::cout << s << std::endl;
std::cout << "Initial mass defect " << std::endl;
for (const auto s : assembler->residual()[massIdx])
std::cout << s << std::endl;
// the linear solver // the linear solver
using LinearSolver = Dumux::UMFPackBackend; using LinearSolver = Dumux::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>(); auto linearSolver = std::make_shared<LinearSolver>();
......
...@@ -80,7 +80,7 @@ public: ...@@ -80,7 +80,7 @@ public:
{ {
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0); kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
useNeumann_ = getParam<bool>("Problem.UseNeumann", false); useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
useQuad_ = getParam<bool>("Problem.UseQuad", false); interpolateExactVelocity_ = getParam<bool>("Problem.InterpolateExactVelocity", false);
} }
/*! /*!
...@@ -144,17 +144,8 @@ public: ...@@ -144,17 +144,8 @@ public:
// If only Dirichlet BCs are set for the momentum balance, fix the pressure at some cells such that the solution is fully defined. // If only Dirichlet BCs are set for the momentum balance, fix the pressure at some cells such that the solution is fully defined.
if (!useNeumann_) if (!useNeumann_)
{ {
auto fvGeometry = localView(this->gridGeometry()); const auto fvGeometry = localView(this->gridGeometry()).bindElement(element);
fvGeometry.bindElement(element);
const auto isAtBoundary = [&](const FVElementGeometry& fvGeometry)
{
if (fvGeometry.hasBoundaryScvf()) if (fvGeometry.hasBoundaryScvf())
return true;
return false;
};
if (isAtBoundary(fvGeometry))
values.set(Indices::pressureIdx); values.set(Indices::pressureIdx);
} }
...@@ -182,7 +173,7 @@ public: ...@@ -182,7 +173,7 @@ public:
*/ */
PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{ {
if (ParentType::isMomentumProblem() && useQuad_) if (ParentType::isMomentumProblem() && interpolateExactVelocity_)
return velocityDirichlet_(scvf); return velocityDirichlet_(scvf);
else else
return analyticalSolution(scvf.center(), time_); return analyticalSolution(scvf.center(), time_);
...@@ -222,7 +213,7 @@ public: ...@@ -222,7 +213,7 @@ public:
const Scalar t = time_; const Scalar t = time_;
momentumFlux[0][0] = M_PI*M_PI *(-4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (4.0*sin(2*M_PI*y)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*x) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t); momentumFlux[0][0] = M_PI*M_PI *(-4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (4.0*sin(2*M_PI*y)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*x) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
momentumFlux[0][1] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y); momentumFlux[0][1] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y);
momentumFlux[1][0] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y); momentumFlux[1][0] = momentumFlux[0][1];
momentumFlux[1][1] = M_PI*M_PI *( 4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (sin(M_PI*x)*sin(M_PI*x)*cos(2*M_PI*y)*cos(2*M_PI*y) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t); momentumFlux[1][1] = M_PI*M_PI *( 4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (sin(M_PI*x)*sin(M_PI*x)*cos(2*M_PI*y)*cos(2*M_PI*y) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
const auto normal = scvf.unitOuterNormal(); const auto normal = scvf.unitOuterNormal();
...@@ -280,7 +271,7 @@ public: ...@@ -280,7 +271,7 @@ public:
// because the discrete mass balance equation is not fulfilled when just taking // because the discrete mass balance equation is not fulfilled when just taking
// point-wise values of the analytical solution. // point-wise values of the analytical solution.
if (!useQuad_) if (!interpolateExactVelocity_)
return analyticalSolution(scv.dofPosition(), 0.0); return analyticalSolution(scv.dofPosition(), 0.0);
// Get the element intersection/facet corresponding corresponding to the dual grid scv // Get the element intersection/facet corresponding corresponding to the dual grid scv
...@@ -288,13 +279,11 @@ public: ...@@ -288,13 +279,11 @@ public:
// on that facet. // on that facet.
const auto& element = this->gridGeometry().element(scv.elementIndex()); const auto& element = this->gridGeometry().element(scv.elementIndex());
int idx = 0;
for (const auto& intersection : intersections(this->gridGeometry().gridView(), element)) for (const auto& intersection : intersections(this->gridGeometry().gridView(), element))
{ {
if (idx == scv.indexInElement()) if (intersection.indexInInside() == scv.indexInElement())
return velocityDirichlet_(intersection); return velocityDirichlet_(intersection);
++idx;
} }
DUNE_THROW(Dune::InvalidStateException, "No intersection found"); DUNE_THROW(Dune::InvalidStateException, "No intersection found");
} }
...@@ -339,7 +328,7 @@ private: ...@@ -339,7 +328,7 @@ private:
Scalar kinematicViscosity_; Scalar kinematicViscosity_;
Scalar time_ = 0.0; Scalar time_ = 0.0;
bool useNeumann_; bool useNeumann_;
bool useQuad_; bool interpolateExactVelocity_;
}; };
} // end namespace Dumux } // end namespace Dumux
......
...@@ -26,7 +26,6 @@ ...@@ -26,7 +26,6 @@
#define DUMUX_ANGELI_TEST_PROPERTIES_HH #define DUMUX_ANGELI_TEST_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh> #include <dune/grid/yaspgrid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dumux/freeflow/navierstokes/momentum/model.hh> #include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1p/model.hh> #include <dumux/freeflow/navierstokes/mass/1p/model.hh>
......
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