diff --git a/changelog-entries/678.md b/changelog-entries/678.md new file mode 100644 index 000000000..7b6077542 --- /dev/null +++ b/changelog-entries/678.md @@ -0,0 +1 @@ +- Added new [free-flow-over-porous-media tutorial](https://precice.org/tutorials-free-flow-over-porous-media-2d.html) from [the example case in dumux-adapter](https://github.com/precice/dumux-adapter/tree/77e0fe5ca0dc6a1414d6cce5813ca914f0904259/examples/ff-pm) and updated to follow tutorials structure. diff --git a/free-flow-over-porous-media/README.md b/free-flow-over-porous-media/README.md new file mode 100644 index 000000000..2a08750c4 --- /dev/null +++ b/free-flow-over-porous-media/README.md @@ -0,0 +1,71 @@ +--- +title: Free flow over porous media 2D +permalink: tutorials-free-flow-over-porous-media-2d.html +keywords: DuMux, porous media +summary: Flow-flow coupling example with porous media field and free flow field. +--- + +{% note %} +Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/free-flow-over-porous-media-2d). Read how in the [tutorials introduction](https://precice.org/tutorials.html). +{% endnote %} + +## Setup + +This tutorial solves a coupled system consisting of a one-phase free flow and a one-phase flow in a porous media. + +A pressure gradient is applied to the free flow domain from left to right. The top edge of the free-flow is a non-permeable wall with no-slip boundary conditions. In the porous media, there is a no-flow condition across the domain boundaries (left, bottom, and right boundaries). At the interface, a no-slip condition applies. The case is stationary (solved to a steady-state solution). + +The setting is illustrated in the following figure: + +![Free flow over porous media setup](images/tutorials-free-flow-over-porous-media-setup.png) + +## Configuration + +preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)): + +![preCICE configuration visualization](images/tutorials-free-flow-over-porous-media-precice-config-visualization.png) + +## Available solvers + +Both the participants are computed using the simulation code [DuMux](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/). + +## Solver setup + +To solve the flows with the DuMux framework, the necessary DUNE modules need to be downloaded and set up. This is done by running `sh setup-dumux.sh` in the tutorial folder. + +Note that if an existing installation of DUNE modules is detected in a default location, this may lead to problems in running the `setup-dumux.sh` script. The script suppresses the environment variable `DUNE_CONTROL_PATH`. + +To only recompile the participants, run `sh compile-dumux-cases.sh` in the tutorial folder. + +## Running the simulation + +Each participant has a `run.sh` script. + +To run the free-flow participant, run: + +```bash +cd free-flow-dumux +./run.sh +``` + +To run the porous-media participant, run: + +```bash +cd porous-media-dumux +./run.sh +``` + +Participants can be executed only in serial. Parallel execution is not supported. The case takes approximately two minutes to finish. + +## Post-processing + +Both participants write VTU outputs, which can be viewed using ParaView. + +## Further information + +The results of the pressure and the velocity fields are as follows: + +![Free flow over porous media results - pressure](images/tutorials-free-flow-over-porous-media-result-pressure.png) +![Free flow over porous media results - velocity](images/tutorials-free-flow-over-porous-media-result-ux.png) + +Each solver folder contains an input file (`params.input`) that will be passed to the solver executables. This is a DuMUX input file describing the simulation setting, e.g., pressure, mesh size, time stepping, etc. diff --git a/free-flow-over-porous-media/clean-tutorial.sh b/free-flow-over-porous-media/clean-tutorial.sh new file mode 120000 index 000000000..4713f5092 --- /dev/null +++ b/free-flow-over-porous-media/clean-tutorial.sh @@ -0,0 +1 @@ +../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/free-flow-over-porous-media/compile-dumux-cases.sh b/free-flow-over-porous-media/compile-dumux-cases.sh new file mode 100755 index 000000000..c43bd618c --- /dev/null +++ b/free-flow-over-porous-media/compile-dumux-cases.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env sh +set -e -u + +cd free-flow-dumux/build-cmake/appl +make free_flow_dumux +cd ../../../porous-media-dumux/build-cmake/appl +make porous_media_dumux +cd ../../../ + +# Move free-flow-dumux and porous-media-dumux executables to the participant folder level +mv free-flow-dumux/build-cmake/appl/free_flow_dumux free-flow-dumux/ +mv porous-media-dumux/build-cmake/appl/porous_media_dumux porous-media-dumux/ diff --git a/free-flow-over-porous-media/free-flow-dumux/CMakeLists.txt b/free-flow-over-porous-media/free-flow-dumux/CMakeLists.txt new file mode 100644 index 000000000..7d3d9b426 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/CMakeLists.txt @@ -0,0 +1,26 @@ +cmake_minimum_required(VERSION 3.13) +project(free-flow-dumux CXX) + +if(NOT (dune-common_DIR OR dune-common_ROOT OR + "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*")) + string(REPLACE ${PROJECT_NAME} dune-common dune-common_DIR + ${PROJECT_BINARY_DIR}) +endif() + +#find dune-common and set the module path +find_package(dune-common REQUIRED) +list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules" + ${dune-common_MODULE_PATH}) + +#include the dune macros +include(DuneMacros) + +# start a dune project with information from dune.module +dune_project() + +dune_enable_all_packages() + +add_subdirectory(appl) + +# finalize the dune project, e.g. generating config.h etc. +finalize_dune_project(GENERATE_CONFIG_H_CMAKE) diff --git a/free-flow-over-porous-media/free-flow-dumux/appl/CMakeLists.txt b/free-flow-over-porous-media/free-flow-dumux/appl/CMakeLists.txt new file mode 100644 index 000000000..e9c64417a --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/appl/CMakeLists.txt @@ -0,0 +1,5 @@ +add_executable(free_flow_dumux main.cc) +target_compile_definitions(free_flow_dumux PUBLIC "ENABLEMONOLITHIC=0") +target_link_libraries(free_flow_dumux PRIVATE dumux-precice) + +add_input_file_links() \ No newline at end of file diff --git a/free-flow-over-porous-media/free-flow-dumux/appl/main.cc b/free-flow-over-porous-media/free-flow-dumux/appl/main.cc new file mode 100644 index 000000000..12a2e4648 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/appl/main.cc @@ -0,0 +1,515 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief A test problem for the coupled Stokes/Darcy problem (1p) + */ +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "properties.hh" + +#include "dumux-precice/couplingadapter.hh" + +#include +#include +#include + +// TODO +// Helper function to put pressure on interface + +template +auto velocityAtInterface(const ElementVolumeVariables &elemVolVars, + const SubControlVolumeFace &scvf) +{ + assert(scvf.isFrontal()); + const double scalarVelocity = elemVolVars[scvf.insideScvIdx()].velocity(); + auto velocity = scvf.unitOuterNormal(); + velocity[scvf.normalAxis()] = scalarVelocity; + return velocity; +} + +template +auto pressureAtInterface(const Problem *problem, + const Element &element, + const SubControlVolumeFace &scvf, + const FVElementGeometry &fvGeometry, + const ElementVolumeVariables &elemVolVars, + const ElementFluxVariablesCache &elemFluxVarsCache) +{ + using ElementBoundaryTypes = + Dumux::GetPropType; + using LocalResidual = + Dumux::GetPropType; + using NumEqVector = Dumux::NumEqVector< + Dumux::GetPropType>; +#if DUMUX_VERSION_MAJOR >= 3 && DUMUX_VERSION_MINOR >= 9 + using SlipVelocityPolicy = Dumux::NavierStokesSlipVelocity< + typename FVElementGeometry::GridGeometry::DiscretizationMethod, + Dumux::NavierStokes::SlipConditions::BJ>; + using FluxHelper = Dumux::NavierStokesMomentumBoundaryFlux< + typename FVElementGeometry::GridGeometry::DiscretizationMethod, + SlipVelocityPolicy>; +#else + using FluxHelper = Dumux::NavierStokesMomentumBoundaryFluxHelper; +#endif + + ElementBoundaryTypes elemBcTypes; + auto localResidual = LocalResidual(problem); + elemBcTypes.update(*problem, element, fvGeometry); + + NumEqVector flux(0.0); + NumEqVector neumannFlux(0.0); + const auto &scv = fvGeometry.scv(scvf.insideScvIdx()); + for (const auto &otherScvf : scvfs(fvGeometry, scv)) { + if (otherScvf.index() == scvf.index()) + continue; + flux += localResidual.maybeHandleNeumannBoundary( + *problem, element, fvGeometry, elemVolVars, elemBcTypes, + elemFluxVarsCache, otherScvf); + } + flux += FluxHelper::fixedPressureMomentumFlux( + *problem, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, 0.0, + /*zeroNormalVelocityGradient=*/false)[scvf.normalAxis()] * + scvf.area(); + return -1 * scvf.directionSign() * flux / scvf.area(); +} + +template +void setInterfacePressures(const std::shared_ptr problem, + const GridVariables &gridVars, + const SolutionVector &sol, + const std::string meshName, + const std::string dataName) +{ + const auto &gridGeometry = problem->gridGeometry(); + auto fvGeometry = localView(gridGeometry); + auto elemVolVars = localView(gridVars.curGridVolVars()); + auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache()); + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + + for (const auto &element : elements(gridGeometry.gridView())) { + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, sol); + elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (couplingParticipant.isCoupledEntity(scvf.index())) { + // TODO: What to do here? + const auto p = pressureAtInterface( + problem.get(), element, scvf, fvGeometry, elemVolVars, + elemFluxVarsCache); + couplingParticipant.writeScalarQuantityOnFace( + meshName, dataName, scvf.index(), p); + } + } + } +} + +template +void setInterfaceVelocities(const MomentumProblem &problem, + const GridVariables &gridVars, + const SolutionVector &sol, + const std::string meshName, + const std::string dataName) +{ + const auto &gridGeometry = problem.gridGeometry(); + auto fvGeometry = localView(gridGeometry); + auto elemVolVars = localView(gridVars.curGridVolVars()); + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + + for (const auto &element : elements(gridGeometry.gridView())) { + fvGeometry.bindElement(element); + elemVolVars.bindElement(element, fvGeometry, sol); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (couplingParticipant.isCoupledEntity(scvf.index())) { + // TODO: What to do here? + const auto v = + velocityAtInterface(elemVolVars, scvf)[scvf.normalAxis()]; + couplingParticipant.writeScalarQuantityOnFace( + meshName, dataName, scvf.index(), v); + } + } + } +} + +template +std::tuple writeVelocitiesOnInterfaceToFile( + const std::string &meshName, + const std::string &filename, + const std::shared_ptr problem, + const GridVariables &gridVars, + const SolutionVector &sol) +{ + const auto &gridGeometry = problem->gridGeometry(); + auto fvGeometry = localView(gridGeometry); + auto elemVolVars = localView(gridVars.curGridVolVars()); + + const auto &couplingParticipant = + Dumux::Precice::CouplingAdapter::getInstance(); + + std::ofstream ofs(filename + ".csv", + std::ofstream::out | std::ofstream::trunc); + ofs << "x,y,"; + if (couplingParticipant.getMeshDimensions(meshName) == 3) + ofs << "z,"; + ofs << "velocityY" + << "\n"; + + double min = std::numeric_limits::max(); + double max = std::numeric_limits::min(); + double sum = 0.; + for (const auto &element : elements(gridGeometry.gridView())) { + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, sol); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (couplingParticipant.isCoupledEntity(scvf.index())) { + const auto &pos = scvf.center(); + for (int i = 0; + i < couplingParticipant.getMeshDimensions(meshName); ++i) { + ofs << pos[i] << ","; + } + const double v = + velocityAtInterface(elemVolVars, scvf)[scvf.normalAxis()]; + max = std::max(v, max); + min = std::min(v, min); + sum += v; + const int prec = ofs.precision(); + ofs << std::setprecision(std::numeric_limits::digits10 + + 1) + << v << "\n"; + ofs.precision(prec); + } + } + } + + ofs.close(); + + return std::make_tuple(min, max, sum); +} + +template +void writePressuresOnInterfaceToFile(const std::string &meshName, + const std::string &filename, + const std::shared_ptr problem, + const GridVariables &gridVars, + const SolutionVector &sol) +{ + const auto &gridGeometry = problem->gridGeometry(); + auto fvGeometry = localView(gridGeometry); + auto elemVolVars = localView(gridVars.curGridVolVars()); + auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache()); + + const auto &couplingParticipant = + Dumux::Precice::CouplingAdapter::getInstance(); + + std::ofstream ofs(filename + ".csv", + std::ofstream::out | std::ofstream::trunc); + ofs << "x,y,"; + if (couplingParticipant.getMeshDimensions(meshName) == 3) + ofs << "z,"; + ofs << "pressure" + << "\n"; + for (const auto &element : elements(gridGeometry.gridView())) { + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, sol); + elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (couplingParticipant.isCoupledEntity(scvf.index())) { + const auto &pos = scvf.center(); + for (int i = 0; + i < couplingParticipant.getMeshDimensions(meshName); ++i) { + ofs << pos[i] << ","; + } + const double p = pressureAtInterface( + problem.get(), element, scvf, fvGeometry, elemVolVars, + elemFluxVarsCache); + ofs << p << "\n"; + } + } + } + + ofs.close(); +} + +int main(int argc, char **argv) +{ + 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 MomentumTypeTag = Properties::TTag::FreeFlowSubMomentum; + using MassTypeTag = Properties::TTag::FreeFlowSubMass; + + // try to create a grid (from the given grid file or the input file) + using FreeFlowGridManager = + Dumux::GridManager>; + FreeFlowGridManager freeFlowGridManager; + freeFlowGridManager.init("FreeFlow"); // pass parameter group + + // we compute on the leaf grid view + const auto &freeFlowGridView = freeFlowGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using MomentumGridGeometry = + GetPropType; + auto momentumGridGeometry = + std::make_shared(freeFlowGridView); + using MassGridGeometry = GetPropType; + auto massGridGeometry = + std::make_shared(freeFlowGridView); + + // create the coupling manager to couple the two subproblems of the freeflow participant + using CouplingManager = + GetPropType; + auto couplingManager = std::make_shared(); + constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex; + constexpr auto massIdx = CouplingManager::freeFlowMassIndex; + + // the problem (initial and boundary conditions) + using MomentumProblem = GetPropType; + auto momentumProblem = std::make_shared( + momentumGridGeometry, couplingManager); + using MassProblem = GetPropType; + auto massProblem = + std::make_shared(massGridGeometry, couplingManager); + + // the solution vector + using Traits = MultiDomainTraits; + using SolutionVector = Traits::SolutionVector; + SolutionVector sol; + + // Initialize preCICE.Tell preCICE about: + std::string preciceConfigFilename = "../precice-config.xml"; + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + couplingParticipant.announceSolver("FreeFlow", preciceConfigFilename, + mpiHelper.rank(), mpiHelper.size()); + + const std::string meshName("Free-Flow-Mesh"); // mesh name + const int dim = couplingParticipant.getMeshDimensions(meshName); + std::cout << dim << " " << int(MassGridGeometry::GridView::dimension) + << std::endl; + if (dim != int(MassGridGeometry::GridView::dimension)) + DUNE_THROW(Dune::InvalidStateException, "Dimensions do not match"); + + // GET mesh corodinates + const double xMin = + getParamFromGroup>("Darcy", "Grid.LowerLeft")[0]; + const double xMax = + getParamFromGroup>("Darcy", "Grid.UpperRight")[0]; + std::vector coords; //( dim * vertexSize ); + std::vector coupledScvfIndices; + + for (const auto &element : elements(freeFlowGridView)) { + auto fvGeometry = localView(*momentumGridGeometry); + fvGeometry.bindElement(element); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (!scvf.isFrontal()) + continue; + static constexpr auto eps = 1e-7; + const auto &pos = scvf.center(); + if (pos[1] < momentumGridGeometry->bBoxMin()[1] + eps) { + if (pos[0] > xMin - eps && pos[0] < xMax + eps) { + coupledScvfIndices.push_back(scvf.index()); + for (const auto p : pos) + coords.push_back(p); + } + } + } + } + + couplingParticipant.setMesh(meshName, coords); + couplingParticipant.createIndexMapping(coupledScvfIndices); + + const std::string dataNameV("Velocity"); + const std::string dataNameP("Pressure"); + couplingParticipant.announceQuantity(meshName, dataNameV); + couplingParticipant.announceQuantity(meshName, dataNameP); + + // apply initial solution for instationary problems + momentumProblem->applyInitialSolution(sol[momentumIdx]); + massProblem->applyInitialSolution(sol[massIdx]); + + // the grid variables + using MomentumGridVariables = + GetPropType; + auto momentumGridVariables = std::make_shared( + momentumProblem, momentumGridGeometry); + using MassGridVariables = + GetPropType; + auto massGridVariables = + std::make_shared(massProblem, massGridGeometry); + + // initialize the coupling manager and the grid variables of the subproblems + couplingManager->init( + momentumProblem, massProblem, + std::make_tuple(momentumGridVariables, massGridVariables), sol); + momentumGridVariables->init(sol[momentumIdx]); + massGridVariables->init(sol[massIdx]); + + bool writeInterfaceDataToFile = + getParamFromGroup("Output", "EnableCSVWriter", false); + + // intialize the vtk output module + using IOFields = GetPropType; + VtkOutputModule freeFlowVtkWriter(*massGridVariables, sol[massIdx], + massProblem->name()); + IOFields::initOutputModule(freeFlowVtkWriter); + freeFlowVtkWriter.addVelocityOutput( + std::make_shared>()); + freeFlowVtkWriter.addField(massProblem->getAnalyticalVelocityX(), + "analyticalV_x"); + freeFlowVtkWriter.write(0.0); + + if (couplingParticipant.requiresToWriteInitialData()) { + setInterfacePressures( + momentumProblem, *momentumGridVariables, sol[momentumIdx], meshName, + dataNameP); + couplingParticipant.writeQuantityToOtherSolver(meshName, dataNameP); + } + couplingParticipant.initialize(); + couplingParticipant.initializeCheckpoint(sol[momentumIdx], + *momentumGridVariables); + couplingParticipant.initializeCheckpoint(sol[massIdx], *massGridVariables); + + // the assembler for a stationary problem + using Assembler = + MultiDomainFVAssembler; + auto assembler = std::make_shared( + std::make_tuple(momentumProblem, massProblem), + std::make_tuple(momentumGridGeometry, massGridGeometry), + std::make_tuple(momentumGridVariables, massGridVariables), + couplingManager); + + // the linear solver + using LinearSolver = + UMFPackIstlSolver>; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = + MultiDomainNewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); + + double preciceDt = couplingParticipant.getMaxTimeStepSize(); + auto dt = preciceDt; + + double vtkTime = 1.0; + + while (couplingParticipant.isCouplingOngoing()) { + couplingParticipant.writeCheckpointIfRequired(); + + couplingParticipant.readQuantityFromOtherSolver(meshName, dataNameV, + dt); + // solve the non-linear system + nonLinearSolver.solve(sol); + + if (writeInterfaceDataToFile) { + writeVelocitiesOnInterfaceToFile( + meshName, Dumux::Fmt::format("ff_interface_velocities_{}", vtkTime), + momentumProblem, *momentumGridVariables, sol[momentumIdx]); + writePressuresOnInterfaceToFile( + meshName, Dumux::Fmt::format("ff_interface_pressures_{}", vtkTime), + momentumProblem, *momentumGridVariables, sol[momentumIdx]); + } + + setInterfacePressures( + momentumProblem, *momentumGridVariables, sol[momentumIdx], meshName, + dataNameP); + couplingParticipant.writeQuantityToOtherSolver(meshName, dataNameP); + couplingParticipant.advance(dt); + preciceDt = couplingParticipant.getMaxTimeStepSize(); + dt = std::min(preciceDt, dt); + + if (!couplingParticipant.readCheckpointIfRequired()) { + vtkTime += 1.; + freeFlowVtkWriter.write(vtkTime); + } + } + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + couplingParticipant.finalize(); + + // print dumux end message + if (mpiHelper.rank() == 0) { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} diff --git a/free-flow-over-porous-media/free-flow-dumux/appl/problem.hh b/free-flow-over-porous-media/free-flow-dumux/appl/problem.hh new file mode 100644 index 000000000..3173d4127 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/appl/problem.hh @@ -0,0 +1,396 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * \brief The free flow sub problem + */ +#ifndef DUMUX_STOKES_SUBPROBLEM_HH +#define DUMUX_STOKES_SUBPROBLEM_HH + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +namespace Dumux { + +/*! + * \brief The free flow sub problem + */ +template +class StokesSubProblem : public BaseProblem { + using ParentType = BaseProblem; + + using GridGeometry = GetPropType; + using GridView = typename GridGeometry::GridView; + using Scalar = GetPropType; + + using ModelTraits = GetPropType; + using Indices = typename ModelTraits::Indices; + + using BoundaryTypes = typename ParentType::BoundaryTypes; + using InitialValues = typename ParentType::InitialValues; + using Sources = typename ParentType::Sources; + using DirichletValues = typename ParentType::DirichletValues; + using BoundaryFluxes = typename ParentType::BoundaryFluxes; + + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = + typename FVElementGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using PrimaryVariables = GetPropType; + + using NumEqVector = Dumux::NumEqVector; + + using FluidSystem = GetPropType; + + using CouplingManager = GetPropType; + + static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; + using VelocityVector = Dune::FieldVector; + +public: + StokesSubProblem(std::shared_ptr gridGeometry, + std::shared_ptr couplingManager) + : ParentType(gridGeometry, couplingManager, "FreeFlow"), + eps_(1e-6), + couplingParticipant_(Dumux::Precice::CouplingAdapter::getInstance()) + { + deltaP_ = getParamFromGroup(this->paramGroup(), + "Problem.PressureDifference"); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Return the sources within the domain. + * + * \param globalPos The global position + */ + Sources sourceAtPos(const GlobalPosition &globalPos) const + { + return Sources(0.0); + } + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment. + * + * \param element The finite element + * \param scvf The sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, + const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + + const auto &globalPos = scvf.center(); + + if constexpr (ParentType::isMomentumProblem()) { + if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) { + values.setAllNeumann(); + } + // slip boundary with coupling interface + else if (onLowerBoundary_(globalPos)) { + values.setAllNeumann(); + // TODO: Check the handling of the corners + if (!onLeftBoundary_(scvf.ipGlobal()) && + !onRightBoundary_(scvf.ipGlobal())) + values.setDirichlet(Indices::velocityYIdx); + } else { + values.setAllDirichlet(); + } + } else // mass subproblem + { + if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) { + values.setNeumann(Indices::conti0EqIdx); + } else { + values.setNeumann(Indices::conti0EqIdx); + } + } + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Dirichlet control volume. + * + * \param globalPos The global position + */ + DirichletValues dirichlet(const Element &element, + const SubControlVolumeFace &scvf) const + { + DirichletValues values(0.0); + values = initialAtPos(scvf.ipGlobal()); + + if constexpr (ParentType::isMomentumProblem()) { + const auto faceId = scvf.index(); + if (couplingParticipant_.isCoupledEntity(faceId)) { + values[Indices::velocityYIdx] = + couplingParticipant_.getScalarQuantityOnFace( + "Free-Flow-Mesh", "Velocity", faceId); + } + } else { + auto pressure = onLeftBoundary_(scvf.ipGlobal()) ? deltaP_ : 0.0; + values[Indices::pressureIdx] = pressure; + } + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFaceVars The element face variables + * \param scvf The boundary sub control volume face + */ + template + BoundaryFluxes neumann(const Element &element, + const FVElementGeometry &fvGeometry, + const ElementVolumeVariables &elemVolVars, + const ElementFluxVariablesCache &elemFluxVarsCache, + const SubControlVolumeFace &scvf) const + { + BoundaryFluxes values(0.0); + + const auto &globalPos = scvf.ipGlobal(); + if constexpr (ParentType::isMomentumProblem()) { +#if DUMUX_VERSION_MAJOR >= 3 & DUMUX_VERSION_MINOR >= 9 + using SlipVelocityPolicy = NavierStokesSlipVelocity< + typename GridGeometry::DiscretizationMethod, + NavierStokes::SlipConditions::BJ>; + using FluxHelper = NavierStokesMomentumBoundaryFlux< + typename GridGeometry::DiscretizationMethod, + SlipVelocityPolicy>; +#else + using FluxHelper = NavierStokesMomentumBoundaryFluxHelper; +#endif + if (onSlipBoundary(fvGeometry, scvf)) { + values += FluxHelper::slipVelocityMomentumFlux( + *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache); + } else if (onLeftBoundary_(globalPos) || + onRightBoundary_(globalPos)) { + auto pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0; + values = FluxHelper::fixedPressureMomentumFlux( + *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, + pressure, /* zeroNormalVelocityGradient = */ true); + } + } else { + using FluxHelper = NavierStokesScalarBoundaryFluxHelper< + AdvectiveFlux>; + if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) { + values = FluxHelper::scalarOutflowFlux( + *this, element, fvGeometry, scvf, elemVolVars); + } else if (onSlipBoundary(fvGeometry, scvf)) { + const Scalar density = + 1000; // TODO how to handle compressible fluids? + // TODO: Use flux helper with outside data? + // TODO: remove hard-coded values index y=1 and dirsign = -1.0 + values[Indices::conti0EqIdx] = + density * this->faceVelocity(element, fvGeometry, scvf)[1] * + -1.0; // scvf.directionSign(); + } + } + return values; + } + + bool onSlipBoundary(const FVElementGeometry &fvGeometry, + const SubControlVolumeFace &scvf) const + { + return onLowerBoundary_(scvf.ipGlobal()); + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param globalPos The global position + */ + InitialValues initialAtPos(const GlobalPosition &globalPos) const + { + InitialValues values(0.0); + + if constexpr (ParentType::isMomentumProblem()) { + // values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]); + } else { + if (onLeftBoundary_(globalPos)) + values[Indices::pressureIdx] = deltaP_; + if (onRightBoundary_(globalPos)) + values[Indices::pressureIdx] = 0.0; + } + + return values; + } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar permeability(const Element &element, + const SubControlVolumeFace &scvf) const + { + return 1e-10; // TODO transfer information or just use constant value + } + + /*! + * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar alphaBJ(const FVElementGeometry &fvGeometry, + const SubControlVolumeFace &scvf) const + { + return 1.0; // TODO transfer information or just use constant value + } + + /*! + * \brief Returns the beta value required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar betaBJ(const FVElementGeometry &fvGeometry, + const SubControlVolumeFace &scvf, + const GlobalPosition &tangentialVector) const + { + return 1e+5; // TODO transfer information or just use constant value + } + + /*! + * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffman). + */ + VelocityVector porousMediumVelocity(const FVElementGeometry &fvGeometry, + const SubControlVolumeFace &scvf) const + { + VelocityVector velocity(0.0); + if constexpr (ParentType::isMomentumProblem()) { + const auto faceId = scvf.index(); + if (couplingParticipant_.isCoupledEntity(faceId)) { + velocity[Indices::velocityYIdx] = + couplingParticipant_.getScalarQuantityOnFace( + "Free-Flow-Mesh", "Velocity", faceId); + } + } + return velocity; + } + + /*! + * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967) + */ + void calculateAnalyticalVelocityX() const + { + analyticalVelocityX_.resize(this->gridGeometry().gridView().size(0)); + + using std::sqrt; + const Scalar dPdX = -deltaP_ / (this->gridGeometry().bBoxMax()[0] - + this->gridGeometry().bBoxMin()[0]); + static const Scalar mu = FluidSystem::viscosity(273.15 + 10, 1e5); + static const Scalar alpha = + getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); + static const Scalar K = + getParam("Darcy.SpatialParams.Permeability"); + static const Scalar sqrtK = sqrt(K); + const Scalar sigma = (this->gridGeometry().bBoxMax()[1] - + this->gridGeometry().bBoxMin()[1]) / + sqrtK; + + const Scalar uB = + -K / (2.0 * mu) * + ((sigma * sigma + 2.0 * alpha * sigma) / (1.0 + alpha * sigma)) * + dPdX; + + for (const auto &element : elements(this->gridGeometry().gridView())) { + const auto eIdx = + this->gridGeometry().gridView().indexSet().index(element); + const Scalar y = element.geometry().center()[1] - + this->gridGeometry().bBoxMin()[1]; + + const Scalar u = + uB * (1.0 + alpha / sqrtK * y) + + 1.0 / (2.0 * mu) * (y * y + 2 * alpha * y * sqrtK) * dPdX; + analyticalVelocityX_[eIdx] = u; + } + } + + /*! + * \brief Get the analytical velocity in x direction + */ + const std::vector &getAnalyticalVelocityX() const + { + if (analyticalVelocityX_.empty()) + calculateAnalyticalVelocityX(); + return analyticalVelocityX_; + } + + // \} + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; + } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; + } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; + } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; + } + + Scalar eps_; + Scalar deltaP_; + + Dumux::Precice::CouplingAdapter &couplingParticipant_; + + mutable std::vector analyticalVelocityX_; +}; +} // namespace Dumux + +#endif // DUMUX_STOKES_SUBPROBLEM_HH diff --git a/free-flow-over-porous-media/free-flow-dumux/appl/properties.hh b/free-flow-over-porous-media/free-flow-dumux/appl/properties.hh new file mode 100644 index 000000000..a2f98b525 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/appl/properties.hh @@ -0,0 +1,119 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * \brief The free flow sub problem + */ +#ifndef DUMUX_STOKES_SUBPROPERTIES_HH +#define DUMUX_STOKES_SUBPROPERTIES_HH + +#ifndef DIMWORLD +#define DIMWORLD 2 +#endif + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "problem.hh" + +namespace Dumux { + +namespace Properties { +// Create new type tags +namespace TTag { +struct FreeFlowModel { +}; +struct FreeFlowSubMomentum { + using InheritsFrom = std:: + tuple; +}; +struct FreeFlowSubMass { + using InheritsFrom = + std::tuple; +}; +} // end namespace TTag + +// the fluid system +template +struct FluidSystem { + using Scalar = GetPropType; + using type = + FluidSystems::OnePLiquid>; +}; + +// Set the grid type +template +struct Grid { + using type = Dune::YaspGrid, + DIMWORLD>>; +}; + +// Set the problem property +template +struct Problem { + using type = + Dumux::StokesSubProblem>; +}; +template +struct Problem { + using type = + Dumux::StokesSubProblem>; +}; + +template +struct EnableGridGeometryCache { + static constexpr bool value = true; +}; +template +struct EnableGridFluxVariablesCache { + static constexpr bool value = true; +}; +template +struct EnableGridVolumeVariablesCache { + static constexpr bool value = true; +}; + +// Define the DuMux coupling manager to couple the momentum and mass subproblems +// of the freeflow participant +template +struct CouplingManager { + using Traits = + MultiDomainTraits; + using type = FreeFlowCouplingManager; +}; +} // end namespace Properties + +} // namespace Dumux + +#endif // DUMUX_STOKES_SUBPROPERTIES_HH diff --git a/free-flow-over-porous-media/free-flow-dumux/clean.sh b/free-flow-over-porous-media/free-flow-dumux/clean.sh new file mode 100755 index 000000000..3e9fd4e51 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_dumux . diff --git a/free-flow-over-porous-media/free-flow-dumux/config.h.cmake b/free-flow-over-porous-media/free-flow-dumux/config.h.cmake new file mode 100644 index 000000000..8a40f6ee9 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/config.h.cmake @@ -0,0 +1,45 @@ +/* begin free-flow-dumux + put the definitions for config.h specific to + your project here. Everything above will be + overwritten +*/ + +/* begin private */ +/* Name of package */ +#define PACKAGE "@DUNE_MOD_NAME@" + +/* Define to the address where bug reports for this package should be sent. */ +#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@" + +/* Define to the full name of this package. */ +#define PACKAGE_NAME "@DUNE_MOD_NAME@" + +/* Define to the full name and version of this package. */ +#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@" + +/* Define to the one symbol short name of this package. */ +#define PACKAGE_TARNAME "@DUNE_MOD_NAME@" + +/* Define to the home page for this package. */ +#define PACKAGE_URL "@DUNE_MOD_URL@" + +/* Define to the version of this package. */ +#define PACKAGE_VERSION "@DUNE_MOD_VERSION@" + +/* end private */ + +/* Define to the version of free-flow-dumux */ +#define FREE_FLOW_DUMUX_VERSION "@FREE_FLOW_DUMUX_VERSION@" + +/* Define to the major version of free-flow-dumux */ +#define FREE_FLOW_DUMUX_VERSION_MAJOR @FREE_FLOW_DUMUX_VERSION_MAJOR@ + +/* Define to the minor version of free-flow-dumux */ +#define FREE_FLOW_DUMUX_VERSION_MINOR @FREE_FLOW_DUMUX_VERSION_MINOR@ + +/* Define to the revision of free-flow-dumux */ +#define FREE_FLOW_DUMUX_VERSION_REVISION @FREE_FLOW_DUMUX_VERSION_REVISION@ + +/* end free-flow-dumux + Everything below here will be overwritten +*/ diff --git a/free-flow-over-porous-media/free-flow-dumux/dune.module b/free-flow-over-porous-media/free-flow-dumux/dune.module new file mode 100644 index 000000000..24d5d5c37 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/dune.module @@ -0,0 +1,12 @@ +################################ +# Dune module information file # +################################ + +# Name of the module +Module: free-flow-dumux +Version: 1.0 +Maintainer: jun.chen@ipvs.uni-stuttgart.de +# Required build dependencies +Depends: dumux-precice +# Optional build dependencies +#Suggests: diff --git a/free-flow-over-porous-media/free-flow-dumux/free-flow-dumux.pc.in b/free-flow-over-porous-media/free-flow-dumux/free-flow-dumux.pc.in new file mode 100644 index 000000000..9565b0e68 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/free-flow-dumux.pc.in @@ -0,0 +1,15 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ +CXX=@CXX@ +CC=@CC@ +DEPENDENCIES=@REQUIRES@ + +Name: @PACKAGE_NAME@ +Version: @VERSION@ +Description: free-flow module +URL: http://dune-project.org/ +Requires: dumux-precice +Libs: -L${libdir} +Cflags: -I${includedir} diff --git a/free-flow-over-porous-media/free-flow-dumux/params.input b/free-flow-over-porous-media/free-flow-dumux/params.input new file mode 100644 index 000000000..6988dd85a --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/params.input @@ -0,0 +1,37 @@ +[FreeFlow] +EnableUnsymmetrizedVelocityGradientForBeaversJoseph = false + +[FreeFlow.Grid] +Verbosity = true +LowerLeft = 0 1 +UpperRight = 1 2 +Cells = 40 40 +Grading1 = 1 + +[Darcy.Grid] +Verbosity = true +LowerLeft = 0 0 +UpperRight = 1 1 +Cells = 40 40 +Grading1 = 1 + +[FreeFlow.Problem] +Name = free-flow-dumux +EnableInertiaTerms = false +PressureDifference = 1e-2 + +[Darcy.SpatialParams] +Permeability = 1e-6 # m^2 +Porosity = 0.4 +AlphaBeaversJoseph = 1.0 + +[Problem] +Name = fvca-iterative +EnableGravity = false +CouplingMode = ReconstructFreeFlowNormalStress + +[Vtk] +AddVelocity = 1 + +[Output] +EnableCSVWriter = false diff --git a/free-flow-over-porous-media/free-flow-dumux/run.sh b/free-flow-over-porous-media/free-flow-dumux/run.sh new file mode 100755 index 000000000..58ae26589 --- /dev/null +++ b/free-flow-over-porous-media/free-flow-dumux/run.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +echo "Free flow solver is launched in serial." +./free_flow_dumux params.input + +close_log \ No newline at end of file diff --git a/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-precice-config-visualization.png b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-precice-config-visualization.png new file mode 100644 index 000000000..24a544ca2 Binary files /dev/null and b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-precice-config-visualization.png differ diff --git a/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-pressure.png b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-pressure.png new file mode 100644 index 000000000..669af0cdf Binary files /dev/null and b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-pressure.png differ diff --git a/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-ux.png b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-ux.png new file mode 100644 index 000000000..3a61f7987 Binary files /dev/null and b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-result-ux.png differ diff --git a/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-setup.png b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-setup.png new file mode 100644 index 000000000..3d25b9e1d Binary files /dev/null and b/free-flow-over-porous-media/images/tutorials-free-flow-over-porous-media-setup.png differ diff --git a/free-flow-over-porous-media/porous-media-dumux/CMakeLists.txt b/free-flow-over-porous-media/porous-media-dumux/CMakeLists.txt new file mode 100644 index 000000000..3e6b44404 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/CMakeLists.txt @@ -0,0 +1,26 @@ +cmake_minimum_required(VERSION 3.13) +project(porous-media-dumux CXX) + +if(NOT (dune-common_DIR OR dune-common_ROOT OR + "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*")) + string(REPLACE ${PROJECT_NAME} dune-common dune-common_DIR + ${PROJECT_BINARY_DIR}) +endif() + +#find dune-common and set the module path +find_package(dune-common REQUIRED) +list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules" + ${dune-common_MODULE_PATH}) + +#include the dune macros +include(DuneMacros) + +# start a dune project with information from dune.module +dune_project() + +dune_enable_all_packages() + +add_subdirectory(appl) + +# finalize the dune project, e.g. generating config.h etc. +finalize_dune_project(GENERATE_CONFIG_H_CMAKE) diff --git a/free-flow-over-porous-media/porous-media-dumux/appl/CMakeLists.txt b/free-flow-over-porous-media/porous-media-dumux/appl/CMakeLists.txt new file mode 100644 index 000000000..ddc9c7f7f --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/appl/CMakeLists.txt @@ -0,0 +1,5 @@ +add_executable(porous_media_dumux main.cc) +target_compile_definitions(porous_media_dumux PUBLIC "ENABLEMONOLITHIC=0") +target_link_libraries(porous_media_dumux PRIVATE dumux-precice) + +add_input_file_links() \ No newline at end of file diff --git a/free-flow-over-porous-media/porous-media-dumux/appl/main.cc b/free-flow-over-porous-media/porous-media-dumux/appl/main.cc new file mode 100644 index 000000000..87fa08e1a --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/appl/main.cc @@ -0,0 +1,300 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief A test problem for the coupled Stokes/Darcy problem (1p) + */ +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include "properties.hh" + +#include "dumux-precice/couplingadapter.hh" + +/*! + * \brief Returns the velocity at the interface using Darcy's law for reconstruction + */ +template +auto velocityAtInterface(const Problem &problem, + const Element &element, + const FVElementGeometry &fvGeometry, + const ElementVolumeVariables &elemVolVars, + const SubControlVolumeFace &scvf, + const ElementFluxVariablesCache &elemFluxVarsCache) +{ + const int phaseIdx = 0; + FluxVariables fluxVars; + fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, + elemFluxVarsCache); + auto upwindTerm = [phaseIdx](const auto &volVars) { + return volVars.mobility(phaseIdx); + }; + const auto scalarVelocity = + fluxVars.advectiveFlux(phaseIdx, upwindTerm) / scvf.area(); + return scalarVelocity; +} + +template +void setInterfaceVelocities(const Problem &problem, + const GridVariables &gridVars, + const SolutionVector &sol, + const std::string meshName, + const std::string dataName) +{ + const auto &gridGeometry = problem.gridGeometry(); + auto fvGeometry = localView(gridGeometry); + auto elemVolVars = localView(gridVars.curGridVolVars()); + auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache()); + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + + for (const auto &element : elements(gridGeometry.gridView())) { + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, sol); + elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); + + for (const auto &scvf : scvfs(fvGeometry)) { + if (couplingParticipant.isCoupledEntity(scvf.index())) { + // TODO: What to do here? + const double v = velocityAtInterface( + problem, element, fvGeometry, elemVolVars, scvf, + elemFluxVarsCache); + couplingParticipant.writeScalarQuantityOnFace( + meshName, dataName, scvf.index(), v); + } + } + } +} + +int main(int argc, char **argv) +{ + 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); + + using DarcyTypeTag = Properties::TTag::DarcyOneP; + + using DarcyGridManager = + Dumux::GridManager>; + DarcyGridManager darcyGridManager; + darcyGridManager.init("Darcy"); // pass parameter group + + // we compute on the leaf grid view + const auto &darcyGridView = darcyGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using DarcyGridGeometry = + GetPropType; + auto darcyGridGeometry = std::make_shared(darcyGridView); + darcyGridGeometry->update(darcyGridManager.grid().leafGridView()); + + using DarcyProblem = GetPropType; + auto darcyProblem = std::make_shared(darcyGridGeometry); + + // the solution vector + using SolutionVector = + GetPropType; + SolutionVector sol; + sol.resize(darcyGridGeometry->numDofs()); + + // Initialize preCICE.Tell preCICE about: + std::string preciceConfigFilename = "../precice-config.xml"; + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + couplingParticipant.announceSolver("Darcy", preciceConfigFilename, + mpiHelper.rank(), mpiHelper.size()); + + const std::string meshName("Porous-Media-Mesh"); + const int dim = couplingParticipant.getMeshDimensions(meshName); + std::cout << dim << " " << int(DarcyGridGeometry::GridView::dimension) + << std::endl; + if (dim != int(DarcyGridGeometry::GridView::dimension)) + DUNE_THROW(Dune::InvalidStateException, "Dimensions do not match"); + + // GET mesh coordinates + const double xMin = + getParam>("Grid.LowerLeft")[0]; + const double xMax = + getParam>("Grid.UpperRight")[0]; + std::vector coords; //( dim * vertexSize ); + std::vector coupledScvfIndices; + + for (const auto &element : elements(darcyGridView)) { + auto fvGeometry = localView(*darcyGridGeometry); + fvGeometry.bindElement(element); + + for (const auto &scvf : scvfs(fvGeometry)) { + static constexpr auto eps = 1e-7; + const auto &pos = scvf.center(); + if (pos[1] > darcyGridGeometry->bBoxMax()[1] - eps) { + if (pos[0] > xMin - eps && pos[0] < xMax + eps) { + coupledScvfIndices.push_back(scvf.index()); + for (const auto p : pos) + coords.push_back(p); + } + } + } + } + + couplingParticipant.setMesh(meshName, coords); + couplingParticipant.createIndexMapping(coupledScvfIndices); + + const std::string dataNameV("Velocity"); + const std::string dataNameP("Pressure"); + couplingParticipant.announceQuantity(meshName, dataNameP); + couplingParticipant.announceQuantity(meshName, dataNameV); + + darcyProblem->applyInitialSolution(sol); + + // the grid variables + using DarcyGridVariables = + GetPropType; + auto darcyGridVariables = + std::make_shared(darcyProblem, darcyGridGeometry); + darcyGridVariables->init(sol); + + // intialize the vtk output module + + VtkOutputModule> + darcyVtkWriter(*darcyGridVariables, sol, darcyProblem->name()); + using DarcyVelocityOutput = + GetPropType; + darcyVtkWriter.addVelocityOutput( + std::make_shared(*darcyGridVariables)); + GetPropType::initOutputModule( + darcyVtkWriter); + darcyVtkWriter.write(0.0); + + using FluxVariables = GetPropType; + if (couplingParticipant.requiresToWriteInitialData()) { + setInterfaceVelocities( + *darcyProblem, *darcyGridVariables, sol, meshName, dataNameV); + couplingParticipant.writeQuantityToOtherSolver(meshName, dataNameV); + } + couplingParticipant.initialize(); + + // initialize checkpointing + couplingParticipant.initializeCheckpoint(sol, *darcyGridVariables); + + // the assembler for a stationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared( + darcyProblem, darcyGridGeometry, darcyGridVariables); + + // the linear solver + using LinearSolver = + UMFPackIstlSolver>; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = Dumux::NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + double preciceDt = couplingParticipant.getMaxTimeStepSize(); + auto dt = preciceDt; + + double vtkTime = 1.0; + + while (couplingParticipant.isCouplingOngoing()) { + couplingParticipant.writeCheckpointIfRequired(); + + couplingParticipant.readQuantityFromOtherSolver(meshName, dataNameP, + dt); + + // solve the non-linear system + nonLinearSolver.solve(sol); + + setInterfaceVelocities( + *darcyProblem, *darcyGridVariables, sol, meshName, dataNameV); + couplingParticipant.writeQuantityToOtherSolver(meshName, dataNameV); + + couplingParticipant.advance(dt); + preciceDt = couplingParticipant.getMaxTimeStepSize(); + dt = std::min(preciceDt, dt); + + if (!couplingParticipant.readCheckpointIfRequired()) { + vtkTime += 1.; + darcyVtkWriter.write(vtkTime); + } + } + couplingParticipant.finalize(); + + //////////////////////////////////////////////////////////// + // 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 diff --git a/free-flow-over-porous-media/porous-media-dumux/appl/problem.hh b/free-flow-over-porous-media/porous-media-dumux/appl/problem.hh new file mode 100644 index 000000000..4f80c0c68 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/appl/problem.hh @@ -0,0 +1,208 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief The porous medium flow sub problem + */ +#ifndef DUMUX_DARCY_SUBPROBLEM_HH +#define DUMUX_DARCY_SUBPROBLEM_HH + +#include + +#include + +#include + +namespace Dumux { + +/*! + * \brief The porous medium flow sub problem + */ +template +class DarcySubProblem : public PorousMediumFlowProblem { + using ParentType = PorousMediumFlowProblem; + using GridView = + typename GetPropType::GridView; + using Scalar = GetPropType; + using PrimaryVariables = GetPropType; + + using NumEqVector = Dumux::NumEqVector; + + using BoundaryTypes = Dumux::BoundaryTypes< + GetPropType::numEq()>; + using VolumeVariables = GetPropType; + using FVElementGeometry = + typename GetPropType::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = + typename FVElementGeometry::SubControlVolumeFace; + using GridGeometry = GetPropType; + + using Indices = + typename GetPropType::Indices; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + DarcySubProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry, "Darcy"), + eps_(1e-7), + couplingParticipant_(Dumux::Precice::CouplingAdapter::getInstance()) + { + } + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary control volume. + * + * \param element The element + * \param scvf The boundary sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, + const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + + // set Neumann BCs to all boundaries first + values.setAllNeumann(); + + const auto faceId = scvf.index(); + if (couplingParticipant_.isCoupledEntity(faceId)) + values.setAllDirichlet(); + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Dirichlet control volume. + * + * \param element The element for which the Dirichlet boundary condition is set + * \param scvf The boundary subcontrolvolumeface + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichlet(const Element &element, + const SubControlVolumeFace &scvf) const + { + // set p = 0 at the bottom + PrimaryVariables values(0.0); + values = initial(element); + + const auto faceId = scvf.index(); + if (couplingParticipant_.isCoupledEntity(faceId)) + values = couplingParticipant_.getScalarQuantityOnFace( + "Porous-Media-Mesh", "Pressure", faceId); + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param scvf The boundary sub control volume face + * + * For this method, the \a values variable stores primary variables. + */ + template + NumEqVector neumann(const Element &element, + const FVElementGeometry &fvGeometry, + const ElementVolumeVariables &elemVolVars, + const ElementFluxVarsCache &elemFluxVarsCache, + const SubControlVolumeFace &scvf) const + { + // no-flow everywhere ... + NumEqVector values(0.0); + return values; + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + /*! + * \brief Evaluate the source term for all phases within a given + * sub-control-volume. + * + * \param element The element for which the source term is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param scv The subcontrolvolume + */ + template + NumEqVector source(const Element &element, + const FVElementGeometry &fvGeometry, + const ElementVolumeVariables &elemVolVars, + const SubControlVolume &scv) const + { + return NumEqVector(0.0); + } + + // \} + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param element The element + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initial(const Element &element) const + { + static const Scalar p = + getParamFromGroup(this->paramGroup(), "Problem.InitialP"); + return PrimaryVariables(p); + } + + // \} + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; + } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; + } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; + } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { + return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; + } + + Scalar eps_; + + Dumux::Precice::CouplingAdapter &couplingParticipant_; +}; +} // namespace Dumux + +#endif // DUMUX_DARCY_SUBPROBLEM_HH diff --git a/free-flow-over-porous-media/porous-media-dumux/appl/properties.hh b/free-flow-over-porous-media/porous-media-dumux/appl/properties.hh new file mode 100644 index 000000000..aeb522e87 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/appl/properties.hh @@ -0,0 +1,81 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief The porous medium flow sub problem + */ +#ifndef DUMUX_DARCY_SUBPROPERTIES_HH +#define DUMUX_DARCY_SUBPROPERTIES_HH + +#ifndef DIMWORLD +#define DIMWORLD 2 +#endif + +#include + +#include + +#include + +#include "spatialparams.hh" + +#include +#include + +#include "problem.hh" + +namespace Dumux::Properties { + +// Create new type tags +namespace TTag { +struct DarcyOneP { + using InheritsFrom = std::tuple; +}; +} // end namespace TTag + +// Set the problem property +template +struct Problem { + using type = Dumux::DarcySubProblem; +}; + +// the fluid system +template +struct FluidSystem { + using Scalar = GetPropType; + using type = + FluidSystems::OnePLiquid>; +}; + +// Set the grid type +template +struct Grid { + using type = Dune::YaspGrid; +}; + +template +struct SpatialParams { + using type = OnePSpatialParams, + GetPropType>; +}; + +} // end namespace Dumux::Properties + +#endif // DUMUX_DARCY_SUBPROPERTIES_HH diff --git a/free-flow-over-porous-media/porous-media-dumux/appl/spatialparams.hh b/free-flow-over-porous-media/porous-media-dumux/appl/spatialparams.hh new file mode 100644 index 000000000..07459cb73 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/appl/spatialparams.hh @@ -0,0 +1,109 @@ +// -*- 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 2 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup OnePTests + * \brief The spatial parameters class for the test problem using the 1p cc model + */ +#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH +#define DUMUX_1P_TEST_SPATIALPARAMS_HH + +#include + +namespace Dumux { +/*! + * \ingroup OnePModel + * + * \brief The spatial parameters class for the test problem using the + * 1p cc model + */ +template +class OnePSpatialParams : public FVPorousMediumFlowSpatialParamsOneP< + FVGridGeometry, + Scalar, + OnePSpatialParams> { + using GridView = typename FVGridGeometry::GridView; + using ParentType = FVPorousMediumFlowSpatialParamsOneP< + FVGridGeometry, + Scalar, + OnePSpatialParams>; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + // export permeability type + using PermeabilityType = Scalar; + + OnePSpatialParams(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) + { + permeability_ = getParam("SpatialParams.Permeability"); + porosity_ = getParam("SpatialParams.Porosity"); + alphaBJ_ = getParam("SpatialParams.AlphaBeaversJoseph"); + } + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. + * + * \param globalPos The global position + * \return the intrinsic permeability + */ + PermeabilityType permeabilityAtPos(const GlobalPosition &globalPos) const + { + return permeability_; + } + + /*! \brief Define the porosity in [-]. + * + * \param globalPos The global position + */ + Scalar porosityAtPos(const GlobalPosition &globalPos) const + { + return porosity_; + } + + /*! \brief Define the Beavers-Joseph coefficient in [-]. + * + * \param globalPos The global position + */ + Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const + { + return alphaBJ_; + } + + /*! + * \brief Return the temperature within the domain in [K]. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperatureAtPos(const GlobalPosition &globalPos) const + { + return 273.15 + 10; // 10°C + } + +private: + Scalar permeability_; + Scalar porosity_; + Scalar alphaBJ_; +}; + +} // namespace Dumux + +#endif diff --git a/free-flow-over-porous-media/porous-media-dumux/clean.sh b/free-flow-over-porous-media/porous-media-dumux/clean.sh new file mode 100755 index 000000000..3e9fd4e51 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_dumux . diff --git a/free-flow-over-porous-media/porous-media-dumux/config.h.cmake b/free-flow-over-porous-media/porous-media-dumux/config.h.cmake new file mode 100644 index 000000000..c508ddec4 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/config.h.cmake @@ -0,0 +1,45 @@ +/* begin porous-media-dumux + put the definitions for config.h specific to + your project here. Everything above will be + overwritten +*/ + +/* begin private */ +/* Name of package */ +#define PACKAGE "@DUNE_MOD_NAME@" + +/* Define to the address where bug reports for this package should be sent. */ +#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@" + +/* Define to the full name of this package. */ +#define PACKAGE_NAME "@DUNE_MOD_NAME@" + +/* Define to the full name and version of this package. */ +#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@" + +/* Define to the one symbol short name of this package. */ +#define PACKAGE_TARNAME "@DUNE_MOD_NAME@" + +/* Define to the home page for this package. */ +#define PACKAGE_URL "@DUNE_MOD_URL@" + +/* Define to the version of this package. */ +#define PACKAGE_VERSION "@DUNE_MOD_VERSION@" + +/* end private */ + +/* Define to the version of porous-media-dumux */ +#define POROUS_MEDIA_DUMUX_VERSION "@POROUS_MEDIA_DUMUX_VERSION@" + +/* Define to the major version of porous-media-dumux */ +#define POROUS_MEDIA_DUMUX_VERSION_MAJOR @POROUS_MEDIA_DUMUX_VERSION_MAJOR@ + +/* Define to the minor version of porous-media-dumux */ +#define POROUS_MEDIA_DUMUX_VERSION_MINOR @POROUS_MEDIA_DUMUX_VERSION_MINOR@ + +/* Define to the revision of porous-media-dumux */ +#define POROUS_MEDIA_DUMUX_VERSION_REVISION @POROUS_MEDIA_DUMUX_VERSION_REVISION@ + +/* end porous-media-dumux + Everything below here will be overwritten +*/ diff --git a/free-flow-over-porous-media/porous-media-dumux/dune.module b/free-flow-over-porous-media/porous-media-dumux/dune.module new file mode 100644 index 000000000..752d3ccce --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/dune.module @@ -0,0 +1,12 @@ +################################ +# Dune module information file # +################################ + +# Name of the module +Module: porous-media-dumux +Version: 1.0 +Maintainer: jun.chen@ipvs.uni-stuttgart.de +# Required build dependencies +Depends: dumux-precice +# Optional build dependencies +#Suggests: diff --git a/free-flow-over-porous-media/porous-media-dumux/params.input b/free-flow-over-porous-media/porous-media-dumux/params.input new file mode 100644 index 000000000..a7cddbae6 --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/params.input @@ -0,0 +1,18 @@ +[Grid] +Verbosity = true +LowerLeft = 0 0 +UpperRight = 1 1 +Cells = 40 40 +Grading1 = 1 + +[Problem] +Name = porous-media-dumux +InitialP = 0.0e-9 + +[SpatialParams] +Permeability = 1e-6 # m^2 +Porosity = 0.4 +AlphaBeaversJoseph = 1.0 + +[Vtk] +AddVelocity = 1 diff --git a/free-flow-over-porous-media/porous-media-dumux/porous-media-dumux.pc.in b/free-flow-over-porous-media/porous-media-dumux/porous-media-dumux.pc.in new file mode 100644 index 000000000..d24c6f29f --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/porous-media-dumux.pc.in @@ -0,0 +1,15 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ +CXX=@CXX@ +CC=@CC@ +DEPENDENCIES=@REQUIRES@ + +Name: @PACKAGE_NAME@ +Version: @VERSION@ +Description: porous-media-dumux module +URL: http://dune-project.org/ +Requires: dumux-precice +Libs: -L${libdir} +Cflags: -I${includedir} diff --git a/free-flow-over-porous-media/porous-media-dumux/run.sh b/free-flow-over-porous-media/porous-media-dumux/run.sh new file mode 100755 index 000000000..db42f655b --- /dev/null +++ b/free-flow-over-porous-media/porous-media-dumux/run.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +echo "Porous medium solver is launched in serial." +./porous_media_dumux params.input + +close_log \ No newline at end of file diff --git a/free-flow-over-porous-media/precice-config.xml b/free-flow-over-porous-media/precice-config.xml new file mode 100644 index 000000000..75f02b4ae --- /dev/null +++ b/free-flow-over-porous-media/precice-config.xml @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/free-flow-over-porous-media/setup-dumux.sh b/free-flow-over-porous-media/setup-dumux.sh new file mode 100755 index 000000000..56fc29fc6 --- /dev/null +++ b/free-flow-over-porous-media/setup-dumux.sh @@ -0,0 +1,35 @@ +#!/usr/bin/env sh +set -e -u + +# This script sets up a DUNE environment in the working directory to solve the two-scale-heat-conduction problem with DuMuX on one or both scales + +# Clean any old leftover dumux or dune folders +rm -rfv dumux/ dumux-adapter/ +rm -rfv dune-*/ +rm -rfv install* + +# Get the DuMuX install script and install it +wget https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/raw/releases/3.7/bin/installdumux.py +DUNE_CONTROL_PATH=. python3 installdumux.py +# clear build directories +cd dumux +rm -r dune-common/build-cmake/dune-env/lib/dunecontrol || true +DUNE_CONTROL_PATH=. ./dune-common/bin/dunecontrol exec rm -r build-cmake +cd .. + +# Take out all the module folders from the dumux/ folder and remove the dumux/ folder +mv dumux dumux-install +mv dumux-install/* ./ +rm -r dumux-install + +# Get additional required DUNE modules +# DuMux-preCICE adapter +git clone https://github.com/precice/dumux-adapter.git +# DUNE SPGrid for periodic boundary conditions +DUNE_CONTROL_PATH=. python3 dumux/bin/installexternal.py spgrid + +# Re-build environment +DUNE_CONTROL_PATH=. ./dune-common/bin/dunecontrol --opts=./dumux/cmake.opts all + +# Compile and move macro-dumux and micro-dumux executables to the participant folder level +./compile-dumux-cases.sh