el2pmodel.hh 31.7 KB
Newer Older
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
// -*- 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 <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/

/*!
* \file
*
* \brief Adaption of the fully implicit scheme to the two-phase linear elasticity model.
*/

#ifndef DUMUX_ELASTIC2P_MODEL_HH
#define DUMUX_ELASTIC2P_MODEL_HH

#include "el2pproperties.hh"
#include <dumux/common/eigenvalues.hh>
#include <dune/pdelab/gridfunctionspace/interpolate.hh>

namespace Dumux {
/*!
 * \ingroup ElTwoPBoxModel
 * \brief Adaption of the fully implicit scheme to the two-phase linear elasticity model.
 *
 * This model implements a two-phase flow of compressible immiscible fluids \f$\alpha \in \{ w, n \}\f$.
 * The deformation of the solid matrix is described with a quasi-stationary momentum balance equation.
 * The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941).
 * The total stress acting on a rock is partially supported by the rock matrix and partially supported
 * by the pore fluid. The effective stress represents the share of the total stress which is supported
 * by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
 *
 * As an equation for the conservation of momentum within the fluid phases the standard multiphase Darcy's approach is used:
 \f[
 v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
 \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
 \f]
 *
 * Gravity can be enabled or disabled via the property system.
 * By inserting this into the continuity equation, one gets
\f[
 \frac{\partial \phi_{eff} \varrho_\alpha S_\alpha}{\partial t}
 - \text{div} \left\{ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}
 \mbox{\bf K_{eff}} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
 - \phi_{eff} \varrho_\alpha S_\alpha \frac{\partial \mathbf{u}}{\partial t}
 \right\} - q_\alpha = 0 \;,
 \f]
 *
 *
 * A quasi-stationary momentum balance equation is solved for the changes with respect to the initial conditions (Darcis 2012), note
 * that this implementation assumes the soil mechanics sign convention (i.e. compressive stresses are negative):
 \f[
 \text{div}\left( \boldsymbol{\Delta \sigma'}- \Delta p_{eff} \boldsymbol{I} \right) + \Delta \varrho_b {\textbf g} = 0 \;,
 \f]
 * with the effective stress:
 \f[
  \boldsymbol{\sigma'} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \mathbf{I}.
 \f]
 *
 * and the strain tensor \f$\boldsymbol{\epsilon}\f$ as a function of the solid displacement gradient \f$\textbf{grad} \mathbf{u}\f$:
 \f[
  \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \mathbf{u} + \textbf{grad}^T \mathbf{u}).
 \f]
 *
 * Here, the rock mechanics sign convention is switch off which means compressive stresses are < 0 and tensile stresses are > 0.
 * The rock mechanics sign convention can be switched on for the vtk output via the property system.
 *
 * The effective porosity and the effective permeability are calculated as a function of the solid displacement:
 \f[
      \phi_{eff} = \frac{\phi_{init} + \text{div} \mathbf{u}}{1 + \text{div} \mathbf{u}}
 \f]
 \f[
      K_{eff} = K_{init} \text{exp}\left( 22.2(\phi_{eff}/\phi_{init} -1 )\right)
 \f]
 * The mass balance equations are discretized using a vertex-centered finite volume (box)
 * or cell-centered finite volume scheme as spatial and the implicit Euler method as time discretization.
 * The momentum balance equations are discretized using a standard Galerkin Finite Element method as
 * spatial discretization scheme.
 *
 *
 * The primary variables are the wetting phase pressure \f$p_w\f$, the nonwetting phase saturation \f$S_n\f$ and the solid
 * displacement vector \f$\mathbf{u}\f$ (changes in solid displacement with respect to initial conditions).
 */

namespace Properties {
NEW_PROP_TAG(InitialDisplacement); //!< The initial displacement function
NEW_PROP_TAG(InitialPressSat); //!< The initial pressure and saturation function
}

template<class TypeTag>
class ElTwoPModel: public GET_PROP_TYPE(TypeTag, BaseModel)
{
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;

    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    enum {
        numEq = GET_PROP_VALUE(TypeTag, NumEq),
        nPhaseIdx = Indices::nPhaseIdx,
114
        wPhaseIdx = Indices::wPhaseIdx
115
116
117
118
119
120
121
122
123
    };

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
124
125
126
127
128
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
            typedef typename Element::Geometry::JacobianInverseTransposed JacobianInverseTransposed;
#else
            typedef typename Element::Geometry::Jacobian JacobianInverseTransposed;
#endif
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

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef Dune::FieldVector<Scalar, dim> DimVector;
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
    typedef Dune::PDELab::LocalFunctionSpace<GridFunctionSpace> LocalFunctionSpace;

public:

