Commit db869d96 authored by Timo Koch's avatar Timo Koch Committed by Timo Koch
Browse files

[test][ff] Improve output of Donea momentum test

parent 89d1f879
......@@ -46,6 +46,42 @@
#include "properties_momentum.hh"
template<class GridGeometry, class GridVariables, class SolutionVector>
void updateVelocities(
std::vector<Dune::FieldVector<double, 2>>& velocity,
std::vector<Dune::FieldVector<double, 2>>& faceVelocity,
const GridGeometry& gridGeometry,
const GridVariables& gridVariables,
const SolutionVector& x
){
auto fvGeometry = localView(gridGeometry);
auto elemVolVars = localView(gridVariables.curGridVolVars());
for (const auto& element : elements(gridGeometry.gridView()))
{
fvGeometry.bind(element);
elemVolVars.bind(element, fvGeometry, x);
for (const auto& scv : scvs(fvGeometry))
{
const auto& vars = elemVolVars[scv];
velocity[scv.elementIndex()][scv.dofAxis()] += 0.5*vars.velocity();
faceVelocity[scv.dofIndex()][scv.dofAxis()] = vars.velocity();
}
}
}
template<class GridGeometry>
void updateRank(
std::vector<int>& rank,
const GridGeometry& gridGeometry
){
for (const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
{
const auto eIdxGlobal = gridGeometry.elementMapper().index(element);
rank[eIdxGlobal] = gridGeometry.gridView().comm().rank();
}
}
int main(int argc, char** argv)
{
using namespace Dumux;
......@@ -82,75 +118,56 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>>;
const auto dofMapper = LinearSolverTraits<GridGeometry>::DofMapper(gridGeometry->gridView());
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), dofMapper);
std::vector<Dune::FieldVector<double, 2>> velocity(gridGeometry->gridView().size(0));
std::vector<Dune::FieldVector<double, 2>> faceVelocity(x.size());
updateVelocities(velocity, faceVelocity, *gridGeometry, *gridVariables, x);
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
nonLinearSolver.solve(x);
////////////////////////////////////////////////////////////
// write VTK output
////////////////////////////////////////////////////////////
std::vector<int> rank(gridGeometry->gridView().size(0));
updateRank(rank, *gridGeometry);
std::vector<Dune::FieldVector<double,2>> velocity(gridGeometry->gridView().size(0));
std::vector<Dune::FieldVector<double,2>> faceVelocityVector(x.size());
for (const auto& element : elements(gridGeometry->gridView()))
std::vector<std::size_t> dofIdx(x.size());
for (const auto& facet : facets(gridGeometry->gridView()))
{
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bind(element);
auto elemVolVars = localView(gridVariables->curGridVolVars());
elemVolVars.bind(element, fvGeometry, x);
const auto eIdx = gridGeometry->elementMapper().index(element);
for (const auto& scv : scvs(fvGeometry))
{
velocity[eIdx][scv.dofAxis()] += 0.5*elemVolVars[scv].velocity();
faceVelocityVector[scv.dofIndex()][scv.dofAxis()] = elemVolVars[scv].velocity();
}
const auto idx = gridGeometry->gridView().indexSet().index(facet);
dofIdx[idx] = idx;
}
Dune::VTKWriter<std::decay_t<decltype(gridGeometry->gridView())>> writer(gridGeometry->gridView());
using Field = Vtk::template Field<std::decay_t<decltype(gridGeometry->gridView())>>;
Dune::VTKWriter<typename GridGeometry::GridView> writer(gridGeometry->gridView());
using Field = Vtk::template Field<typename GridGeometry::GridView>;
writer.addCellData(Field(
gridGeometry->gridView(), gridGeometry->elementMapper(), velocity,
"velocity", /*numComp*/2, /*codim*/0
).get());
std::vector<int> rank(gridGeometry->gridView().size(0));
for (const auto& element : elements(gridGeometry->gridView(), Dune::Partitions::interior))
{
const auto eIdxGlobal = gridGeometry->elementMapper().index(element);
rank[eIdxGlobal] = gridGeometry->gridView().comm().rank();
}
writer.addCellData(rank, "rank");
writer.write("donea_new_momentum");
writer.write("donea_momentum_0");
ConformingIntersectionWriter faceVtk(gridGeometry->gridView());
std::vector<std::size_t> dofIdx(x.size());
for (const auto& facet : facets(gridGeometry->gridView()))
{
const auto idx = gridGeometry->gridView().indexSet().index(facet);
dofIdx[idx] = idx;
}
faceVtk.addField(dofIdx, "dofIdx");
faceVtk.addField(faceVelocityVector, "velocityVector");
const auto partionType = [&](const auto& is, const auto idx)
{
faceVtk.addField(faceVelocity, "velocityVector");
faceVtk.addField([&](const auto& is, const auto idx) {
const auto& facet = is.inside().template subEntity <1> (is.indexInInside());
return facet.partitionType();
};
}, "partitionType");
faceVtk.write("donea_momentum_face_0_" + std::to_string(gridGeometry->gridView().comm().rank()), Dune::VTK::ascii);
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>>;
const auto dofMapper = LinearSolverTraits<GridGeometry>::DofMapper(gridGeometry->gridView());
auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), dofMapper);
faceVtk.addField(partionType, "partitionType");
faceVtk.write("facedata_" + std::to_string(gridGeometry->gridView().comm().rank()), Dune::VTK::ascii);
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
nonLinearSolver.solve(x);
////////////////////////////////////////////////////////////
// write VTK output
////////////////////////////////////////////////////////////
updateVelocities(velocity, faceVelocity, *gridGeometry, *gridVariables, x);
writer.write("donea_momentum_1");
faceVtk.write("donea_momentum_face_1_" + std::to_string(gridGeometry->gridView().comm().rank()), Dune::VTK::ascii);
////////////////////////////////////////////////////////////
// finalize, print parameters and Dumux message to say goodbye
......
......@@ -23,7 +23,7 @@ MaxRelativeShift = 1e-5
WriteFaceData = false
[LinearSolver]
Type = bicgstabsolver
Type = cgsolver
#Type = restartedgmressolver
Verbosity = 1
ResidualReduction = 1e-13
......
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