3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
fvtransportmultiphysics.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_FVTRANSPORT2P2C_MULTIPHYSICS_HH
25#define DUMUX_FVTRANSPORT2P2C_MULTIPHYSICS_HH
26
29
30namespace Dumux {
54template<class TypeTag>
56{
60
62
64
66
67 enum
68 {
69 dim = GridView::dimension
70 };
71 enum
72 {
73 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
74 wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
75 };
76
77 using PhaseVector = Dune::FieldVector<Scalar, 2>;
79
81 Problem& problem()
82 { return this->problem_; }
83
84 using LocalTimesteppingData = typename FVTransport2P2C<TypeTag>::LocalTimesteppingData;
85
86public:
87 virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet = false);
88
90
93 FVTransport2P2CMultiPhysics(Problem& problem) : FVTransport2P2C<TypeTag>(problem)
94 {}
95
97 {}
98};
99
118template<class TypeTag>
119void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet)
120{
121 // initialize dt very large
122 dt = 1E100;
123
124 unsigned int size = problem().gridView().size(0);
125 if (this->enableLocalTimeStepping())
126 {
127 if (this->timeStepData_.size() != size)
128 this->timeStepData_.resize(size);
129 }
130 // store if we do update Estimate for flux functions
131 this->impet_ = impet;
132 this->averagedFaces_ = 0.;
133
134 // resize update vector and set to zero
135 updateVec.resize(getPropValue<TypeTag, Properties::NumComponents>());
136 updateVec[wCompIdx].resize(problem().gridView().size(0));
137 updateVec[nCompIdx].resize(problem().gridView().size(0));
138 updateVec[wCompIdx] = 0;
139 updateVec[nCompIdx] = 0;
140
141 // Cell which restricts time step size
142 int restrictingCell = -1;
143
144 PhaseVector entries(0.), timestepFlux(0.);
145 // compute update vector
146 for (const auto& element : elements(problem().gridView()))
147 {
148 // get cell infos
149 int globalIdxI = problem().variables().index(element);
150 CellData& cellDataI = problem().variables().cellData(globalIdxI);
151
152 if (impet || cellDataI.subdomain() == 2) // estimate only necessary in subdomain
153 {
154 // some variables for time step calculation
155 double sumfactorin = 0;
156 double sumfactorout = 0;
157
158 // run through all intersections with neighbors and boundary
159 for (const auto& intersection : intersections(problem().gridView(), element))
160 {
161 int indexInInside = intersection.indexInInside();
162
163 /****** interior face *****************/
164 if (intersection.neighbor())
165 this->getFlux(entries, timestepFlux, intersection, cellDataI);
166
167 /****** Boundary Face *****************/
168 if (intersection.boundary())
169 this->getFluxOnBoundary(entries, timestepFlux, intersection, cellDataI);
170
171 if (this->enableLocalTimeStepping())
172 {
173 LocalTimesteppingData& localData = this->timeStepData_[globalIdxI];
174
175 if (localData.faceTargetDt[indexInInside] < this->accumulatedDt_ + this->dtThreshold_)
176 {
177 localData.faceFluxes[indexInInside] = entries;
178 }
179 }
180 else
181 {
182 // add to update vector
183 updateVec[wCompIdx][globalIdxI] += entries[wCompIdx];
184 updateVec[nCompIdx][globalIdxI] += entries[nCompIdx];
185 }
186
187 // for time step calculation
188 sumfactorin += timestepFlux[0];
189 sumfactorout += timestepFlux[1];
190
191 }// end all intersections
192
193 if (this->enableLocalTimeStepping())
194 {
195 LocalTimesteppingData& localData = this->timeStepData_[globalIdxI];
196 for (int i=0; i < 2*dim; i++)
197 {
198 updateVec[wCompIdx][globalIdxI] += localData.faceFluxes[i][wCompIdx];
199 updateVec[nCompIdx][globalIdxI] += localData.faceFluxes[i][nCompIdx];
200 }
201 }
202
203 /*********** Handle source term ***************/
204 PrimaryVariables q(NAN);
205 problem().source(q, element);
206 updateVec[wCompIdx][globalIdxI] += q[Indices::contiWEqIdx];
207 updateVec[nCompIdx][globalIdxI] += q[Indices::contiNEqIdx];
208
209 // account for porosity in fluxes for time-step
210 using std::max;
211 sumfactorin = max(sumfactorin,sumfactorout)
212 / problem().spatialParams().porosity(element);
213
214 //calculate time step
215 if (this->enableLocalTimeStepping())
216 {
217 this->timeStepData_[globalIdxI].dt = 1./sumfactorin;
218 if ( 1./sumfactorin < dt)
219 {
220 dt = 1./sumfactorin;
221 restrictingCell= globalIdxI;
222 }
223 }
224 else
225 {
226 if ( 1./sumfactorin < dt)
227 {
228 dt = 1./sumfactorin;
229 restrictingCell= globalIdxI;
230 }
231 }
232 }
233 } // end grid traversal
234
235#if HAVE_MPI
236 // communicate updated values
238 using ElementMapper = typename SolutionTypes::ElementMapper;
240 for (int i = 0; i < updateVec.size(); i++)
241 {
242 DataHandle dataHandle(problem().variables().elementMapper(), updateVec[i]);
243 problem().gridView().template communicate<DataHandle>(dataHandle,
244 Dune::InteriorBorder_All_Interface,
245 Dune::ForwardCommunication);
246 }
247 dt = problem().gridView().comm().min(dt);
248#endif
249
250 if(impet)
251 {
252 Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell <<
253 " leads to dt = "<<dt * getParam<Scalar>("Impet.CFLFactor")<< std::endl;
254 if(this->averagedFaces_ != 0)
255 Dune::dinfo << " Averageing done for " << this->averagedFaces_ << " faces. "<< std::endl;
256 }
257 return;
258}
259} // end namespace Dumux
260#endif
Contains a class to exchange entries of a vector.
Finite volume discretization of the component transport equation.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
A data handle class to exchange entries of a vector.
Definition: vectorcommdatahandle.hh:78
Compositional transport step in a Finite Volume discretization.
Definition: fvtransport.hh:62
Problem & problem_
Definition: fvtransport.hh:302
Definition: fvtransport.hh:114
Compositional transport step in a finite volume discretization.
Definition: fvtransportmultiphysics.hh:56
virtual ~FVTransport2P2CMultiPhysics()
Definition: fvtransportmultiphysics.hh:96
FVTransport2P2CMultiPhysics(Problem &problem)
Constructs a FVTransport2P2CMultiPhysics object.
Definition: fvtransportmultiphysics.hh:93
virtual void update(const Scalar t, Scalar &dt, TransportSolutionType &updateVec, bool impet=false)
Calculate the update vector and determine timestep size.
Definition: fvtransportmultiphysics.hh:119