    /*!
     * \brief Write the current solution to a restart file.
     *
     * \param outStream The output stream of one vertex for the restart file
     * \param entity The Entity
     *
     * Due to the mixed discretization schemes which are combined via pdelab for this model
     * the solution vector has a different form than in the pure box models
     * it sorts the primary variables in the following way:
     * p_vertex0 S_vertex0 p_vertex1 S_vertex1 p_vertex2 ....p_vertexN S_vertexN
     * ux_vertex0 uy_vertex0 uz_vertex0 ux_vertex1 uy_vertex1 uz_vertex1 ...
     *
     * Therefore, the serializeEntity function has to be modified.
     */
    template <class Entity>
155
    void serializeEntity(std::ostream &outStream,
156
157
158
159
160
161
                         const Entity &entity)
    {
        // vertex index
        int dofIdx = this->dofMapper().map(entity);

        // write phase state
162
        if (!outStream.good()) {
163
164
165
166
167
168
169
            DUNE_THROW(Dune::IOError,
                       "Could not serialize vertex "
                       << dofIdx);
        }
        int numScv = this->gridView().size(dim);
        // get p and S entries for this vertex
        for (int eqIdx = 0; eqIdx < numEq-dim; ++eqIdx) {
170
            outStream << this->curSol()[dofIdx*(numEq-dim) + eqIdx][0]<<" ";
171
172
173
        }
        // get ux, uy, uz entries for this vertex
        for (int j = 0; j< dim; ++j)
174
            outStream << this->curSol()[numScv*(numEq-dim) + dofIdx*dim + j][0] <<" ";
175

176
        int vIdx = this->dofMapper().map(entity);
177
        if (!outStream.good())
178
            DUNE_THROW(Dune::IOError, "Could not serialize vertex " << vIdx);
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    }

    /*!
     * \brief Reads the current solution for a vertex from a restart
     *        file.
     *
     * \param inStream The input stream of one vertex from the restart file
     * \param entity The Entity
     *
     * Due to the mixed discretization schemes which are combined via pdelab for this model
     * the solution vector has a different form than in the pure box models
     * it sorts the primary variables in the following way:
     * p_vertex0 S_vertex0 p_vertex1 S_vertex1 p_vertex2 ....p_vertexN S_vertexN
     * ux_vertex0 uy_vertex0 uz_vertex0 ux_vertex1 uy_vertex1 uz_vertex1 ...
     *
     * Therefore, the deserializeEntity function has to be modified.
     */
    template<class Entity>
197
    void deserializeEntity(std::istream &inStream, const Entity &entity)
198
199
200
    {
        int dofIdx = this->dofMapper().map(entity);

201
        if (!inStream.good()){
202
203
204
205
206
207
208
                DUNE_THROW(Dune::IOError,
                           "Could not deserialize vertex "
                           << dofIdx);
        }
        int numScv = this->gridView().size(dim);
        for (int eqIdx = 0; eqIdx < numEq-dim; ++eqIdx) {
        // read p and S entries for this vertex
209
        inStream >> this->curSol()[dofIdx*(numEq-dim) + eqIdx][0];}
210
211
        for (int j = 0; j< dim; ++j){
            // read ux, uy, uz entries for this vertex
212
            inStream >> this->curSol()[numScv*(numEq-dim) + dofIdx*dim + j][0];}
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
    }


