File: fvbasediscretizationfemadapt.hh

package info (click to toggle)
opm-simulators 2024.10%2Bds-6
  • links: PTS, VCS
  • area: main
  • in suites: trixie
  • size: 19,416 kB
  • sloc: cpp: 165,337; sh: 1,285; lisp: 1,108; python: 355; makefile: 24; awk: 10
file content (174 lines) | stat: -rw-r--r-- 6,921 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM 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.

  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
/*!
 * \file
 *
 * \copydoc Opm::FvBaseDiscretization
 */
#ifndef EWOMS_FV_BASE_DISCRETIZATION_FEMADAPT_HH
#define EWOMS_FV_BASE_DISCRETIZATION_FEMADAPT_HH

#include <opm/models/discretization/common/fvbasediscretization.hh>

#include <dune/fem/space/common/adaptationmanager.hh>
#include <dune/fem/space/common/restrictprolongtuple.hh>
#include <dune/fem/function/blockvectorfunction.hh>
#include <dune/fem/misc/capabilities.hh>

namespace Opm {

template<class TypeTag>
class FvBaseDiscretizationFemAdapt;

namespace Properties {

template<class TypeTag>
struct BaseDiscretizationType<TypeTag,TTag::FvBaseDiscretization> {
    using type = FvBaseDiscretizationFemAdapt<TypeTag>;
};

template<class TypeTag>
struct DiscreteFunction<TypeTag, TTag::FvBaseDiscretization> {
    using DiscreteFunctionSpace  = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
    using PrimaryVariables  = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using type = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;
};

} // namespace Properties

/*!
 * \ingroup FiniteVolumeDiscretizations
 *
 * \brief The base class for the finite volume discretization schemes.
 */

template <class TypeTag>
class FvBaseDiscretizationFemAdapt : public FvBaseDiscretization<TypeTag>
{
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using ParentType = FvBaseDiscretization<TypeTag>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;

    static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();

    using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
    // discrete function storing solution data
    using DiscreteFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;

    // problem restriction and prolongation operator for adaptation
    using ProblemRestrictProlongOperator = typename Problem::RestrictProlongOperator;

    // discrete function restriction and prolongation operator for adaptation
    using DiscreteFunctionRestrictProlong = Dune::Fem::RestrictProlongDefault<DiscreteFunction>;
    using RestrictProlong
        = Dune::Fem::RestrictProlongTuple<DiscreteFunctionRestrictProlong, ProblemRestrictProlongOperator>;

    // adaptation classes
    using AdaptationManager = Dune::Fem::AdaptationManager<Grid, RestrictProlong>;

public:
    template<class Serializer>
    struct SerializeHelper {
        template<class SolutionType>
        static void serializeOp(Serializer& serializer,
                                SolutionType& solution)
        {
            for (auto& sol : solution) {
                serializer(sol->blockVector());
            }
        }
    };

    FvBaseDiscretizationFemAdapt(Simulator& simulator)
        : ParentType(simulator)
        , space_(simulator.vanguard().gridPart())
    {
        if (this->enableGridAdaptation_ && !Dune::Fem::Capabilities::isLocallyAdaptive<Grid>::v) {
            throw std::invalid_argument("Grid adaptation enabled, but chosen Grid is not capable"
                                        " of adaptivity");
        }

        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", space_);
        }
    }

    void adaptGrid()
    {
        // adapt the grid if enabled and if all dependencies are available
        // adaptation is only done if markForGridAdaptation returns true
        if (this->enableGridAdaptation_) {
            // check if problem allows for adaptation and cells were marked
            if (this->simulator_.problem().markForGridAdaptation()) {
                // adapt the grid and load balance if necessary
                adaptationManager().adapt();

                // if the grid has potentially changed, we need to re-create the
                // supporting data structures.
                this->elementMapper_.update(this->gridView_);
                this->vertexMapper_.update(this->gridView_);
                this->resetLinearizer();
                // this is a bit hacky because it supposes that Problem::finishInit()
                // works fine multiple times in a row.
                //
                // TODO: move this to Problem::gridChanged()
                this->finishInit();

                // notify the problem that the grid has changed
                //
                // TODO: come up with a mechanism to access the unadapted data structures
                // outside of the problem (i.e., grid, mappers, solutions)
                this->simulator_.problem().gridChanged();

                // notify the modules for visualization output
                auto outIt = this->outputModules_.begin();
                auto outEndIt = this->outputModules_.end();
                for (; outIt != outEndIt; ++outIt)
                    (*outIt)->allocBuffers();
            }
        }
    }

    AdaptationManager& adaptationManager()
    {
        if (!adaptationManager_) {
            // create adaptation objects here, because when doing so in constructor
            // problem is not yet intialized, aka seg fault
            restrictProlong_ = std::make_unique<RestrictProlong>(DiscreteFunctionRestrictProlong(*(this->solution_[/*timeIdx=*/0])),
                                                                 this->simulator_.problem().restrictProlongOperator());
            adaptationManager_ = std::make_unique<AdaptationManager>(this->simulator_.vanguard().grid(), *restrictProlong_);
        }
        return *adaptationManager_;
    }

private:
    DiscreteFunctionSpace space_;
    std::unique_ptr<RestrictProlong> restrictProlong_;
    std::unique_ptr<AdaptationManager> adaptationManager_;
};

} // namespace Opm

#endif