Merge branch 'feature/structural-dynamics' into 'master'

Feature/structural dynamics

See merge request !3878

dumux/geomechanics/hyperelastic/localresidual.hh

@@ -14,6 +14,8 @@
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/std/type_traits.hh>
+#include <dune/common/exceptions.hh>
 
 #include <dumux/common/math.hh>
 #include <dumux/common/properties.hh>
@@ -22,6 +24,26 @@
 
 namespace Dumux {
 
+#ifndef DOXYGEN
+namespace Detail::Hyperelastic {
+// helper struct detecting if the user-defined problem class has a solidDensity method
+template <typename T, typename ...Ts>
+using SolidDensityDetector = decltype(std::declval<T>().solidDensity(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool hasSolidDensity()
+{ return Dune::Std::is_detected<SolidDensityDetector, T, Args...>::value; }
+
+// helper struct detecting if the user-defined problem class has a acceleration method
+template <typename T, typename ...Ts>
+using AccelerationDetector = decltype(std::declval<T>().acceleration(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool hasAcceleration()
+{ return Dune::Std::is_detected<AccelerationDetector, T, Args...>::value; }
+} // end namespace Detail::Hyperelastic
+#endif
+
 /*!
  * \ingroup Hyperelastic
  * \brief Local residual for the hyperelastic model
@@ -51,15 +73,36 @@ public:
     using ParentType::ParentType;
     using ElementResidualVector = typename ParentType::ElementResidualVector;
 
+    using ParentType::evalStorage;
+
     /*!
-     * \brief Evaluate the rate of change of all conserved quantities
+     * \brief Evaluate the storage contribution to the residual: rho*d^2u/dt^2
      */
-    NumEqVector computeStorage(const Problem& problem,
-                               const SubControlVolume& scv,
-                               const VolumeVariables& volVars) const
+    void evalStorage(ElementResidualVector& residual,
+                     const Problem& problem,
+                     const Element& element,
+                     const FVElementGeometry& fvGeometry,
+                     const ElementVolumeVariables& prevElemVolVars,
+                     const ElementVolumeVariables& curElemVolVars,
+                     const SubControlVolume& scv) const
     {
-        // only the static model is implemented so far
-        return NumEqVector(0.0);
+        const auto& curVolVars = curElemVolVars[scv];
+
+        if constexpr (Detail::Hyperelastic::hasSolidDensity<Problem, const Element&, const SubControlVolume&>()
+                      && Detail::Hyperelastic::hasAcceleration<Problem, const Element&, const SubControlVolume&, Scalar, decltype(curVolVars.displacement())>())
+        {
+            const auto& curVolVars = curElemVolVars[scv];
+            NumEqVector storage = problem.solidDensity(element, scv)*problem.acceleration(
+                element, scv, this->timeLoop().timeStepSize(), curVolVars.displacement()
+            );
+            storage *= curVolVars.extrusionFactor()*Extrusion::volume(fvGeometry, scv);
+            residual[scv.localDofIndex()] += storage;
+        }
+        else
+            DUNE_THROW(Dune::InvalidStateException,
+                "Problem class must implement solidDensity and acceleration"
+                " methods to evaluate storage term"
+            );
     }
 
     /*!

dumux/geomechanics/hyperelastic/spatialparams.hh

@@ -51,6 +51,24 @@ private:
     Scalar E_, nu_, mu_, K_, lambda_;
 };
 
+template<class GridGeometry, class Scalar>
+class DefaultDynamicHyperelasticSpatialParams
+: public DefaultHyperelasticSpatialParams<GridGeometry, Scalar>
+{
+    using ParentType = DefaultHyperelasticSpatialParams<GridGeometry, Scalar>;
+public:
+    DefaultDynamicHyperelasticSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+    : ParentType(gridGeometry)
+    , rho_(getParam<Scalar>("SpatialParams.SolidDensity"))
+    {}
+
+    Scalar solidDensity() const
+    { return rho_; }
+
+private:
+    Scalar rho_;
+};
+
 } // end namespace Dumux
 
 #endif

test/geomechanics/CMakeLists.txt

@@ -3,4 +3,5 @@
 
 add_subdirectory(elastic)
 add_subdirectory(hyperelastic)
+add_subdirectory(hyperelastic_dynamic)
 add_subdirectory(poroelastic)

test/geomechanics/hyperelastic_dynamic/CMakeLists.txt

+# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+dune_symlink_to_source_files(FILES params.input solid.msh)
+dumux_add_test(
+    NAME test_saint_venant_kirchhoff_dynamic
+    LABELS geomechanics elastic
+    SOURCES main.cc
+    CMAKE_GUARD "( HAVE_UMFPACK AND dune-alugrid_FOUND )"
+    COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+    CMD_ARGS  --script fuzzy
+              --files ${CMAKE_SOURCE_DIR}/test/references/test_saint_venant_kirchhoff_dynamic-reference.vtu
+                      ${CMAKE_CURRENT_BINARY_DIR}/test_saint_venant_kirchhoff_dynamic-00002.vtu
+              --command "${CMAKE_CURRENT_BINARY_DIR}/test_saint_venant_kirchhoff_dynamic params.input -TimeLoop.TEnd 1 -VTKOutput.Every 50"
+)

test/geomechanics/hyperelastic_dynamic/main.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+#include <config.h>
+#include <iostream>
+
+#include <dumux/common/initialize.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/istlsolvers.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/experimental/timestepping/newmarkbeta.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_alu.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = Properties::TTag::DynamicHyperelasticityTest;
+
+    // maybe initialize MPI and/or multithreading backend
+    Dumux::initialize(argc, argv);
+
+    // initialize parameter tree
+    Parameters::init(argc, argv);
+
+    using Grid = GetPropType<TypeTag, Properties::Grid>;
+    GridManager<Grid> gridManager;
+    gridManager.init();
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+
+    // the problem (initial and boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(gridGeometry);
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x(gridGeometry->numDofs());
+    x = 0.0;
+
+    // the grid variables
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+    gridVariables->init(x);
+
+    // initialize the vtk output module
+    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+    vtkWriter.addVolumeVariable([](const auto& v){
+        return Dune::FieldVector<double, 2>{v.displacement(0), v.displacement(1)};
+    }, "d");
+    vtkWriter.write(0.0);
+
+    // solution at previous time
+    auto xOld = x;
+
+    // get some time loop parameters
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
+
+    // the time stepping scheme
+    auto newmarkBeta = std::make_shared<Experimental::NewmarkBeta<Scalar, SolutionVector>>();
+    newmarkBeta->initialize(x);
+    problem->setNewmarkScheme(newmarkBeta);
+
+    // the assembler with time loop for a transient problem
+    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
+
+    // the linear solver
+    using LAT = LinearAlgebraTraitsFromAssembler<Assembler>;
+    using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LAT>;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    auto nonLinearSolver = std::make_shared<NewtonSolver>(assembler, linearSolver);
+
+    const int vtkInterval = getParam<int>("VTKOutput.Every", 5);
+
+    // time loop
+    timeLoop->start(); do
+    {
+        // linearize & solve
+        nonLinearSolver->solve(x);
+
+        // update the solution in the time stepping scheme
+        newmarkBeta->update(dt, x);
+
+        // 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
+        if (timeLoop->timeStepIndex() % vtkInterval == 0 || timeLoop->finished())
+            vtkWriter.write(timeLoop->time());
+
+        // report statistics of this time step
+        timeLoop->reportTimeStep();
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    return 0;
+}

test/geomechanics/hyperelastic_dynamic/params.input

+[TimeLoop]
+TEnd = 10
+DtInitial = 0.01
+
+[Grid]
+File = solid.msh
+
+[Problem]
+Name = test_saint_venant_kirchhoff_dynamic
+Gravity = 2.0
+
+[SpatialParams]
+YoungsModulus = 1.4e6
+PoissonRatio = 0.4
+SolidDensity = 1000
+
+[Vtk]
+Precision = Float64
+CoordPrecision = Float64
+AddProcessRank = false
+
+[Newton]
+EnableDynamicOutput = false
+MaxRelativeShift = 1e-6

test/geomechanics/hyperelastic_dynamic/problem.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+#ifndef DUMUX_DYNAMIC_HYPERELASTICITY_TEST_PROBLEM_HH
+#define DUMUX_DYNAMIC_HYPERELASTICITY_TEST_PROBLEM_HH
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/fvproblemwithspatialparams.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+// This test case corresponds to the CSM benchmark problem
+// of a hyperelastic solid material under gravity given in
+// Turek, S., Hron, J. (2006). "Proposal for Numerical Benchmarking of Fluid-Structure Interaction
+// between an Elastic Object and Laminar Incompressible Flow."
+// https://doi.org/10.1007/3-540-34596-5_15
+template<class TypeTag>
+class DynamicHyperelasticityProblem : public FVProblemWithSpatialParams<TypeTag>
+{
+    using ParentType = FVProblemWithSpatialParams<TypeTag>;
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using SubControlVolume = typename GridGeometry::SubControlVolume;
+    using GridView = typename GridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+    static constexpr int dim = GridView::dimension;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+    using Tensor = Dune::FieldMatrix<Scalar, dim, dim>;
+
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+
+public:
+    DynamicHyperelasticityProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+    : ParentType(gridGeometry)
+    , gravity_(getParam<Scalar>("Problem.Gravity"))
+    {}
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param globalPos The position of the center of the finite volume
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+        const auto r = std::hypot(globalPos[0]-0.2, globalPos[1]-0.2);
+        if (r < 0.05 + eps_)
+        {
+            values.setDirichlet(0);
+            values.setDirichlet(1);
+        }
+        else
+            values.setAllNeumann();
+        return values;
+    }
+
+    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+    { return PrimaryVariables(0.0); }
+
+    NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+    {
+        // gravity forcing
+        return {0.0, -this->spatialParams().solidDensity()*gravity_};
+    }
+
+    // Saint-Venant Kirchhoff material
+    Tensor firstPiolaKirchhoffStressTensor(const Tensor& F) const
+    {
+        // material parameters
+        const auto mu = this->spatialParams().shearModulus();
+        const auto lambda = this->spatialParams().firstLameParameter();
+
+        // Lagrangian Green strain E = 1/2*(F^T F - I)
+        auto E = multiplyMatrices(transpose(F), F);
+        E *= 0.5;
+        Scalar trace = 0.0;
+        for (int i = 0; i < dim; ++i)
+        {
+            E[i][i] -= 0.5;
+            trace += E[i][i];
+        }
+
+        // 2nd Piola Kirchhoff stress tensor S = λtr(E)I + 2µE
+        auto& S = E;
+        S *= 2*mu;
+        for (int i = 0; i < dim; ++i)
+            S[i][i] += lambda*trace;
+
+        // 1st Piola Kirchhoff stress tensor P = FS
+        return multiplyMatrices(F, S);
+    }
+
+    // the following methods are needed to solve structural dynamics
+    // with the Newmark-beta time integration scheme
+
+    // we use the Newmark scheme for time integration
+    void setNewmarkScheme(std::shared_ptr<const Experimental::NewmarkBeta<Scalar, SolutionVector>> newmark)
+    { newmark_ = std::move(newmark); }
+
+    // the effective density of the solid material
+    Scalar solidDensity(const Element&, const SubControlVolume&) const
+    { return this->spatialParams().solidDensity(); }
+
+    // the Newmark scheme is used for time integration and this
+    // computes the acceleration at the current time step for us
+    auto acceleration(const Element& element,
+                      const SubControlVolume& scv,
+                      const Scalar dt,
+                      const PrimaryVariables& d) const
+    {
+        const auto dofIndex = scv.dofIndex();
+        return newmark_->acceleration(dofIndex, dt, d);
+    }
+
+private:
+    static constexpr Scalar eps_ = 1e-7;
+    Scalar gravity_;
+
+    std::shared_ptr<const Experimental::NewmarkBeta<Scalar, SolutionVector>> newmark_;
+};
+
+} // end namespace Dumux
+
+#endif

test/geomechanics/hyperelastic_dynamic/properties.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+#ifndef DUMUX_DYNAMIC_HYPERELASTICITY_TEST_PROPERTIES_HH
+#define DUMUX_DYNAMIC_HYPERELASTICITY_TEST_PROPERTIES_HH
+
+#include <type_traits>
+
+#include <dune/alugrid/grid.hh>
+#include <dumux/discretization/box.hh>
+#include <dumux/geomechanics/hyperelastic/model.hh>
+
+#include "problem.hh"
+
+namespace Dumux::Properties {
+
+namespace TTag {
+struct DynamicHyperelasticityTest
+{
+    using InheritsFrom = std::tuple<Hyperelastic, BoxModel>;
+    using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
+    using EnableGridVolumeVariablesCache = std::true_type;
+    using EnableGridFluxVariablesCache = std::true_type;
+    using EnableGridGeometryCache = std::true_type;
+};
+} // end namespace TTag
+
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DynamicHyperelasticityTest>
+{ using type = DynamicHyperelasticityProblem<TypeTag>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DynamicHyperelasticityTest>
+{
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = DefaultDynamicHyperelasticSpatialParams<GridGeometry, Scalar>;
+};
+
+
+} // end namespace Dumux::Properties
+
+#endif

test/geomechanics/hyperelastic_dynamic/solid.geo

+//+
+SetFactory("OpenCASCADE");
+Rectangle(1) = {0.2, 0.19, 0, 0.4, 0.02, 0};
+//+
+Disk(2) = {0.2, 0.2, 0, 0.05, 0.05};
+//+
+BooleanDifference{ Surface{1}; Delete; }{ Surface{2}; Delete; }