    /*!
     * \brief \copybrief ImplicitModel::addOutputVtkFields
     *
     * Specialization for the ElOnePTwoCBoxModel, add one-phase two-component
     * properties, solid displacement, stresses, effective properties and the
     * process rank to the VTK writer.
     */
    template<class MultiWriter>
    void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) {
        // check whether compressive stresses are defined to be positive
        // (rockMechanicsSignConvention_ == true) or negative
        rockMechanicsSignConvention_ =  GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, RockMechanicsSignConvention);

        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
        typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;

        // create the required scalar and vector fields
        unsigned numScv = this->gridView_().size(dim);
        unsigned numElements = this->gridView_().size(0);

        // create the required fields for vertex data
        ScalarField &pw = *writer.allocateManagedBuffer(numScv);
        ScalarField &pn = *writer.allocateManagedBuffer(numScv);
        ScalarField &pc = *writer.allocateManagedBuffer(numScv);
        ScalarField &sw = *writer.allocateManagedBuffer(numScv);
        ScalarField &sn = *writer.allocateManagedBuffer(numScv);
        VectorField &displacement = *writer.template allocateManagedBuffer<Scalar, dim>(numScv);
        ScalarField &rhoW = *writer.allocateManagedBuffer(numScv);
        ScalarField &rhoN = *writer.allocateManagedBuffer(numScv);
        ScalarField &Te = *writer.allocateManagedBuffer(numScv);

        // create the required fields for element data
        // effective stresses
        VectorField &deltaEffStressX = *writer.template allocateManagedBuffer<Scalar,
                dim>(numElements);
        VectorField &deltaEffStressY = *writer.template allocateManagedBuffer<Scalar,
                dim>(numElements);
        VectorField &deltaEffStressZ = *writer.template allocateManagedBuffer<Scalar,
                dim>(numElements);
        // total stresses
        VectorField &totalStressX = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        VectorField &totalStressY = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        VectorField &totalStressZ = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        // initial stresses
        VectorField &initStressX = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        VectorField &initStressY = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        VectorField &initStressZ = *writer.template allocateManagedBuffer<
                Scalar, dim>(numElements);
        // principal stresses
        ScalarField &principalStress1 = *writer.allocateManagedBuffer(
                numElements);
        ScalarField &principalStress2 = *writer.allocateManagedBuffer(
                numElements);
        ScalarField &principalStress3 = *writer.allocateManagedBuffer(
                numElements);


        ScalarField &effKx = *writer.allocateManagedBuffer(numElements);
        ScalarField &effPorosity = *writer.allocateManagedBuffer(numElements);
        ScalarField &effectivePressure = *writer.allocateManagedBuffer(numElements);
        ScalarField &deltaEffPressure = *writer.allocateManagedBuffer(numElements);


        ScalarField &Pcrtens = *writer.allocateManagedBuffer(numElements);
        ScalarField &Pcrshe = *writer.allocateManagedBuffer(numElements);

