Commit 6f9f5faf authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[md][stokesdarcy] Add convergence test with analytical solution

parent 13af9d39
add_subdirectory(convergencetest)
add_input_file_links()
add_executable(test_md_boundary_darcy1p_stokes1p EXCLUDE_FROM_ALL main.cc)
......
add_input_file_links()
dune_symlink_to_source_files(FILES "convergencetest.py")
dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_convtest
SOURCES main.cc
LABELS multidomain multidomain_boundary stokesdarcy
TIMEOUT 1000
CMAKE_GUARD HAVE_UMFPACK
COMMAND ./convergencetest.py
CMD_ARGS test_md_boundary_darcy1p_stokes1p_convtest params.input)
#!/usr/bin/env python
from math import *
import subprocess
import sys
if len(sys.argv) < 2:
sys.stderr.write('Please provide a single argument <testname> to the script\n')
sys.exit(1)
testname = str(sys.argv[1])
testargs = [str(i) for i in sys.argv][2:]
# remove the old log files
subprocess.call(['rm', testname + '_stokes.log'])
print("Removed old log file ({})!".format(testname + '_stokes.log'))
subprocess.call(['rm', testname + '_darcy.log'])
print("Removed old log file ({})!".format(testname + '_darcy.log'))
# do the runs with different refinement
for i in [0, 1, 2]:
subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
'-Vtk.OutputName', testname])
def checkRatesStokes():
# check the rates and append them to the log file
logfile = open(testname + '_stokes.log', "r+")
errorP = []
errorVx = []
errorVy = []
for line in logfile:
line = line.strip("\n")
line = line.strip("\[ConvergenceTest\]")
line = line.split()
errorP.append(float(line[2]))
errorVx.append(float(line[5]))
errorVy.append(float(line[8]))
resultsP = []
resultsVx = []
resultsVy = []
logfile.truncate(0)
logfile.write("n\terrorP\t\trateP\t\terrorVx\t\trateVx\t\terrorVy\t\trateVy\n")
logfile.write("-"*50 + "\n")
for i in range(len(errorP)-1):
if isnan(errorP[i]) or isinf(errorP[i]):
continue
if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12) and (errorVx[i] < 1e-12 or errorVx[i+1] < 1e-12) and (errorVy[i] < 1e-12 or errorVy[i+1] < 1e-12)):
rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
rateVx = (log(errorVx[i])-log(errorVx[i+1]))/log(2)
rateVy = (log(errorVy[i])-log(errorVy[i+1]))/log(2)
message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP, errorVx[i], rateVx, errorVy[i], rateVy)
logfile.write(message)
resultsP.append(rateP)
resultsVx.append(rateVx)
resultsVy.append(rateVy)
else:
logfile.write("error: exact solution!?")
i = len(errorP)-1
message = "{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\n".format(i, errorP[i], "", errorVx[i], "", errorVy[i], "")
logfile.write(message)
logfile.close()
print("\nComputed the following convergence rates for {}:\n".format(testname))
subprocess.call(['cat', testname + '_stokes.log'])
return {"p" : resultsP, "v_x" : resultsVx, "v_y" : resultsVy}
def checkRatesDarcy():
# check the rates and append them to the log file
logfile = open(testname + '_darcy.log', "r+")
errorP = []
for line in logfile:
line = line.strip("\n")
line = line.strip("\[ConvergenceTest\]")
line = line.split()
errorP.append(float(line[2]))
resultsP = []
logfile.truncate(0)
logfile.write("n\terrorP\t\trateP\n")
logfile.write("-"*50 + "\n")
for i in range(len(errorP)-1):
if isnan(errorP[i]) or isinf(errorP[i]):
continue
if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)):
rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP)
logfile.write(message)
resultsP.append(rateP)
else:
logfile.write("error: exact solution!?")
i = len(errorP)-1
message = "{}\t{:0.4e}\n".format(i, errorP[i], "")
logfile.write(message)
logfile.close()
print("\nComputed the following convergence rates for {}:\n".format(testname))
subprocess.call(['cat', testname + '_darcy.log'])
return {"p" : resultsP}
def checkRatesStokesAndDarcy():
resultsStokes = checkRatesStokes()
resultsDarcy = checkRatesDarcy()
def mean(numbers):
return float(sum(numbers)) / len(numbers)
# check the rates, we expect rates around 2
if mean(resultsStokes["p"]) < 2.05 and mean(resultsStokes["p"]) < 1.84:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
if mean(resultsStokes["v_x"]) < 2.05 and mean(resultsStokes["v_x"]) < 1.95:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for x-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
if mean(resultsStokes["v_y"]) < 2.05 and mean(resultsStokes["v_y"]) < 1.95:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for y-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
if mean(resultsDarcy["p"]) < 2.05 and mean(resultsDarcy["p"]) < 1.95:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
checkRatesStokesAndDarcy()
// -*- 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 BoundaryTests
* \brief A test problem for the coupled Stokes/Darcy problem (1p).
* The analytical solution is given in Shiue et al., 2018:
* "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/partial.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <test/freeflow/navierstokes/l2error.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "problem_darcy.hh"
#include "problem_stokes.hh"
namespace Dumux {
namespace Properties {
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::StokesOneP>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::DarcyOneP>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
} // end namespace Properties
} // end namespace Dumux
/*!
* \brief Creates analytical solution.
* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Scalar, class Problem>
auto createStokesAnalyticalSolution(const Problem& problem)
{
const auto& gridGeometry = problem.gridGeometry();
using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
static constexpr auto dim = GridView::dimension;
static constexpr auto dimWorld = GridView::dimensionworld;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
std::vector<Scalar> analyticalPressure;
std::vector<VelocityVector> analyticalVelocity;
std::vector<Scalar> analyticalVelocityOnFace;
analyticalPressure.resize(gridGeometry.numCellCenterDofs());
analyticalVelocity.resize(gridGeometry.numCellCenterDofs());
analyticalVelocityOnFace.resize(gridGeometry.numFaceDofs());
using Indices = typename Problem::Indices;
for (const auto& element : elements(gridGeometry.gridView()))
{
auto fvGeometry = localView(gridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto ccDofIdx = scv.dofIndex();
auto ccDofPosition = scv.dofPosition();
auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
// 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 = problem.analyticalSolution(faceDofPosition);
analyticalVelocityOnFace[faceDofIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace);
}
/*!
* \brief Creates analytical solution.
* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Scalar, class Problem>
auto createDarcyAnalyticalSolution(const Problem& problem)
{
const auto& gridGeometry = problem.gridGeometry();
using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
static constexpr auto dim = GridView::dimension;
static constexpr auto dimWorld = GridView::dimensionworld;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
std::vector<Scalar> analyticalPressure;
std::vector<VelocityVector> analyticalVelocity;
analyticalPressure.resize(gridGeometry.numDofs());
analyticalVelocity.resize(gridGeometry.numDofs());
for (const auto& element : elements(gridGeometry.gridView()))
{
auto fvGeometry = localView(gridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
const auto ccDofIdx = scv.dofIndex();
const auto ccDofPosition = scv.dofPosition();
const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim];
for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx];
}
}
return std::make_tuple(analyticalPressure, analyticalVelocity);
}
template<class Problem, class SolutionVector>
void printStokesL2Error(const Problem& problem, const SolutionVector& x)
{
using namespace Dumux;
using Scalar = double;
using TypeTag = Properties::TTag::StokesOneP;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
const auto l2error = L2Error::calculateL2Error(problem, x);
const int numCellCenterDofs = problem.gridGeometry().numCellCenterDofs();
const int numFaceDofs = problem.gridGeometry().numFaceDofs();
std::ostream tmpOutputObject(std::cout.rdbuf()); // create temporary output with fixed formatting without affecting std::cout
tmpOutputObject << 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;
// write the norm into a log file
std::ofstream logFile;
logFile.open(problem.name() + ".log", std::ios::app);
logFile << "[ConvergenceTest] L2(p) = " << l2error.first[Indices::pressureIdx] << " L2(vx) = " << l2error.first[Indices::velocityXIdx] << " L2(vy) = " << l2error.first[Indices::velocityYIdx] << std::endl;
logFile.close();
}
template<class Problem, class SolutionVector>
void printDarcyL2Error(const Problem& problem, const SolutionVector& x)
{
using namespace Dumux;
using Scalar = double;
Scalar l2error = 0.0;
for (const auto& element : elements(problem.gridGeometry().gridView()))
{
auto fvGeometry = localView(problem.gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
const auto dofIdx = scv.dofIndex();
const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/];
l2error += scv.volume()*(delta*delta);
}
}
using std::sqrt;
l2error = sqrt(l2error);
const auto numDofs = problem.gridGeometry().numDofs();
std::ostream tmp(std::cout.rdbuf());
tmp << std::setprecision(8) << "** L2 error (abs) for "
<< std::setw(6) << numDofs << " cc dofs "
<< std::scientific
<< "L2 error = " << l2error
<< std::endl;
// write the norm into a log file
std::ofstream logFile;
logFile.open(problem.name() + ".log", std::ios::app);
logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl;
logFile.close();
}
int main(int argc, char** argv) try
{
using namespace Dumux;
// 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);
// Define the sub problem type tags
using StokesTypeTag = Properties::TTag::StokesOneP;
using DarcyTypeTag = Properties::TTag::DarcyOneP;
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
// we compute on the leaf grid view
const auto& darcyGridView = darcyGridManager.grid().leafGridView();
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
auto stokesGridGeometry = std::make_shared<StokesGridGeometry>(stokesGridView);
stokesGridGeometry->update();
using DarcyGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
auto darcyGridGeometry = std::make_shared<DarcyGridGeometry>(darcyGridView);
darcyGridGeometry->update();
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesGridGeometry, darcyGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesGridGeometry, couplingManager);
using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager);
// the solution vector
Traits::SolutionVector sol;
sol[stokesCellCenterIdx].resize(stokesGridGeometry->numCellCenterDofs());
sol[stokesFaceIdx].resize(stokesGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyGridGeometry->numDofs());
// get a solution vector storing references to the two Stokes solution vectors
auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
couplingManager->init(stokesProblem, darcyProblem, sol);
// the grid variables
using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesGridGeometry);
stokesGridVariables->init(stokesSol);
using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
using Scalar = typename Traits::Scalar;
StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
const auto stokesAnalyticalSolution = createStokesAnalyticalSolution<Scalar>(*stokesProblem);
stokesVtkWriter.addField(std::get<0>(stokesAnalyticalSolution), "pressureExact");
stokesVtkWriter.addField(std::get<1>(stokesAnalyticalSolution), "velocityExact");
stokesVtkWriter.addFaceField(std::get<2>(stokesAnalyticalSolution), "faceVelocityExact");
stokesVtkWriter.write(0.0);
VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
const auto darcyAnalyticalSolution = createDarcyAnalyticalSolution<Scalar>(*darcyProblem);
darcyVtkWriter.addField(std::get<0>(darcyAnalyticalSolution), "pressureExact");
darcyVtkWriter.addField(std::get<1>(darcyAnalyticalSolution), "velocityExact");
darcyVtkWriter.write(0.0);
// the assembler for a stationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
std::make_tuple(stokesGridGeometry->faceFVGridGeometryPtr(),
stokesGridGeometry->cellCenterFVGridGeometryPtr(),
darcyGridGeometry),
std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
stokesGridVariables->cellCenterGridVariablesPtr(),
darcyGridVariables),
couplingManager);
// the linear solver
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// solve the non-linear system
nonLinearSolver.solve(sol);
// write vtk output
stokesVtkWriter.write(1.0);
darcyVtkWriter.write(1.0);
printStokesL2Error(*stokesProblem, stokesSol);
printDarcyL2Error(*darcyProblem, sol[darcyIdx]);
////////////////////////////////////////////////////////////
// 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;
}
import sympy as sp
import numpy as np
from sympy import *
# divergence of a matrix
def div(M):
return np.array([sp.diff(M[0,0],x) + sp.diff(M[0,1],y) , sp.diff(M[1,0],x) + sp.diff(M[1,1],y) ])