diff --git a/dumux/flux/porenetwork/advection.hh b/dumux/flux/porenetwork/advection.hh index 14dc8b0da376d6a50b853ed625f86c1d25f590a0..171ad5dc5edb4c1cf6ef635fb772e29d542ad203 100644 --- a/dumux/flux/porenetwork/advection.hh +++ b/dumux/flux/porenetwork/advection.hh @@ -150,6 +150,203 @@ public: } }; +/*! + * \file + * \ingroup PoreNetworkFlux + * \brief Non-creeping flow flux law to describe the advective flux for pore-network models based on + * El-Zehairy et al.(2019). + * https://doi.org/10.1016/j.advwatres.2019.103378 + */ +template +class NonCreepingFlow +{ +public: + //! Export the Scalar type + using Scalar = ScalarT; + + //! Export the creeping flow transmissibility law + using TransmissibilityCreepingFlow = Detail::Transmissibility; + + //! Inherit transmissibility from creeping flow transmissibility law to cache non-creeping flow-related parameters + class Transmissibility: public TransmissibilityCreepingFlow + { + public: + class SinglePhaseCache : public TransmissibilityCreepingFlow::SinglePhaseCache + { + using ParentType = typename TransmissibilityCreepingFlow::SinglePhaseCache; + public: + using Scalar = ScalarT; + + template + void fill(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxVariablesCache& fluxVarsCache, + const int phaseIdx) + { + ParentType::fill(problem, element, fvGeometry,scvf, elemVolVars, fluxVarsCache, phaseIdx); + const auto elemSol = elementSolution(element, elemVolVars, fvGeometry); + + for (const auto& scv : scvs(fvGeometry)) + { + const auto localIdx = scv.indexInElement(); + const Scalar throatToPoreAreaRatio = fluxVarsCache.throatCrossSectionalArea() / problem.spatialParams().poreCrossSectionalArea(element, scv, elemSol); + + // dimensionless momentum coefficient + const Scalar momentumCoefficient = problem.spatialParams().momentumCoefficient(element, scv, elemSol); + + // dimensionless kinetik-energy coefficient + const Scalar kineticEnergyCoefficient = problem.spatialParams().kineticEnergyCoefficient(element, scv, elemSol); + + expansionCoefficient_[localIdx] = getExpansionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient); + contractionCoefficient_[localIdx] = getContractionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient); + } + } + + Scalar expansionCoefficient(int downstreamIdx) const + { return expansionCoefficient_[downstreamIdx]; } + + Scalar contractionCoefficient(int upstreamIdx) const + { return contractionCoefficient_[upstreamIdx]; } + + private: + + Scalar getExpansionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const + { + Scalar expansionCoefficient = throatToPoreAreaRatio * throatToPoreAreaRatio * (2 * momentumCoefficient - kineticEnergyCoefficient) + + kineticEnergyCoefficient - 2 * momentumCoefficient * throatToPoreAreaRatio + - (1 - throatToPoreAreaRatio * throatToPoreAreaRatio); + + return expansionCoefficient; + } + + Scalar getContractionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const + { + const Scalar contractionAreaRatio = getContractionAreaRatio_(throatToPoreAreaRatio); + Scalar contractionCoefficient = (1 - (throatToPoreAreaRatio * throatToPoreAreaRatio * kineticEnergyCoefficient - 2 * momentumCoefficient /*+1-throatToPoreAreaRatio*throatToPoreAreaRatio*/) + * contractionAreaRatio * contractionAreaRatio - 2 * contractionAreaRatio) / (contractionAreaRatio * contractionAreaRatio); + + return contractionCoefficient; + } + + Scalar getContractionAreaRatio_(const Scalar throatToPoreAreaRatio) const + { + return 1-(1-throatToPoreAreaRatio)/(2.08*(1-throatToPoreAreaRatio)+0.5371); + } + + std::array expansionCoefficient_; + std::array contractionCoefficient_; + }; + }; + + /*! + * \brief Returns the advective flux of a fluid phase + * across the given sub-control volume face (corresponding to a pore throat). + * \note The flux is given in N*m, and can be converted + * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme + * for the mobility (1/viscosity) or the product of density and mobility, respectively. + */ + template + static Scalar flux(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const int phaseIdx, + const ElemFluxVarsCache& elemFluxVarsCache) + { + const auto& fluxVarsCache = elemFluxVarsCache[scvf]; + + // Get the inside and outside volume variables + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); + const auto& insideVolVars = elemVolVars[insideScv]; + const auto& outsideVolVars = elemVolVars[outsideScv]; + + // calculate the pressure difference + Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx); + + // add gravity term + static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); + if (enableGravity) + { + const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx); + + // projected distance between pores in gravity field direction + const auto poreToPoreVector = fluxVarsCache.poreToPoreDistance() * scvf.unitOuterNormal(); + const auto& gravityVector = problem.spatialParams().gravity(scvf.center()); + + deltaP += rho * (gravityVector * poreToPoreVector); + } + const Scalar creepingFlowTransmissibility = fluxVarsCache.transmissibility(phaseIdx); + const Scalar throatCrossSectionalArea = fluxVarsCache.throatCrossSectionalArea(); + + assert(scvf.insideScvIdx() == 0); + assert(scvf.outsideScvIdx() == 1); + + // determine the flow direction to predict contraction and expansion + const auto [upstreamIdx, downstreamIdx] = deltaP > 0 ? std::pair(0, 1) : std::pair(1, 0); + + const Scalar contractionCoefficient = fluxVarsCache.singlePhaseFlowVariables().contractionCoefficient(upstreamIdx); + const Scalar expansionCoefficient = fluxVarsCache.singlePhaseFlowVariables().expansionCoefficient(downstreamIdx); + const Scalar mu = elemVolVars[upstreamIdx].viscosity(); + + //! El-Zehairy et al.(2019): Eq.(14): + //! A0Q^2 + B0Q + C0 = 0 where Q is volumetric flowrate, A0 = [Kc + Ke] * rho /[2aij^2], Kc and Ke are contraction and expansion coefficient respectively, + //! aij is cross sectional area of the throat, B0 = 1/K0 (K0 is trnasmissibility used in paper) and C0 = -deltaP + //! In our implementation viscosity of fluid will be included using upwinding later. Thus, creepingFlowTransmissibility = mu * K0 and + //! the volumetric flowrate calculated here q = mu * Q. Substitution of new variables into El-Zehairy et al.(2019), Eq.(14) gives: + //! Aq^2 + Bq + C = 0 where + //! A = A0 / mu^2, B = B0/mu = 1/creepingFlowTransmissibility and C = C0 + //! attention: the q, volumetric flowrate, calculated here is always positive and its sign needs to be determined based on flow direction + //! this approach is taken to prevent the term under the square root (discriminant) becoming negative + const Scalar A = (contractionCoefficient * elemVolVars[upstreamIdx].density() + expansionCoefficient * elemVolVars[downstreamIdx].density()) + / (2.0 * mu * mu * throatCrossSectionalArea * throatCrossSectionalArea); + const Scalar B = 1/ creepingFlowTransmissibility; + const Scalar C = -std::abs(deltaP); + + using std::sqrt; + const auto discriminant = B*B - 4*A*C; + const auto q = (-B + sqrt(discriminant)) / (2*A); + + //! give the volume flowrate proper sign based on flow direction. + if (upstreamIdx == 0) + return q; + else + return -q; + + } + + template + static Scalar calculateTransmissibility(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const ElementVolumeVariables& elemVolVars, + const FluxVariablesCache& fluxVarsCache, + const int phaseIdx) + { + static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1); + return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx); + } + + template + static std::array calculateTransmissibilities(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const FluxVariablesCache& fluxVarsCache) + { + static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1); + const Scalar t = calculateTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, 0); + return {t, -t}; + } + +}; } // end namespace Dumux diff --git a/test/porenetwork/1p/CMakeLists.txt b/test/porenetwork/1p/CMakeLists.txt index 3db08be27a333786a604255f8dae2902e9f984ae..5d5bb3df63a3e6908df287ed13e8f62d46655574 100644 --- a/test/porenetwork/1p/CMakeLists.txt +++ b/test/porenetwork/1p/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(nonisothermal) +add_subdirectory(noncreepingflow) add_input_file_links() dune_symlink_to_source_files(FILES grids) diff --git a/test/porenetwork/1p/noncreepingflow/CMakeLists.txt b/test/porenetwork/1p/noncreepingflow/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..68122a5f5774970018431f00037a4026a6c88ad3 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/CMakeLists.txt @@ -0,0 +1,14 @@ +add_subdirectory(nonisothermal) + +add_input_file_links() + +dune_add_test(NAME test_pnm_1p_noncreeping_flow + SOURCES main.cc + LABELS porenetworkflow + COMPILE_DEFINITIONS ISOTHERMAL=1 + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1p_noncreepingflow-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreepingflow-00000.vtp + --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreeping_flow") diff --git a/test/porenetwork/1p/noncreepingflow/main.cc b/test/porenetwork/1p/noncreepingflow/main.cc new file mode 100644 index 0000000000000000000000000000000000000000..1f9327c481d13e1e07cb1370874aa4c1a8e55658 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/main.cc @@ -0,0 +1,133 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief test for the pore-network model with non-creeping flow + */ +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include "properties.hh" + +int main(int argc, char** argv) +{ + using namespace Dumux; + + using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem; + + // 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 the command line arguments and input file + //////////////////////////////////////////////////////////// + + // parse command line arguments + Parameters::init(argc, argv); + + ////////////////////////////////////////////////////////////////////// + // try to create a grid (from the given grid file or the input file) + ///////////////////////////////////////////////////////////////////// + + using GridManager = Dumux::PoreNetwork::GridManager<3>; + GridManager gridManager; + gridManager.init(); + + // we compute on the leaf grid view + const auto& leafGridView = gridManager.grid().leafGridView(); + auto gridData = gridManager.getGridData(); + + // create the finite volume grid geometry + using GridGeometry = GetPropType; + auto gridGeometry = std::make_shared(leafGridView); + gridGeometry->update(*gridData); + + // the spatial parameters + using SpatialParams = GetPropType; + auto spatialParams = std::make_shared(gridGeometry); + + // the problem (boundary conditions) + using Problem = GetPropType; + auto problem = std::make_shared(gridGeometry, spatialParams); + + // the solution vector + using GridView = typename GridGeometry::GridView; + using SolutionVector = GetPropType; + SolutionVector x(leafGridView.size(GridView::dimension)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = GetPropType; + auto gridVariables = std::make_shared(problem, gridGeometry); + gridVariables->init(x); + + // the assembler + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, gridGeometry, gridVariables); + + // helper class to calculate the boundary flux + const auto boundaryFlux = Dumux::PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x); + + // the linear solver + using LinearSolver = ILU0RestartedGMResBackend; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); + nonLinearSolver.solve(x); + + // output result to vtk + using IOFields = GetPropType; + Dumux::PoreNetwork::VtkOutputModule, SolutionVector> vtkWriter(*gridVariables, x, problem->name()); + IOFields::initOutputModule(vtkWriter); + + vtkWriter.write(0.0); + + if (mpiHelper.rank() == 0) + { + std::cout << "cumulative outflux is: " << boundaryFlux.getFlux("max", 0, true) << std::endl; + DumuxMessage::print(); + Parameters::print(); + } + + return 0; +} diff --git a/test/porenetwork/1p/noncreepingflow/nonisothermal/CMakeLists.txt b/test/porenetwork/1p/noncreepingflow/nonisothermal/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4e1b0aa1d445d76557ccb0cc6b47e3453a20b998 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/nonisothermal/CMakeLists.txt @@ -0,0 +1,12 @@ +add_input_file_links() + +dune_add_test(NAME test_pnm_1pni_noncreeping_flow + SOURCES main.cc + LABELS porenetworkflow + COMPILE_DEFINITIONS ISOTHERMAL=0 + CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )" + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1pni_noncreepingflow-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreepingflow-00017.vtp + --command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreeping_flow") diff --git a/test/porenetwork/1p/noncreepingflow/nonisothermal/main.cc b/test/porenetwork/1p/noncreepingflow/nonisothermal/main.cc new file mode 100644 index 0000000000000000000000000000000000000000..4df31d9f67ddf0216126ea50fa07e4c5a349c179 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/nonisothermal/main.cc @@ -0,0 +1,175 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief test for the non-isothermal pore-network model with non-creeping flow + */ +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "../properties.hh" +#include "../problem.hh" + +int main(int argc, char** argv) +{ + using namespace Dumux; + + using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem; + + // 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 the command line arguments and input file + //////////////////////////////////////////////////////////// + + // parse command line arguments + Parameters::init(argc, argv); + + ////////////////////////////////////////////////////////////////////// + // try to create a grid (from the given grid file or the input file) + ///////////////////////////////////////////////////////////////////// + + using GridManager = Dumux::PoreNetwork::GridManager<3>; + GridManager gridManager; + gridManager.init(); + + // we compute on the leaf grid view + const auto& leafGridView = gridManager.grid().leafGridView(); + auto gridData = gridManager.getGridData(); + + // create the finite volume grid geometry + using GridGeometry = GetPropType; + auto gridGeometry = std::make_shared(leafGridView); + gridGeometry->update(*gridData); + + // the spatial parameters + using SpatialParams = GetPropType; + auto spatialParams = std::make_shared(gridGeometry); + + // the problem (boundary conditions) + using Problem = GetPropType; + auto problem = std::make_shared(gridGeometry, spatialParams); + + // the solution vector + using GridView = typename GridGeometry::GridView; + using SolutionVector = GetPropType; + SolutionVector x(leafGridView.size(GridView::dimension)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = GetPropType; + auto gridVariables = std::make_shared(problem, gridGeometry); + gridVariables->init(x); + + // get some time loop parameters + using Scalar = GetPropType; + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("TimeLoop.DtInitial"); + + // intialize the vtk output module + using IOFields = GetPropType; + Dumux::PoreNetwork::VtkOutputModule, SolutionVector> vtkWriter(*gridVariables, x, problem->name()); + IOFields::initOutputModule(vtkWriter); //! Add model specific output fields + + vtkWriter.write(0.0); + + // instantiate time loop + auto timeLoop = std::make_shared>(0.0, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); + + // the linear solver + using LinearSolver = ILU0BiCGSTABBackend; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // solve the non-linear system + nonLinearSolver.solve(x, *timeLoop); + + // make the new solution the old solution + xOld = x; + gridVariables->advanceTimeStep(); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + // write vtk output + vtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton solver + timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); + + } while (!timeLoop->finished()); + + timeLoop->finalize(leafGridView.comm()); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} diff --git a/test/porenetwork/1p/noncreepingflow/nonisothermal/params.input b/test/porenetwork/1p/noncreepingflow/nonisothermal/params.input new file mode 100644 index 0000000000000000000000000000000000000000..f88b2c64006dc5c265443e6ec75889a07fd80695 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/nonisothermal/params.input @@ -0,0 +1,26 @@ +[TimeLoop] +DtInitial = 1e-3 # [s] +TEnd = 1 # [s] + +[Grid] +File = ./../../grids/pnm_converted.dgf +PoreGeometry = Cube +ThroatCrossSectionShape = Square + +[Problem] +Name = test_pnm_1pni_noncreepingflow # name passed to the output routines +Verbose = true +EnableGravity = false + +[Vtk] +AddVelocity = 0 + +[Component] +SolidHeatCapacity = 1.0 # for compatibility +SolidDensity = 1.0 # for compatibility +SolidThermalConductivity = 1.0 # for compatibility + +[Newton] +MaxSteps = 500 +TargetSteps = 6 +MaxRelativeShift = 1e-8 \ No newline at end of file diff --git a/test/porenetwork/1p/noncreepingflow/params.input b/test/porenetwork/1p/noncreepingflow/params.input new file mode 100644 index 0000000000000000000000000000000000000000..9efa77f6747f2ee97dfdf927dadaf38f40e95b45 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/params.input @@ -0,0 +1,16 @@ +[Grid] +File = ./../grids/pnm_converted.dgf +PoreGeometry = Cube +ThroatCrossSectionShape = Square + +[Problem] +Name = test_pnm_1p_noncreepingflow # name passed to the output routines +EnableGravity = false + +[Vtk] +AddVelocity = 0 + +[Newton] +MaxSteps = 500 +TargetSteps = 6 +MaxRelativeShift = 1e-8 diff --git a/test/porenetwork/1p/noncreepingflow/problem.hh b/test/porenetwork/1p/noncreepingflow/problem.hh new file mode 100644 index 0000000000000000000000000000000000000000..cc6d463f6b5035ef8d576ae158680ca62d27f8b5 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/problem.hh @@ -0,0 +1,147 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief A test problem for the one-phase pore-network model with non-creeping flow. + */ +#ifndef DUMUX_PNM1P_NONCREEPING_PROBLEM_HH +#define DUMUX_PNM1P_NONCREEPING_PROBLEM_HH + +#include +#include + +namespace Dumux { + +template +class PNMOnePNonCreepingProblem : public PorousMediumFlowProblem +{ + using ParentType = PorousMediumFlowProblem; + using Scalar = GetPropType; + using PrimaryVariables = GetPropType; + using BoundaryTypes = Dumux::BoundaryTypes::numEq()>; + using GridGeometry = GetPropType; + using GridView = typename GridGeometry::GridView; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using Indices = typename GetPropType::Indices; + using Labels = GetPropType; + + using Element = typename GridView::template Codim<0>::Entity; + using Vertex = typename GridView::template Codim::Entity; + +public: + template + PNMOnePNonCreepingProblem(std::shared_ptr gridGeometry, std::shared_ptr spatialParams) + : ParentType(gridGeometry, spatialParams) + {} + + /*! + * \name Simulation steering + */ + // \{ + +#if ISOTHERMAL + /*! + * \brief Return the temperature within the domain in [K]. + * + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + // \} +#endif + /*! + * \name Boundary conditions + */ + // \{ + //! Specifies which kind of boundary condition should be used for + //! which equation for a finite volume on the boundary. + BoundaryTypes boundaryTypes(const Element& element, const SubControlVolume& scv) const + { + BoundaryTypes bcTypes; + if (isInletPore_(scv) || isOutletPore_(scv)) + bcTypes.setAllDirichlet(); + else // neuman for the remaining boundaries + bcTypes.setAllNeumann(); +#if !ISOTHERMAL + bcTypes.setDirichlet(Indices::temperatureIdx); +#endif + return bcTypes; + } + + /*! + * \brief Evaluate the boundary conditions for a dirichlet + * control volume. + * + * \param values The dirichlet values for the primary variables + * \param vertex The vertex (pore body) for which the condition is evaluated + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichlet(const Element& element, + const SubControlVolume& scv) const + { + PrimaryVariables values(0.0); + if(isInletPore_(scv)) + values[Indices::pressureIdx] = 1.0e5; + else + values[Indices::pressureIdx] = 0.9e5; +#if !ISOTHERMAL + if(isInletPore_(scv)) + values[Indices::temperatureIdx] = 273.15 +25; + else + values[Indices::temperatureIdx] = 273.15 +20; +#endif + return values; + } + + /*! + * \brief Evaluate the initial value for a control volume. + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initial(const Vertex& vertex) const + { + PrimaryVariables values(0.0); + values[Indices::pressureIdx] = 1e5; +#if !ISOTHERMAL + values[Indices::temperatureIdx] = 273.15 +20; +#endif + return values; + } + + // \} + +private: + + bool isInletPore_(const SubControlVolume& scv) const + { + return this->gridGeometry().poreLabel(scv.dofIndex()) == Labels::inlet; + } + + bool isOutletPore_(const SubControlVolume& scv) const + { + return this->gridGeometry().poreLabel(scv.dofIndex()) == Labels::outlet; + } + +}; +} //end namespace Dumux + +#endif diff --git a/test/porenetwork/1p/noncreepingflow/properties.hh b/test/porenetwork/1p/noncreepingflow/properties.hh new file mode 100644 index 0000000000000000000000000000000000000000..6a46dc4e2baa9276142ff31a1d537dba869c1e50 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/properties.hh @@ -0,0 +1,99 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * + * \brief The properties for the one-phase pore-network model with non-creeping flow. + */ +#ifndef DUMUX_PNM1P_NONCREEPING_PROPERTIES_HH +#define DUMUX_PNM1P_NONCREEPING_PROPERTIES_HH + +#include + +#include + +// Pore network model +#include + +#include +#include +#include + +#include "spatialparams.hh" +#include "problem.hh" + +////////// +// Specify the properties +////////// +namespace Dumux::Properties { + +namespace TTag { +#if ISOTHERMAL +struct PNMOnePNonCreepingProblem { using InheritsFrom = std::tuple; }; +#else +struct PNMOnePNonCreepingProblem { using InheritsFrom = std::tuple; }; +#endif +} + +// Set the problem property +template +struct Problem +{ using type = PNMOnePNonCreepingProblem; }; + +// the fluid system +template +struct FluidSystem +{ + using Scalar = GetPropType; + using type = FluidSystems::OnePLiquid > ; +}; + +// the grid +template +struct Grid +{ using type = Dune::FoamGrid<1, 3>; }; + +template +struct SpatialParams +{ +private: + using GridGeometry = GetPropType; + using Scalar = GetPropType; +public: + using type = Dumux::PoreNetwork::NonCreepingSpatialParams; +}; + +//! The advection type +template +struct AdvectionType +{ +private: + using Scalar = GetPropType; + using TransmissibilityLaw = Dumux::PoreNetwork::TransmissibilityPatzekSilin; +public: + using type = Dumux::PoreNetwork::NonCreepingFlow; +}; + +// use the incompressible local residual (provides analytic jacobian) +template +struct LocalResidual { using type = OnePIncompressibleLocalResidual; }; + +} //end namespace Dumux::Properties + +#endif diff --git a/test/porenetwork/1p/noncreepingflow/spatialparams.hh b/test/porenetwork/1p/noncreepingflow/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..11509dc73cf1367163b27f42c242b482a547f9b7 --- /dev/null +++ b/test/porenetwork/1p/noncreepingflow/spatialparams.hh @@ -0,0 +1,90 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * + * \ingroup SpatialParameters + * \brief Spatial parameters for a pore-network model with non-creeping flow. + */ +#ifndef DUMUX_PNM_NONCREEPING_SPATIAL_PARAMS_1P_HH +#define DUMUX_PNM_NONCREEPING_SPATIAL_PARAMS_1P_HH + +#include +#include +#include + +namespace Dumux::PoreNetwork { + +template +class NonCreepingSpatialParams : public BaseSpatialParams> +{ + using ParentType = BaseSpatialParams>; + + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + +public: + using PermeabilityType = Scalar; + using ParentType::ParentType; + + template + auto poreShapeFactor(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { + return Throat::shapeFactorCircle(); + } + + template + Scalar poreCrossSectionalArea(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { + const auto poreRadius = this->gridGeometry().poreInscribedRadius(scv.dofIndex()); + return poreRadius*poreRadius*M_PI; + } + + template + Scalar poreLength(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { + return this->gridGeometry().poreInscribedRadius(scv.dofIndex()); + } + + // dimensionless kinetic-energy coeffiecient which for non-creeping flow is equal to 1.0 + template + Scalar kineticEnergyCoefficient(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { return 1.0; } + + // dimensionless momentum coeffiecient which for non-creeping flow is equal to 1.0 + template + Scalar momentumCoefficient(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { return 1.0; } + +}; +} // end namespace Dumux::PoreNetwork + +#endif diff --git a/test/references/test_pnm_1p_noncreepingflow-reference.vtp b/test/references/test_pnm_1p_noncreepingflow-reference.vtp new file mode 100644 index 0000000000000000000000000000000000000000..95e7df9ee09faf4d8de40c9c9071810523f216ba --- /dev/null +++ b/test/references/test_pnm_1p_noncreepingflow-reference.vtp @@ -0,0 +1,203 @@ + + + + + + + 98363.1 100000 98013 100000 97885.8 100000 100000 100000 97566.5 98397.5 97540.4 100000 + 97803.4 97137.9 97137.9 97137.9 97919.9 95980.5 95980.5 97137.9 97137.9 95980.5 96677.3 96426.2 + 96097.4 96658.6 96658.6 96039.1 97280.2 97141.1 97788 97185 97137.9 97134.6 96811.5 97072.7 + 97493.4 97597 97559.6 96812.3 95893.8 95937 96658.6 95659.2 93643.4 96079.3 97117.8 97157.5 + 97030.4 96995.5 96835.2 96648.9 97117.3 95717.9 97197.5 96783 96491.7 94963.4 95186.9 95495.2 + 96730.2 95718.8 96065.5 92382.2 95387.1 94739.4 93529.7 95880.2 94623.4 93148.8 91202.5 96379.5 + 96405.5 96397.9 96402.3 96096.2 96888.6 96648.9 96648.9 95955.3 95319 97117.3 95880.2 90000 + 90000 93817.6 93148.8 93148.8 90000 93529.7 + + + 0.00010842 0.0002263 0.00018872 8.6711e-05 0.00014246 0.00017024 0.00014006 5.8137e-05 0.00017043 0.00014712 0.00018289 0.00011213 + 8.001e-05 0.00026652 0.00010984 0.00018149 0.0001191 0.00010862 6.7563e-05 0.00013353 8.8091e-05 0.00011873 8.6177e-05 0.00014165 + 9.7691e-05 0.00014071 8.5166e-05 0.00010965 8.1041e-05 0.00015832 6.1496e-05 0.00011323 0.00033347 6.7184e-05 0.00010274 0.00016438 + 0.00011759 8.0701e-05 0.00015113 0.00013784 9.8161e-05 0.00013 8.1223e-05 0.00014997 0.00013153 0.00020701 0.00010946 8.5023e-05 + 7.8721e-05 0.00017262 0.00012575 0.00010082 0.00015856 0.00012022 0.00019603 9.1747e-05 0.0001146 0.00025862 0.00012908 0.00025963 + 0.00010958 6.7231e-05 0.00012907 0.00015315 9.4275e-05 0.00019972 7.1623e-05 0.00014025 0.00011092 0.00015744 5.0362e-05 0.00022549 + 0.00010874 9.8868e-05 0.00010793 0.00015434 8.0185e-05 6.8354e-05 0.00014551 0.00014525 0.00017218 0.00012177 9.0671e-05 0.00015715 + 0.00012027 0.00015784 0.00022979 0.00010367 5.3866e-05 0.00024136 + + + 6 1 4 1 3 1 5 2 6 3 4 1 + 2 4 2 3 2 6 1 1 1 1 4 2 + 3 2 1 2 3 7 2 4 3 4 6 2 + 3 4 3 5 5 3 4 6 4 5 4 4 + 4 4 5 7 3 3 2 3 2 5 2 4 + 3 2 2 2 3 4 5 4 4 6 2 5 + 4 3 2 2 2 1 1 4 2 1 1 3 + 1 2 1 1 1 1 + + + -1 2 -1 2 -1 2 2 -1 -1 -1 -1 2 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 3 + 3 -1 -1 -1 3 -1 + + + + + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 + + + 2 2 2 2 2 2 2 2 2 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 3 -1 3 -1 -1 -1 -1 3 -1 + -1 -1 3 -1 3 + + + 3.3304e-05 3.384e-05 3.3615e-05 3.3259e-05 3.3423e-05 3.3717e-05 3.3301e-05 3.3963e-05 3.3418e-05 3.4247e-05 3.3272e-05 3.37e-05 + 3.3125e-05 3.3391e-05 3.3434e-05 3.3768e-05 3.3951e-05 3.3734e-05 3.3508e-05 3.3908e-05 3.3903e-05 3.3674e-05 3.3826e-05 3.3772e-05 + 3.3303e-05 3.389e-05 3.3413e-05 3.3431e-05 3.3693e-05 3.353e-05 3.3589e-05 3.3834e-05 3.3987e-05 3.3567e-05 3.3743e-05 3.3553e-05 + 3.371e-05 3.3781e-05 3.3426e-05 3.3726e-05 3.3836e-05 3.3664e-05 3.422e-05 3.3429e-05 3.3386e-05 3.305e-05 3.3941e-05 3.3926e-05 + 3.3606e-05 3.3937e-05 3.3695e-05 3.3731e-05 3.3671e-05 3.3781e-05 3.3667e-05 3.4584e-05 3.3243e-05 3.3015e-05 3.3252e-05 3.3267e-05 + 3.3488e-05 3.3576e-05 3.356e-05 3.3816e-05 3.3764e-05 3.3887e-05 3.4237e-05 3.4047e-05 3.3204e-05 3.2859e-05 3.3938e-05 3.3021e-05 + 3.3656e-05 3.3645e-05 3.4391e-05 3.361e-05 3.3463e-05 3.3712e-05 3.3716e-05 3.3657e-05 3.3429e-05 3.3224e-05 3.374e-05 3.3184e-05 + 3.3288e-05 3.4084e-05 3.3492e-05 3.3586e-05 3.3937e-05 3.3532e-05 3.3981e-05 3.3517e-05 3.3975e-05 3.3845e-05 3.3562e-05 3.3343e-05 + 3.3284e-05 3.3542e-05 3.3486e-05 3.3496e-05 3.3966e-05 3.3533e-05 3.402e-05 3.3455e-05 3.3962e-05 3.3458e-05 3.3693e-05 3.3961e-05 + 3.3494e-05 3.3486e-05 3.4314e-05 3.3954e-05 3.3488e-05 3.3852e-05 3.3512e-05 3.3871e-05 3.4113e-05 3.3098e-05 3.3981e-05 3.4127e-05 + 3.3649e-05 3.4226e-05 3.3685e-05 3.322e-05 3.2903e-05 3.3521e-05 3.3874e-05 3.374e-05 3.3772e-05 3.344e-05 3.3682e-05 3.3086e-05 + 3.3378e-05 3.3359e-05 3.3808e-05 3.3628e-05 3.3247e-05 + + + 6.6609e-05 3.384e-05 3.3615e-05 9.9777e-05 6.6846e-05 0.00010115 9.9902e-05 6.7925e-05 0.00010025 0.00010274 9.9816e-05 3.37e-05 + 9.9376e-05 6.6781e-05 0.0001003 6.7537e-05 6.7903e-05 0.0001012 3.3508e-05 6.7817e-05 0.00010171 6.7347e-05 0.00010148 0.00010132 + 6.6606e-05 3.389e-05 6.6826e-05 3.3431e-05 6.7385e-05 6.706e-05 0.00010077 0.0001015 0.00010196 0.0001007 6.7487e-05 0.00010066 + 6.7421e-05 3.3781e-05 3.3426e-05 0.00010118 6.7672e-05 0.00010099 3.422e-05 6.6858e-05 6.6772e-05 3.305e-05 6.7881e-05 6.7852e-05 + 3.3606e-05 6.7873e-05 3.3695e-05 3.3731e-05 3.3671e-05 3.3781e-05 3.3667e-05 3.4584e-05 6.6486e-05 3.3015e-05 6.6503e-05 9.9802e-05 + 6.6975e-05 6.7153e-05 6.7119e-05 3.3816e-05 6.7529e-05 0.00010166 0.00010271 3.4047e-05 3.3204e-05 6.5718e-05 6.7876e-05 6.6043e-05 + 3.3656e-05 0.00010094 3.4391e-05 3.361e-05 3.3463e-05 3.3712e-05 6.7432e-05 6.7313e-05 3.3429e-05 9.9673e-05 0.00010122 9.9553e-05 + 3.3288e-05 0.00010225 0.00010048 0.00010076 0.00010181 0.0001006 6.7962e-05 3.3517e-05 6.7951e-05 3.3845e-05 3.3562e-05 0.00010003 + 3.3284e-05 3.3542e-05 3.3486e-05 6.6992e-05 6.7931e-05 3.3533e-05 0.00010206 6.691e-05 0.00010189 6.6916e-05 6.7385e-05 3.3961e-05 + 3.3494e-05 6.6972e-05 3.4314e-05 0.00010186 3.3488e-05 0.00010156 6.7024e-05 0.00010161 0.00010234 3.3098e-05 6.7962e-05 6.8254e-05 + 0.00010095 0.00010268 3.3685e-05 3.322e-05 6.5807e-05 3.3521e-05 6.7748e-05 6.748e-05 3.3772e-05 3.344e-05 0.00010104 9.9259e-05 + 6.6757e-05 0.00010008 3.3808e-05 0.00010088 9.974e-05 + + + 1.03853e-14 2.17901e-14 2.13583e-14 6.89564e-15 1.04972e-14 7.18452e-15 6.92186e-15 1.10144e-14 6.99529e-15 7.52868e-15 6.90373e-15 2.15208e-14 + 6.81256e-15 1.04673e-14 7.0052e-15 1.08255e-14 1.10024e-14 7.19547e-15 2.1155e-14 1.09607e-14 7.30394e-15 1.07357e-14 7.25421e-15 7.21938e-15 + 1.03846e-14 2.18868e-14 1.04878e-14 2.10095e-14 1.07538e-14 1.05984e-14 7.10273e-15 7.25965e-15 7.3585e-15 7.08906e-15 1.08015e-14 7.08005e-15 + 1.07698e-14 2.16763e-14 2.10001e-14 7.19006e-15 1.08912e-14 7.15077e-15 2.25324e-14 1.05029e-14 1.04624e-14 2.02994e-14 1.09931e-14 1.09783e-14 + 2.13412e-14 1.09892e-14 2.15112e-14 2.15802e-14 2.14653e-14 2.16763e-14 2.14576e-14 2.32591e-14 1.03285e-14 2.02349e-14 1.03371e-14 6.90054e-15 + 1.05587e-14 1.06419e-14 1.0627e-14 2.17438e-14 1.08216e-14 7.29374e-15 7.52208e-15 2.21924e-14 2.05845e-14 9.97473e-15 1.099e-14 1.01228e-14 + 2.14366e-14 7.13817e-15 2.28719e-14 2.13488e-14 2.10699e-14 2.15438e-14 1.07757e-14 1.07194e-14 2.10058e-14 6.87382e-15 7.19916e-15 6.84902e-15 + 2.07411e-14 7.42176e-15 7.0413e-15 7.10089e-15 7.32608e-15 7.06656e-15 1.10318e-14 2.11721e-14 1.10258e-14 2.17998e-14 2.12575e-14 6.94795e-15 + 2.07336e-14 2.12195e-14 2.11134e-14 1.05662e-14 1.10174e-14 2.12024e-14 7.37989e-15 1.05274e-14 7.34192e-15 1.05302e-14 1.07538e-14 2.20247e-14 + 2.11285e-14 1.05567e-14 2.27186e-14 7.33716e-15 2.11172e-14 7.27081e-15 1.05813e-14 7.28356e-15 7.4405e-15 2.03879e-14 1.10318e-14 1.11746e-14 + 7.14086e-15 7.51462e-15 2.1492e-14 2.06142e-14 1.00147e-14 2.11797e-14 1.09279e-14 1.07987e-14 2.1659e-14 2.10265e-14 7.16253e-15 6.78852e-15 + 1.04547e-14 6.95781e-15 2.17283e-14 7.12799e-15 6.88824e-15 + + + 5.99015e-09 7.53841e-09 7.62382e-09 5.39852e-09 0 6.93452e-09 5.20387e-09 7.67762e-09 0 1.29019e-09 1.29019e-09 0 + 0 2.16052e-09 2.16052e-09 0 0 0 0 1.7922e-09 1.7922e-09 0 3.93211e-10 3.93211e-10 + 3.31915e-09 4.30467e-09 1.68595e-09 4.56227e-09 3.46089e-11 3.46089e-11 0 2.17966e-09 3.27572e-09 2.22224e-09 3.37608e-10 5.39567e-09 + 4.17008e-09 4.86626e-09 5.80559e-10 3.98173e-09 3.03968e-10 0 4.46822e-09 6.22097e-09 3.6823e-09 1.49277e-09 0 1.69068e-09 + 1.12927e-09 7.55509e-09 1.89552e-09 1.39899e-09 6.62929e-10 1.55494e-09 1.76421e-09 5.08805e-10 2.17724e-09 1.47599e-09 1.21006e-09 1.76634e-09 + 2.28712e-09 1.68595e-09 6.84715e-11 4.07889e-10 2.89192e-09 2.45672e-09 1.99284e-09 2.22224e-09 5.25451e-09 8.06972e-12 2.26651e-09 2.28274e-09 + 1.51394e-09 3.97066e-09 2.80652e-09 1.74541e-09 4.12753e-09 2.0229e-09 2.0229e-09 2.26651e-09 7.68804e-10 2.76164e-09 1.99284e-09 1.9949e-09 + 2.99875e-09 2.09168e-09 8.24509e-09 1.50283e-09 4.24535e-09 4.74139e-09 4.82964e-09 4.1571e-09 4.34523e-09 2.99875e-09 4.50158e-09 2.28899e-09 + 9.51145e-09 1.37076e-09 5.50772e-10 3.13159e-09 8.01346e-10 1.49561e-10 3.29597e-11 1.82521e-10 2.20564e-09 3.29597e-11 2.02312e-09 1.72656e-09 + 2.26745e-09 1.72656e-09 0 1.06653e-22 4.09644e-09 2.80652e-09 4.04945e-09 3.1957e-09 2.05881e-09 0 2.09168e-09 2.69012e-09 + 3.48988e-09 1.09428e-22 2.02312e-09 1.16226e-08 2.05881e-09 8.24509e-09 3.92003e-09 5.45577e-09 6.03731e-09 3.92003e-09 8.13293e-09 9.85568e-23 + 1.52274e-22 1.92391e-09 9.51145e-09 0 8.8548e-09 + + + + + 2 2 3 1 1 4 2 1 4 1 1 5 + 2 1 1 1 2 1 1 2 3 2 1 3 + 2 3 4 2 2 4 2 2 2 1 2 4 + 1 2 5 2 3 2 1 3 2 2 3 3 + 1 3 4 2 4 4 1 3 5 1 4 1 + 1 4 2 1 4 4 2 3 5 1 4 5 + 2 5 4 2 4 2 1 5 1 1 5 4 + 3 1 1 3 2 1 3 1 4 3 1 3 + 2 2 1 3 3 1 3 3 3 3 2 3 + 3 2 2 2 2 5 3 2 5 3 2 4 + 2 4 3 3 4 5 3 3 2 3 4 3 + 3 5 4 3 4 2 4 1 2 3 1 2 + 4 2 1 4 1 1 4 2 4 4 2 2 + 4 3 1 4 2 5 3 3 5 4 4 1 + 4 2 3 4 3 5 3 3 4 4 3 3 + 3 4 1 3 5 3 4 3 2 4 4 2 + 3 5 2 3 4 4 4 5 4 4 3 4 + 4 4 3 4 4 5 4 5 3 5 2 2 + 4 1 4 4 1 3 5 2 3 5 2 5 + 5 3 2 5 1 1 5 1 3 5 3 3 + 5 2 4 5 4 1 5 3 4 5 4 5 + 5 5 1 4 4 4 5 3 5 5 4 4 + 5 5 4 5 4 3 + + + + + 0 1 2 3 4 5 0 6 7 6 8 6 + 9 6 10 6 7 11 2 12 8 12 13 14 + 15 14 8 16 0 16 17 18 13 19 15 20 + 17 21 22 23 24 23 25 26 24 27 17 27 + 28 4 29 4 30 2 31 2 29 32 33 32 + 13 32 33 10 34 10 35 10 9 0 34 0 + 36 0 37 9 38 37 22 37 8 37 15 13 + 39 8 40 8 41 22 39 22 42 25 43 40 + 17 40 44 40 45 40 24 17 41 17 29 28 + 46 28 31 47 36 47 48 47 49 47 50 31 + 39 31 38 30 33 29 46 29 51 29 42 29 + 52 36 50 35 53 39 34 39 54 38 55 33 + 55 42 43 42 56 34 51 34 45 34 57 58 + 59 58 50 54 55 60 45 60 52 60 43 45 + 61 45 62 43 63 43 64 43 44 65 66 65 + 67 65 59 65 57 41 61 64 68 64 69 44 + 70 44 46 49 48 49 71 49 48 46 72 73 + 74 73 71 73 50 72 74 72 75 72 76 48 + 71 51 76 51 77 51 78 51 79 51 67 56 + 79 50 57 53 80 53 81 52 71 62 79 59 + 71 59 82 67 75 67 83 57 80 57 84 63 + 85 68 79 68 66 68 69 85 83 69 86 69 + 87 69 66 69 88 70 89 66 83 66 + + + 2 4 6 8 10 12 14 16 18 20 22 24 + 26 28 30 32 34 36 38 40 42 44 46 48 + 50 52 54 56 58 60 62 64 66 68 70 72 + 74 76 78 80 82 84 86 88 90 92 94 96 + 98 100 102 104 106 108 110 112 114 116 118 120 + 122 124 126 128 130 132 134 136 138 140 142 144 + 146 148 150 152 154 156 158 160 162 164 166 168 + 170 172 174 176 178 180 182 184 186 188 190 192 + 194 196 198 200 202 204 206 208 210 212 214 216 + 218 220 222 224 226 228 230 232 234 236 238 240 + 242 244 246 248 250 252 254 256 258 260 262 264 + 266 268 270 272 274 + + + + + diff --git a/test/references/test_pnm_1pni_noncreepingflow-reference.vtp b/test/references/test_pnm_1pni_noncreepingflow-reference.vtp new file mode 100644 index 0000000000000000000000000000000000000000..de7e5ef8f770d07db4cc51cd1a9606d8778321a2 --- /dev/null +++ b/test/references/test_pnm_1pni_noncreepingflow-reference.vtp @@ -0,0 +1,213 @@ + + + + + + + 98363.1 100000 98013 100000 97885.8 100000 100000 100000 97566.5 98397.5 97540.4 100000 + 97803.4 97137.9 97137.9 97137.9 97919.9 95980.5 95980.5 97137.9 97137.9 95980.5 96677.3 96426.2 + 96097.4 96658.6 96658.6 96039.1 97280.2 97141.1 97788 97185 97137.9 97134.6 96811.5 97072.7 + 97493.4 97597 97559.6 96812.3 95893.8 95937 96658.6 95659.2 93643.4 96079.3 97117.8 97157.5 + 97030.4 96995.5 96835.2 96648.9 97117.3 95717.9 97197.5 96783 96491.7 94963.4 95186.9 95495.2 + 96730.2 95718.8 96065.5 92382.2 95387.1 94739.4 93529.7 95880.2 94623.4 93148.8 91202.5 96379.5 + 96405.5 96397.9 96402.3 96096.2 96888.6 96648.9 96648.9 95955.3 95319 97117.3 95880.2 90000 + 90000 93817.6 93148.8 93148.8 90000 93529.7 + + + 0.00010842 0.0002263 0.00018872 8.6711e-05 0.00014246 0.00017024 0.00014006 5.8137e-05 0.00017043 0.00014712 0.00018289 0.00011213 + 8.001e-05 0.00026652 0.00010984 0.00018149 0.0001191 0.00010862 6.7563e-05 0.00013353 8.8091e-05 0.00011873 8.6177e-05 0.00014165 + 9.7691e-05 0.00014071 8.5166e-05 0.00010965 8.1041e-05 0.00015832 6.1496e-05 0.00011323 0.00033347 6.7184e-05 0.00010274 0.00016438 + 0.00011759 8.0701e-05 0.00015113 0.00013784 9.8161e-05 0.00013 8.1223e-05 0.00014997 0.00013153 0.00020701 0.00010946 8.5023e-05 + 7.8721e-05 0.00017262 0.00012575 0.00010082 0.00015856 0.00012022 0.00019603 9.1747e-05 0.0001146 0.00025862 0.00012908 0.00025963 + 0.00010958 6.7231e-05 0.00012907 0.00015315 9.4275e-05 0.00019972 7.1623e-05 0.00014025 0.00011092 0.00015744 5.0362e-05 0.00022549 + 0.00010874 9.8868e-05 0.00010793 0.00015434 8.0185e-05 6.8354e-05 0.00014551 0.00014525 0.00017218 0.00012177 9.0671e-05 0.00015715 + 0.00012027 0.00015784 0.00022979 0.00010367 5.3866e-05 0.00024136 + + + 298.15 298.15 298.15 298.15 298.15 298.15 298.15 297.873 298.15 298.15 298.15 298.15 + 298.15 293.152 293.15 293.149 298.15 298.097 296.515 293.15 293.149 294.97 298.15 298.15 + 298.149 293.829 293.376 298.13 298.15 298.141 298.15 298.15 293.447 298.057 298.15 298.149 + 298.15 298.15 298.15 298.15 298.139 298.141 298.097 298.116 298.131 298.128 298.148 298.15 + 298.149 298.149 298.143 298.139 298.114 298.15 298.135 298.058 298.15 298.124 298.105 298.114 + 298.098 298.128 298.123 298.102 298.124 298.118 298.12 298.131 298.133 298.114 298.119 298.126 + 298.131 297.764 296.673 298.121 298.149 297.6 293.608 298.141 298.146 294.937 294.669 293.15 + 293.15 298.131 293.258 294.615 293.15 293.245 + + + 6 1 4 1 3 1 5 2 6 3 4 1 + 2 4 2 3 2 6 1 1 1 1 4 2 + 3 2 1 2 3 7 2 4 3 4 6 2 + 3 4 3 5 5 3 4 6 4 5 4 4 + 4 4 5 7 3 3 2 3 2 5 2 4 + 3 2 2 2 3 4 5 4 4 6 2 5 + 4 3 2 2 2 1 1 4 2 1 1 3 + 1 2 1 1 1 1 + + + -1 2 -1 2 -1 2 2 -1 -1 -1 -1 2 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 3 + 3 -1 -1 -1 3 -1 + + + + + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 + + + 2 2 2 2 2 2 2 2 2 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + -1 -1 -1 3 -1 3 -1 -1 -1 -1 3 -1 + -1 -1 3 -1 3 + + + 3.3304e-05 3.384e-05 3.3615e-05 3.3259e-05 3.3423e-05 3.3717e-05 3.3301e-05 3.3963e-05 3.3418e-05 3.4247e-05 3.3272e-05 3.37e-05 + 3.3125e-05 3.3391e-05 3.3434e-05 3.3768e-05 3.3951e-05 3.3734e-05 3.3508e-05 3.3908e-05 3.3903e-05 3.3674e-05 3.3826e-05 3.3772e-05 + 3.3303e-05 3.389e-05 3.3413e-05 3.3431e-05 3.3693e-05 3.353e-05 3.3589e-05 3.3834e-05 3.3987e-05 3.3567e-05 3.3743e-05 3.3553e-05 + 3.371e-05 3.3781e-05 3.3426e-05 3.3726e-05 3.3836e-05 3.3664e-05 3.422e-05 3.3429e-05 3.3386e-05 3.305e-05 3.3941e-05 3.3926e-05 + 3.3606e-05 3.3937e-05 3.3695e-05 3.3731e-05 3.3671e-05 3.3781e-05 3.3667e-05 3.4584e-05 3.3243e-05 3.3015e-05 3.3252e-05 3.3267e-05 + 3.3488e-05 3.3576e-05 3.356e-05 3.3816e-05 3.3764e-05 3.3887e-05 3.4237e-05 3.4047e-05 3.3204e-05 3.2859e-05 3.3938e-05 3.3021e-05 + 3.3656e-05 3.3645e-05 3.4391e-05 3.361e-05 3.3463e-05 3.3712e-05 3.3716e-05 3.3657e-05 3.3429e-05 3.3224e-05 3.374e-05 3.3184e-05 + 3.3288e-05 3.4084e-05 3.3492e-05 3.3586e-05 3.3937e-05 3.3532e-05 3.3981e-05 3.3517e-05 3.3975e-05 3.3845e-05 3.3562e-05 3.3343e-05 + 3.3284e-05 3.3542e-05 3.3486e-05 3.3496e-05 3.3966e-05 3.3533e-05 3.402e-05 3.3455e-05 3.3962e-05 3.3458e-05 3.3693e-05 3.3961e-05 + 3.3494e-05 3.3486e-05 3.4314e-05 3.3954e-05 3.3488e-05 3.3852e-05 3.3512e-05 3.3871e-05 3.4113e-05 3.3098e-05 3.3981e-05 3.4127e-05 + 3.3649e-05 3.4226e-05 3.3685e-05 3.322e-05 3.2903e-05 3.3521e-05 3.3874e-05 3.374e-05 3.3772e-05 3.344e-05 3.3682e-05 3.3086e-05 + 3.3378e-05 3.3359e-05 3.3808e-05 3.3628e-05 3.3247e-05 + + + 6.6609e-05 3.384e-05 3.3615e-05 9.9777e-05 6.6846e-05 0.00010115 9.9902e-05 6.7925e-05 0.00010025 0.00010274 9.9816e-05 3.37e-05 + 9.9376e-05 6.6781e-05 0.0001003 6.7537e-05 6.7903e-05 0.0001012 3.3508e-05 6.7817e-05 0.00010171 6.7347e-05 0.00010148 0.00010132 + 6.6606e-05 3.389e-05 6.6826e-05 3.3431e-05 6.7385e-05 6.706e-05 0.00010077 0.0001015 0.00010196 0.0001007 6.7487e-05 0.00010066 + 6.7421e-05 3.3781e-05 3.3426e-05 0.00010118 6.7672e-05 0.00010099 3.422e-05 6.6858e-05 6.6772e-05 3.305e-05 6.7881e-05 6.7852e-05 + 3.3606e-05 6.7873e-05 3.3695e-05 3.3731e-05 3.3671e-05 3.3781e-05 3.3667e-05 3.4584e-05 6.6486e-05 3.3015e-05 6.6503e-05 9.9802e-05 + 6.6975e-05 6.7153e-05 6.7119e-05 3.3816e-05 6.7529e-05 0.00010166 0.00010271 3.4047e-05 3.3204e-05 6.5718e-05 6.7876e-05 6.6043e-05 + 3.3656e-05 0.00010094 3.4391e-05 3.361e-05 3.3463e-05 3.3712e-05 6.7432e-05 6.7313e-05 3.3429e-05 9.9673e-05 0.00010122 9.9553e-05 + 3.3288e-05 0.00010225 0.00010048 0.00010076 0.00010181 0.0001006 6.7962e-05 3.3517e-05 6.7951e-05 3.3845e-05 3.3562e-05 0.00010003 + 3.3284e-05 3.3542e-05 3.3486e-05 6.6992e-05 6.7931e-05 3.3533e-05 0.00010206 6.691e-05 0.00010189 6.6916e-05 6.7385e-05 3.3961e-05 + 3.3494e-05 6.6972e-05 3.4314e-05 0.00010186 3.3488e-05 0.00010156 6.7024e-05 0.00010161 0.00010234 3.3098e-05 6.7962e-05 6.8254e-05 + 0.00010095 0.00010268 3.3685e-05 3.322e-05 6.5807e-05 3.3521e-05 6.7748e-05 6.748e-05 3.3772e-05 3.344e-05 0.00010104 9.9259e-05 + 6.6757e-05 0.00010008 3.3808e-05 0.00010088 9.974e-05 + + + 1.03853e-14 2.17901e-14 2.13583e-14 6.89564e-15 1.04972e-14 7.18452e-15 6.92186e-15 1.10144e-14 6.99529e-15 7.52868e-15 6.90373e-15 2.15208e-14 + 6.81256e-15 1.04673e-14 7.0052e-15 1.08255e-14 1.10024e-14 7.19547e-15 2.1155e-14 1.09607e-14 7.30394e-15 1.07357e-14 7.25421e-15 7.21938e-15 + 1.03846e-14 2.18868e-14 1.04878e-14 2.10095e-14 1.07538e-14 1.05984e-14 7.10273e-15 7.25965e-15 7.3585e-15 7.08906e-15 1.08015e-14 7.08005e-15 + 1.07698e-14 2.16763e-14 2.10001e-14 7.19006e-15 1.08912e-14 7.15077e-15 2.25324e-14 1.05029e-14 1.04624e-14 2.02994e-14 1.09931e-14 1.09783e-14 + 2.13412e-14 1.09892e-14 2.15112e-14 2.15802e-14 2.14653e-14 2.16763e-14 2.14576e-14 2.32591e-14 1.03285e-14 2.02349e-14 1.03371e-14 6.90054e-15 + 1.05587e-14 1.06419e-14 1.0627e-14 2.17438e-14 1.08216e-14 7.29374e-15 7.52208e-15 2.21924e-14 2.05845e-14 9.97473e-15 1.099e-14 1.01228e-14 + 2.14366e-14 7.13817e-15 2.28719e-14 2.13488e-14 2.10699e-14 2.15438e-14 1.07757e-14 1.07194e-14 2.10058e-14 6.87382e-15 7.19916e-15 6.84902e-15 + 2.07411e-14 7.42176e-15 7.0413e-15 7.10089e-15 7.32608e-15 7.06656e-15 1.10318e-14 2.11721e-14 1.10258e-14 2.17998e-14 2.12575e-14 6.94795e-15 + 2.07336e-14 2.12195e-14 2.11134e-14 1.05662e-14 1.10174e-14 2.12024e-14 7.37989e-15 1.05274e-14 7.34192e-15 1.05302e-14 1.07538e-14 2.20247e-14 + 2.11285e-14 1.05567e-14 2.27186e-14 7.33716e-15 2.11172e-14 7.27081e-15 1.05813e-14 7.28356e-15 7.4405e-15 2.03879e-14 1.10318e-14 1.11746e-14 + 7.14086e-15 7.51462e-15 2.1492e-14 2.06142e-14 1.00147e-14 2.11797e-14 1.09279e-14 1.07987e-14 2.1659e-14 2.10265e-14 7.16253e-15 6.78852e-15 + 1.04547e-14 6.95781e-15 2.17283e-14 7.12799e-15 6.88824e-15 + + + 5.99015e-09 7.53841e-09 7.62382e-09 5.39852e-09 0 6.93452e-09 5.20387e-09 7.67762e-09 0 1.29019e-09 1.29019e-09 0 + 0 2.16052e-09 2.16052e-09 0 0 0 0 1.7922e-09 1.7922e-09 0 3.93211e-10 3.93211e-10 + 3.31915e-09 4.30467e-09 1.68595e-09 4.56227e-09 3.46089e-11 3.46089e-11 0 2.17966e-09 3.27572e-09 2.22224e-09 3.37608e-10 5.39567e-09 + 4.17008e-09 4.86626e-09 5.80559e-10 3.98173e-09 3.03968e-10 0 4.46822e-09 6.22097e-09 3.6823e-09 1.49277e-09 0 1.69068e-09 + 1.12927e-09 7.55509e-09 1.89552e-09 1.39899e-09 6.62929e-10 1.55494e-09 1.76421e-09 5.08805e-10 2.17724e-09 1.47599e-09 1.21006e-09 1.76634e-09 + 2.28712e-09 1.68595e-09 6.84715e-11 4.07889e-10 2.89192e-09 2.45672e-09 1.99284e-09 2.22224e-09 5.25451e-09 8.06972e-12 2.26651e-09 2.28274e-09 + 1.51394e-09 3.97066e-09 2.80652e-09 1.74541e-09 4.12753e-09 2.0229e-09 2.0229e-09 2.26651e-09 7.68804e-10 2.76164e-09 1.99284e-09 1.9949e-09 + 2.99875e-09 2.09168e-09 8.24509e-09 1.50283e-09 4.24535e-09 4.74139e-09 4.82964e-09 4.1571e-09 4.34523e-09 2.99875e-09 4.50158e-09 2.28899e-09 + 9.51145e-09 1.37076e-09 5.50772e-10 3.13159e-09 8.01346e-10 1.49561e-10 3.29597e-11 1.82521e-10 2.20564e-09 3.29597e-11 2.02312e-09 1.72656e-09 + 2.26745e-09 1.72656e-09 0 0 4.09644e-09 2.80652e-09 4.04945e-09 3.1957e-09 2.05881e-09 0 2.09168e-09 2.69012e-09 + 3.48988e-09 0 2.02312e-09 1.16226e-08 2.05881e-09 8.24509e-09 3.92003e-09 5.45577e-09 6.03731e-09 3.92003e-09 8.13293e-09 0 + 0 1.92391e-09 9.51145e-09 0 8.8548e-09 + + + + + 2 2 3 1 1 4 2 1 4 1 1 5 + 2 1 1 1 2 1 1 2 3 2 1 3 + 2 3 4 2 2 4 2 2 2 1 2 4 + 1 2 5 2 3 2 1 3 2 2 3 3 + 1 3 4 2 4 4 1 3 5 1 4 1 + 1 4 2 1 4 4 2 3 5 1 4 5 + 2 5 4 2 4 2 1 5 1 1 5 4 + 3 1 1 3 2 1 3 1 4 3 1 3 + 2 2 1 3 3 1 3 3 3 3 2 3 + 3 2 2 2 2 5 3 2 5 3 2 4 + 2 4 3 3 4 5 3 3 2 3 4 3 + 3 5 4 3 4 2 4 1 2 3 1 2 + 4 2 1 4 1 1 4 2 4 4 2 2 + 4 3 1 4 2 5 3 3 5 4 4 1 + 4 2 3 4 3 5 3 3 4 4 3 3 + 3 4 1 3 5 3 4 3 2 4 4 2 + 3 5 2 3 4 4 4 5 4 4 3 4 + 4 4 3 4 4 5 4 5 3 5 2 2 + 4 1 4 4 1 3 5 2 3 5 2 5 + 5 3 2 5 1 1 5 1 3 5 3 3 + 5 2 4 5 4 1 5 3 4 5 4 5 + 5 5 1 4 4 4 5 3 5 5 4 4 + 5 5 4 5 4 3 + + + + + 0 1 2 3 4 5 0 6 7 6 8 6 + 9 6 10 6 7 11 2 12 8 12 13 14 + 15 14 8 16 0 16 17 18 13 19 15 20 + 17 21 22 23 24 23 25 26 24 27 17 27 + 28 4 29 4 30 2 31 2 29 32 33 32 + 13 32 33 10 34 10 35 10 9 0 34 0 + 36 0 37 9 38 37 22 37 8 37 15 13 + 39 8 40 8 41 22 39 22 42 25 43 40 + 17 40 44 40 45 40 24 17 41 17 29 28 + 46 28 31 47 36 47 48 47 49 47 50 31 + 39 31 38 30 33 29 46 29 51 29 42 29 + 52 36 50 35 53 39 34 39 54 38 55 33 + 55 42 43 42 56 34 51 34 45 34 57 58 + 59 58 50 54 55 60 45 60 52 60 43 45 + 61 45 62 43 63 43 64 43 44 65 66 65 + 67 65 59 65 57 41 61 64 68 64 69 44 + 70 44 46 49 48 49 71 49 48 46 72 73 + 74 73 71 73 50 72 74 72 75 72 76 48 + 71 51 76 51 77 51 78 51 79 51 67 56 + 79 50 57 53 80 53 81 52 71 62 79 59 + 71 59 82 67 75 67 83 57 80 57 84 63 + 85 68 79 68 66 68 69 85 83 69 86 69 + 87 69 66 69 88 70 89 66 83 66 + + + 2 4 6 8 10 12 14 16 18 20 22 24 + 26 28 30 32 34 36 38 40 42 44 46 48 + 50 52 54 56 58 60 62 64 66 68 70 72 + 74 76 78 80 82 84 86 88 90 92 94 96 + 98 100 102 104 106 108 110 112 114 116 118 120 + 122 124 126 128 130 132 134 136 138 140 142 144 + 146 148 150 152 154 156 158 160 162 164 166 168 + 170 172 174 176 178 180 182 184 186 188 190 192 + 194 196 198 200 202 204 206 208 210 212 214 216 + 218 220 222 224 226 228 230 232 234 236 238 240 + 242 244 246 248 250 252 254 256 258 260 262 264 + 266 268 270 272 274 + + + + +