        // initialize cell stresses, cell-wise hydraulic parameters and cell pressure with zero


290
291
        for(int eIdx = 0; eIdx < numElements; ++eIdx){
            deltaEffStressX[eIdx] = Scalar(0.0);
292
            if (dim >= 2)
293
                deltaEffStressY[eIdx] = Scalar(0.0);
294
            if (dim >= 3)
295
                deltaEffStressZ[eIdx] = Scalar(0.0);
296

297
            totalStressX[eIdx] = Scalar(0.0);
298
            if (dim >= 2)
299
                totalStressY[eIdx] = Scalar(0.0);
300
            if (dim >= 3)
301
                totalStressZ[eIdx] = Scalar(0.0);
302

303
            initStressX[eIdx] = Scalar(0.0);
304
            if (dim >= 2)
305
                initStressY[eIdx] = Scalar(0.0);
306
            if (dim >= 3)
307
                initStressZ[eIdx] = Scalar(0.0);
308

309
            principalStress1[eIdx] = Scalar(0.0);
310
            if (dim >= 2)
311
                principalStress2[eIdx] = Scalar(0.0);
312
            if (dim >= 3)
313
                principalStress3[eIdx] = Scalar(0.0);
314

315
316
317
318
            effPorosity[eIdx] = Scalar(0.0);
            effKx[eIdx] = Scalar(0.0);
            effectivePressure[eIdx] = Scalar(0.0);
            deltaEffPressure[eIdx] = Scalar(0.0);
319

320
321
            Pcrtens[eIdx] = Scalar(0.0);
            Pcrshe[eIdx] = Scalar(0.0);
322
323
324
325
326
327
328
329
330
        }

        ScalarField &rank = *writer.allocateManagedBuffer(numElements);


        FVElementGeometry fvGeometry;
        ElementVolumeVariables elemVolVars;

