Commit c651312e authored by Ned Coltman's avatar Ned Coltman
Browse files

Merge branch 'feature/freeflow-periodic' into 'master'

Feature/freeflow periodic

See merge request !2827
parents db6de99e 57e25265
......@@ -63,7 +63,7 @@ public:
iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
}
decltype(auto) dereference() const
const Intersection& dereference() const
{
if constexpr (std::is_lvalue_reference_v<decltype(*std::declval<IntersectionIterator>())>)
return *iIt_;
......
......@@ -203,6 +203,9 @@ public:
auto& jacRow = (*jacobian_)[domainId];
auto& subRes = (*residual_)[domainId];
this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridGeometry, curSol[domainId]);
});
}
......@@ -551,6 +554,48 @@ private:
domainJ, gridGeometry(domainJ));
}
// build periodic contraints into the system matrix
template<std::size_t i, class JacRow, class Sol, class GG>
void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Sol& res, const GG& gridGeometry, const Sol& curSol)
{
if constexpr (GG::discMethod == DiscretizationMethod::box || GG::discMethod == DiscretizationMethod::fcstaggered)
{
for (const auto& m : gridGeometry.periodicVertexMap())
{
if (m.first < m.second)
{
auto& jac = jacRow[domainI];
// add the second row to the first
res[m.first] += res[m.second];
// enforce the solution of the first periodic DOF to the second one
res[m.second] = curSol[m.second] - curSol[m.first];
const auto end = jac[m.second].end();
for (auto it = jac[m.second].begin(); it != end; ++it)
jac[m.first][it.index()] += (*it);
// enforce constraint in second row
for (auto it = jac[m.second].begin(); it != end; ++it)
(*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
using namespace Dune::Hybrid;
forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
{
auto& jacCoupling = jacRow[couplingDomainId];
for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
jacCoupling[m.first][it.index()] += (*it);
for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
(*it) = 0.0;
});
}
}
}
}
//! pointer to the problem to be solved
ProblemTuple problemTuple_;
......
add_subdirectory(donea)
add_subdirectory(angeli)
add_subdirectory(kovasznay)
add_subdirectory(channel)
add_subdirectory(donea)
add_subdirectory(kovasznay)
add_subdirectory(periodic)
add_subdirectory(sincos)
dune_symlink_to_source_files(FILES convergence.sh)
dune_symlink_to_source_files(FILES "params.input" "periodic.dgf")
dumux_add_test(NAME test_ff_navierstokes_periodic
LABELS freeflow navierstokes
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMPILE_DEFINITIONS UPWINDSCHEMEORDER=1
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy --zeroThreshold {"velocity_liq \(m/s\)":1e-12}
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_periodic-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_periodic-00000.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_periodic params.input
-Problem.Name test_ff_navierstokes_periodic")
// -*- 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 3 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 NavierStokesTests
* \brief Stationary test for the staggered grid Navier-Stokes model with periodic BC
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/vtk/function.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/io/vtk/intersectionwriter.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include "properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
// define the type tag for this problem
using MomentumTypeTag = Properties::TTag::PeriodicTestMomentum;
using MassTypeTag = Properties::TTag::PeriodicTestMass;
// 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);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<MomentumTypeTag, Properties::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 MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
// the coupling manager
using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>();
// the problem (boundary conditions)
using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// the solution vector
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
x[momentumIdx].resize(momentumGridGeometry->numDofs());
x[massIdx].resize(massGridGeometry->numDofs());
// the grid variables
using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
// initialize coupling stencil first and then grid variables (need coupling variables)
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
massGridVariables->init(x[massIdx]);
momentumGridVariables->init(x[momentumIdx]);
// the assembler with time loop for instationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
std::make_tuple(momentumGridGeometry, massGridGeometry),
std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager);
// intialize the vtk output module
using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
ConformingIntersectionWriter faceVtk(momentumGridGeometry->gridView());
std::vector<std::size_t> dofIdx(x[momentumIdx].size());
for (const auto& facet : facets(momentumGridGeometry->gridView()))
{
const auto idx = momentumGridGeometry->gridView().indexSet().index(facet);
dofIdx[idx] = idx;
}
faceVtk.addField(dofIdx, "dofIdx");
using LinearSolver = Dumux::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// linearize & solve
Dune::Timer timer;
nonLinearSolver.solve(x);
std::vector<Dune::FieldVector<double,2>> faceVelocityVector(x[momentumIdx].size());
for (const auto& element : elements(momentumGridGeometry->gridView()))
{
auto fvGeometry = localView(*momentumGridGeometry);
fvGeometry.bind(element);
auto elemVolVars = localView(momentumGridVariables->curGridVolVars());
elemVolVars.bind(element, fvGeometry, x);
for (const auto& scv : scvs(fvGeometry))
faceVelocityVector[scv.dofIndex()][scv.dofAxis()] = elemVolVars[scv].velocity();
}
faceVtk.addField(faceVelocityVector, "velocityVector");
faceVtk.write("facedata", Dune::VTK::ascii);
// write vtk output
vtkWriter.write(1.0);
timer.stop();
const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
<< comm.size() << " processes.\n"
<< "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
[Grid]
File = periodic.dgf
[Problem]
Name = test_kovasznay # name passed to the output routines
EnableGravity = true
PrintL2Error = false
EnableInertiaTerms = true
PrintMatrix = true
UsePressureDifference = false
[Component]
LiquidDensity = 1
LiquidKinematicViscosity = 0.025
[ Newton ]
MaxSteps = 10
MaxRelativeShift = 1e-5
[Vtk]
WriteFaceData = false
[Flux]
UpwindWeight = 0.5
[Assembly]
NumericDifference.BaseEpsilon = 1e-4
DGF
INTERVAL
0 0 % lower left corner
1 1 % upper right corner
10 10 % number of cells in each direction
#
PERIODICFACETRANSFORMATION
1 0, 0 1 + 0 1 % make periodic in y-direction
#
BOUNDARYDOMAIN
DEFAULT 1
#
GRIDPARAMETER
OVERLAP 1
#
#
// -*- 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 3 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 NavierStokesTests
* \brief Stationary test for the staggered grid Navier-Stokes model with periodic BC
*/
#ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_PERIODIC_PROBLEM_HH
#define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_PERIODIC_PROBLEM_HH
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
namespace Dumux {
/*!
* \ingroup NavierStokesTests
* \brief Periodic test problem for the staggered grid
*
* A two-dimensional Navier-Stokes flow with a periodicity in one direction
* is considered.
*/
template <class TypeTag>
class PeriodicTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using NumEqVector = typename ParentType::NumEqVector;
using PrimaryVariables = typename ParentType::PrimaryVariables;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using Element = typename FVElementGeometry::Element;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
public:
PeriodicTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager)
{
usePressureDifference_ = getParam<bool>("Problem.UsePressureDifference", false);
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the temperature within the domain in [K].
*
* This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperature() const
{ return 298.0; }
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param globalPos The position of the center of the finite volume
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllDirichlet();
return values;
}
/*!
* \brief Returns Dirichlet boundary values at a given position.
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
{
return PrimaryVariables(0.0);
}
template<class ElementVolumeVariables>
NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{
NumEqVector source;
if constexpr (ParentType::isMomentumProblem())
{
if (usePressureDifference_ && scv.dofPosition()[1] < this->gridGeometry().bBoxMin()[1] + eps_)
{
const auto& frontalScvf = (*scvfs(fvGeometry, scv).begin());
source[Indices::momentumYBalanceIdx] = 100 * frontalScvf.area() / scv.volume();
}
}
return source;
}
// \}
/*!
* \name Volume terms
*/
// \{
//! Enable internal Dirichlet constraints
static constexpr bool enableInternalDirichletConstraints()
{ return !ParentType::isMomentumProblem(); }
/*!
* \brief Tag a degree of freedom to carry internal Dirichlet constraints.
* If true is returned for a dof, the equation for this dof is replaced
* by the constraint that its primary variable values must match the
* user-defined values obtained from the function internalDirichlet(),
* which must be defined in the problem.
*
* \param element The finite element
* \param scv The sub-control volume
*/
std::bitset<PrimaryVariables::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
std::bitset<PrimaryVariables::dimension> values;
for (const auto& intersection : intersections(this->gridGeometry().gridView(), element))
{
const auto center = intersection.geometry().center();
if (intersection.boundary() && intersection.neighbor() && center[1] > this->gridGeometry().bBoxMax()[1] - eps_)
values.set(0);
}
return values;
}
/*!
* \brief Define the values of internal Dirichlet constraints for a degree of freedom.
* \param element The finite element
* \param scv The sub-control volume
*/
PrimaryVariables internalDirichlet(const Element& element, const SubControlVolume& scv) const
{
return PrimaryVariables(1.0);
}
private:
static constexpr Scalar eps_ = 1e-6;
bool usePressureDifference_;
};
} // end namespace Dumux
#endif
// -*- 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 3 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 NavierStokesTests
* \brief Stationary test for the staggered grid Navier-Stokes model with periodic BC
*/
#ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_PERIODIC_PROPERTIES_HH
#define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_PERIODIC_PROPERTIES_HH
#include <dune/grid/spgrid.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include <dumux/discretization/fcstaggered.hh>
#include <dumux/discretization/cctpfa.hh>
#include "problem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
struct PeriodicTest {};
struct PeriodicTestMomentum { using InheritsFrom = std::tuple<PeriodicTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
struct PeriodicTestMass { using InheritsFrom = std::tuple<PeriodicTest, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::PeriodicTest>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::PeriodicTest> { using type = Dune::SPGrid<double, 2>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::PeriodicTest> { using type = PeriodicTestProblem<TypeTag> ; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::PeriodicTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::PeriodicTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::PeriodicTest> { static constexpr bool value = true; };
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::PeriodicTest>
{
using Traits = MultiDomainTraits<TTag::PeriodicTestMomentum, TTag::PeriodicTestMass>;
using type = StaggeredFreeFlowCouplingManager<Traits>;
};