        // initialize start and end of element iterator
331
332
        ElementIterator eIt = this->gridView_().template begin<0>();
        ElementIterator eEndit = this->gridView_().template end<0>();
333
        // loop over all elements (cells)
334
        for (; eIt != eEndit; ++eIt) {
335
336
337
338
339

            // get FE function spaces to calculate gradients (gradient data of momentum balance
            // equation is not stored in fluxvars since it is not evaluated at box integration point)
            // copy the values of the sol vector to the localFunctionSpace values of the current element
            LocalFunctionSpace localFunctionSpace(this->problem_().model().jacobianAssembler().gridFunctionSpace());
340
            localFunctionSpace.bind(*eIt);
341
342
343
344
345
346
347
348
349
350
351
352
            std::vector<Scalar> values;
            localFunctionSpace.vread(sol, values);

            // local function space for solid displacement
            typedef typename LocalFunctionSpace::template Child<1>::Type DisplacementLFS;
            const DisplacementLFS& displacementLFS =localFunctionSpace.template child<1>();
            const unsigned int dispSize = displacementLFS.child(0).size();
            typedef typename DisplacementLFS::template Child<0>::Type ScalarDispLFS;
            // further types required for gradient calculations
            typedef typename ScalarDispLFS::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
            typedef typename ScalarDispLFS::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType RF;

353
354
            unsigned int eIdx = this->problem_().model().elementMapper().map(*eIt);
            rank[eIdx] = this->gridView_().comm().rank();
355

356
357
            fvGeometry.update(this->gridView_(), *eIt);
            elemVolVars.update(this->problem_(), *eIt, fvGeometry, false);
358
359

            // loop over all local vertices of the cell
360
            int numScv = eIt->template count<dim>();
361
362
363

            for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
            {
364
                unsigned int globalIdx = this->dofMapper().map(*eIt, scvIdx, dim);
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392


                Te[globalIdx] = elemVolVars[scvIdx].temperature();
                pw[globalIdx] = elemVolVars[scvIdx].pressure(wPhaseIdx);
                pn[globalIdx] = elemVolVars[scvIdx].pressure(nPhaseIdx);
                pc[globalIdx] = elemVolVars[scvIdx].capillaryPressure();
                sw[globalIdx] = elemVolVars[scvIdx].saturation(wPhaseIdx);
                sn[globalIdx] = elemVolVars[scvIdx].saturation(nPhaseIdx);
                rhoW[globalIdx] = elemVolVars[scvIdx].density(wPhaseIdx);
                rhoN[globalIdx] = elemVolVars[scvIdx].density(nPhaseIdx);
                // the following lines are correct for rock mechanics sign convention
                // but lead to a very counter-intuitive output therefore, they are commented.
                // in case of rock mechanics sign convention solid displacement is
                // defined to be negative if it points in positive coordinate direction
//                if(rockMechanicsSignConvention_){
//                    DimVector tmpDispl;
//                    tmpDispl = Scalar(0);
//                    tmpDispl -= elemVolVars[scvIdx].displacement();
//                    displacement[globalIdx] = tmpDispl;
//                    }
//
//                else
                    displacement[globalIdx] = elemVolVars[scvIdx].displacement();

                double Keff;
                double exponent;
                exponent = 22.2    * (elemVolVars[scvIdx].effPorosity
                            / elemVolVars[scvIdx].porosity() - 1);
393
                Keff =    this->problem_().spatialParams().intrinsicPermeability(    *eIt, fvGeometry, scvIdx)[0][0];
394
                Keff *= exp(exponent);
395
396
                effKx[eIdx] += Keff/ numScv;
                effectivePressure[eIdx] += (pn[globalIdx] * sn[globalIdx]
397
398
                                            + pw[globalIdx] * sw[globalIdx])
                                            / numScv;
399
                effPorosity[eIdx] +=elemVolVars[scvIdx].effPorosity / numScv;
400
401
            };

402
403
            const GlobalPosition& cellCenter = eIt->geometry().center();
            const GlobalPosition& cellCenterLocal = eIt->geometry().local(cellCenter);
404

405
            deltaEffPressure[eIdx] = effectivePressure[eIdx] + this->problem().pInit(cellCenter, cellCenterLocal, *eIt);
406
407
408
409
410
411
            // determin changes in effective stress from current solution
            // evaluate gradient of displacement shape functions
            std::vector<JacobianType_V> vRefShapeGradient(dispSize);
            displacementLFS.child(0).finiteElement().localBasis().evaluateJacobian(cellCenterLocal, vRefShapeGradient);

            // get jacobian to transform the gradient to physical element
412
            const JacobianInverseTransposed jacInvT = eIt->geometry().jacobianInverseTransposed(cellCenterLocal);
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
            std::vector < Dune::FieldVector<RF, dim> > vShapeGradient(dispSize);
            for (size_t i = 0; i < dispSize; i++) {
                vShapeGradient[i] = 0.0;
                jacInvT.umv(vRefShapeGradient[i][0], vShapeGradient[i]);
            }
            // calculate gradient of current displacement
            typedef Dune::FieldMatrix<RF, dim, dim> DimMatrix;
            DimMatrix uGradient(0.0);
            for (int coordDir = 0; coordDir < dim; ++coordDir) {
                const ScalarDispLFS & scalarDispLFS = displacementLFS.child(coordDir);

                for (size_t i = 0; i < scalarDispLFS.size(); i++)
                    uGradient[coordDir].axpy(values[scalarDispLFS.localIndex(i)],vShapeGradient[i]);
            }

428
            const Dune::FieldVector<Scalar, 2> lameParams =    this->problem_().spatialParams().lameParams(*eIt,fvGeometry, 0);
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
            const Scalar lambda = lameParams[0];
            const Scalar mu = lameParams[1];

            // calculate strain tensor
            Dune::FieldMatrix<RF, dim, dim> epsilon;
            for (int i = 0; i < dim; ++i)
                for (int j = 0; j < dim; ++j)
                    epsilon[i][j] = 0.5 * (uGradient[i][j] + uGradient[j][i]);

            RF traceEpsilon = 0;
            for (int i = 0; i < dim; ++i)
                traceEpsilon += epsilon[i][i];

            // calculate effective stress tensor
            Dune::FieldMatrix<RF, dim, dim> sigma(0.0);
            for (int i = 0; i < dim; ++i) {
                sigma[i][i] = lambda * traceEpsilon;
                for (int j = 0; j < dim; ++j)
                    sigma[i][j] += 2.0 * mu * epsilon[i][j];
            }

            // in case of rock mechanics sign convention compressive stresses
            // are defined to be positive
            if(rockMechanicsSignConvention_){
453
                deltaEffStressX[eIdx] -= sigma[0];
454
                if (dim >= 2) {
455
                    deltaEffStressY[eIdx] -= sigma[1];
456
457
                }
                if (dim >= 3) {
458
                    deltaEffStressZ[eIdx] -= sigma[2];
459
460
461
                }
            }
            else{
462
                deltaEffStressX[eIdx] = sigma[0];
463
                if (dim >= 2) {
464
                    deltaEffStressY[eIdx] = sigma[1];
465
466
                }
                if (dim >= 3) {
467
                    deltaEffStressZ[eIdx] = sigma[2];
468
469
470
471
472
473
                }
            }

            // retrieve prescribed initial stresses from problem file
            DimVector tmpInitStress = this->problem_().initialStress(cellCenter, 0);
            if(rockMechanicsSignConvention_){
474
                initStressX[eIdx][0] = tmpInitStress[0];
475
                if (dim >= 2) {
476
                    initStressY[eIdx][1] = tmpInitStress[1];
477
478
                    }
                if (dim >= 3) {
479
                    initStressZ[eIdx][2] = tmpInitStress[2];
480
481
482
                }
            }
            else{
483
                initStressX[eIdx][0] -= tmpInitStress[0];
484
                if (dim >= 2) {
485
                    initStressY[eIdx][1] -= tmpInitStress[1];
486
487
                    }
                if (dim >= 3) {
488
                    initStressZ[eIdx][2] -= tmpInitStress[2];
489
490
491
492
493
494
495
                }
            }

            // calculate total stresses
            // in case of rock mechanics sign convention compressive stresses
            // are defined to be positive and total stress is calculated by adding the pore pressure
            if(rockMechanicsSignConvention_){
496
497
498
                totalStressX[eIdx][0] = initStressX[eIdx][0] + deltaEffStressX[eIdx][0]    + deltaEffPressure[eIdx];
                totalStressX[eIdx][1] = initStressX[eIdx][1] + deltaEffStressX[eIdx][1];
                totalStressX[eIdx][2] = initStressX[eIdx][2] + deltaEffStressX[eIdx][2];
499
                if (dim >= 2) {
500
501
502
                    totalStressY[eIdx][0] = initStressY[eIdx][0] + deltaEffStressY[eIdx][0];
                    totalStressY[eIdx][1] = initStressY[eIdx][1] + deltaEffStressY[eIdx][1]    + deltaEffPressure[eIdx];
                    totalStressY[eIdx][2] = initStressY[eIdx][2] + deltaEffStressY[eIdx][2];
503
504
                }
                if (dim >= 3) {
505
506
507
                    totalStressZ[eIdx][0] = initStressZ[eIdx][0] + deltaEffStressZ[eIdx][0];
                    totalStressZ[eIdx][1] = initStressZ[eIdx][1] + deltaEffStressZ[eIdx][1];
                    totalStressZ[eIdx][2] = initStressZ[eIdx][2] + deltaEffStressZ[eIdx][2]    + deltaEffPressure[eIdx];
508
509
510
                }
            }
            else{
511
512
513
                totalStressX[eIdx][0] = initStressX[eIdx][0] + deltaEffStressX[eIdx][0]    - deltaEffPressure[eIdx];
                totalStressX[eIdx][1] = initStressX[eIdx][1] + deltaEffStressX[eIdx][1];
                totalStressX[eIdx][2] = initStressX[eIdx][2] + deltaEffStressX[eIdx][2];
514
                if (dim >= 2) {
515
516
517
                    totalStressY[eIdx][0] = initStressY[eIdx][0] + deltaEffStressY[eIdx][0];
                    totalStressY[eIdx][1] = initStressY[eIdx][1] + deltaEffStressY[eIdx][1]    - deltaEffPressure[eIdx];
                    totalStressY[eIdx][2] = initStressY[eIdx][2] + deltaEffStressY[eIdx][2];
518
519
                }
                if (dim >= 3) {
520
521
522
                    totalStressZ[eIdx][0] = initStressZ[eIdx][0] + deltaEffStressZ[eIdx][0];
                    totalStressZ[eIdx][1] = initStressZ[eIdx][1] + deltaEffStressZ[eIdx][1];
                    totalStressZ[eIdx][2] = initStressZ[eIdx][2] + deltaEffStressZ[eIdx][2]    - deltaEffPressure[eIdx];
523
524
525
526
527
528
529
530
531
                }
            }
        }

        // calculate principal stresses i.e. the eigenvalues of the total stress tensor
        Scalar a1, a2, a3;
        DimMatrix totalStress;
        DimVector eigenValues;

532
        for (unsigned int eIdx = 0; eIdx < numElements; eIdx++)
533
534
535
536
        {
            eigenValues = Scalar(0);
            totalStress = Scalar(0);

537
            totalStress[0] = totalStressX[eIdx];
538
            if (dim >= 2)
539
                totalStress[1] = totalStressY[eIdx];
540
            if (dim >= 3)
541
                totalStress[2] = totalStressZ[eIdx];
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557

            calculateEigenValues<dim>(eigenValues, totalStress);


            for (int i = 0; i < dim; i++)
                {
                    if (isnan(eigenValues[i]))
                        eigenValues[i] = 0.0;
                }

            // sort principal stresses: principalStress1 >= principalStress2 >= principalStress3
            if (dim == 2) {
                a1 = eigenValues[0];
                a2 = eigenValues[1];

                if (a1 >= a2) {
558
559
                    principalStress1[eIdx] = a1;
                    principalStress2[eIdx] = a2;
560
                } else {
561
562
                    principalStress1[eIdx] = a2;
                    principalStress2[eIdx] = a1;
563
564
565
566
567
568
569
570
571
572
                }
            }

            if (dim == 3) {
                a1 = eigenValues[0];
                a2 = eigenValues[1];
                a3 = eigenValues[2];

                if (a1 >= a2) {
                    if (a1 >= a3) {
573
                        principalStress1[eIdx] = a1;
574
                        if (a2 >= a3) {
575
576
                            principalStress2[eIdx] = a2;
                            principalStress3[eIdx] = a3;
577
578
579
                        }
                        else //a3 > a2
                        {
580
581
                            principalStress2[eIdx] = a3;
                            principalStress3[eIdx] = a2;
582
583
584
585
                        }
                    }
                    else // a3 > a1
                    {
586
587
588
                        principalStress1[eIdx] = a3;
                        principalStress2[eIdx] = a1;
                        principalStress3[eIdx] = a2;
589
590
591
592
                    }
                } else // a2>a1
                {
                    if (a2 >= a3) {
593
                        principalStress1[eIdx] = a2;
594
                        if (a1 >= a3) {
595
596
                            principalStress2[eIdx] = a1;
                            principalStress3[eIdx] = a3;
597
598
599
                        }
                        else //a3>a1
                        {
600
601
                            principalStress2[eIdx] = a3;
                            principalStress3[eIdx] = a1;
602
603
604
605
                        }
                    }
                    else //a3>a2
                    {
606
607
608
                        principalStress1[eIdx] = a3;
                        principalStress2[eIdx] = a2;
                        principalStress3[eIdx] = a1;
609
610
611
612
613
                    }
                }
            }
            Scalar taum  = 0.0;
            Scalar sigmam = 0.0;
614
            Scalar Peff = effectivePressure[eIdx];
615
616
617

            Scalar theta = M_PI / 6;
            Scalar S0 = 0.0;
618
619
            taum = (principalStress1[eIdx] - principalStress3[eIdx]) / 2;
            sigmam = (principalStress1[eIdx] + principalStress3[eIdx]) / 2;
620
621
622
            Scalar Psc = -fabs(taum) / sin(theta) + S0 * cos(theta) / sin(theta)
                    + sigmam;
            // Pressure margins according to J. Rutqvist et al. / International Journal of Rock Mecahnics & Mining Sciences 45 (2008), 132-143
623
624
            Pcrtens[eIdx] = Peff - principalStress3[eIdx];
            Pcrshe[eIdx] = Peff - Psc;
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687

        }

        writer.attachVertexData(Te, "T");
        writer.attachVertexData(pw, "pW");
        writer.attachVertexData(pn, "pN");
        writer.attachVertexData(pc, "pC");
        writer.attachVertexData(sw, "SW");
        writer.attachVertexData(sn, "SN");
        writer.attachVertexData(rhoW, "rhoW");
        writer.attachVertexData(rhoN, "rhoN");
        writer.attachVertexData(displacement, "u", dim);

        writer.attachCellData(deltaEffStressX, "effective stress changes X", dim);
        if (dim >= 2)
            writer.attachCellData(deltaEffStressY, "effective stress changes Y",    dim);
        if (dim >= 3)
            writer.attachCellData(deltaEffStressZ, "effective stress changes Z",    dim);

        writer.attachCellData(principalStress1, "principal stress 1");
        if (dim >= 2)
            writer.attachCellData(principalStress2, "principal stress 2");
        if (dim >= 3)
            writer.attachCellData(principalStress3, "principal stress 3");

        writer.attachCellData(totalStressX, "total stresses X", dim);
        if (dim >= 2)
            writer.attachCellData(totalStressY, "total stresses Y", dim);
        if (dim >= 3)
            writer.attachCellData(totalStressZ, "total stresses Z", dim);

        writer.attachCellData(initStressX, "initial stresses X", dim);
        if (dim >= 2)
            writer.attachCellData(initStressY, "initial stresses Y", dim);
        if (dim >= 3)
            writer.attachCellData(initStressZ, "initial stresses Z", dim);

        writer.attachCellData(deltaEffPressure, "delta pEff");
         writer.attachCellData(effectivePressure, "effectivePressure");
        writer.attachCellData(Pcrtens, "Pcr_tensile");
        writer.attachCellData(Pcrshe, "Pcr_shear");
        writer.attachCellData(effKx, "effective Kxx");
        writer.attachCellData(effPorosity, "effective Porosity");


    }

    /*!
     * \brief Applies the initial solution for all vertices of the grid.
     */
    void applyInitialSolution_() {
        typedef typename GET_PROP_TYPE(TypeTag, PTAG(InitialPressSat)) InitialPressSat;
        InitialPressSat initialPressSat(this->problem_().gridView());
        std::cout << "el2pmodel calls: initialPressSat" << std::endl;
        initialPressSat.setPressure(this->problem_().pInit());

        typedef typename GET_PROP_TYPE(TypeTag, PTAG(InitialDisplacement)) InitialDisplacement;
        InitialDisplacement initialDisplacement(this->problem_().gridView());

        typedef Dune::PDELab::CompositeGridFunction<InitialPressSat,
                InitialDisplacement> InitialSolution;
        InitialSolution initialSolution(initialPressSat, initialDisplacement);

688
689
690
691
        int numDofs = this->jacobianAssembler().gridFunctionSpace().size();
        this->curSol().resize(numDofs);
        this->prevSol().resize(numDofs);
        std::cout << "numDofs = " << numDofs << std::endl;
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709

        Dune::PDELab::interpolate(initialSolution,
                this->jacobianAssembler().gridFunctionSpace(), this->curSol());
        Dune::PDELab::interpolate(initialSolution,
                this->jacobianAssembler().gridFunctionSpace(), this->prevSol());
    }

    const Problem& problem() const {
        return this->problem_();
    }

private:
    bool rockMechanicsSignConvention_;

};
}
#include "el2ppropertydefaults.hh"
#endif