plotoverline2d.hh 8.65 KB
Newer Older
Philipp Nuske's avatar
Philipp Nuske committed
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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
/*!
 * \file
 * \brief Plot variables over a line specified by two arguments.
 *          These output files are meant for visualization with another program (matlab, gnuplot...)
 *
 */
#ifndef DUMUX_PLOTOVERLINE_2D_HH
#define DUMUX_PLOTOVERLINE_2D_HH

#include <dumux/common/valgrind.hh>
#include <dumux/common/propertysystem.hh>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <iostream>
#include <fstream>

namespace Dumux
{
namespace Properties
{
    NEW_PROP_TAG(Scalar);
    NEW_PROP_TAG(Problem);
    NEW_PROP_TAG(GridView);
    NEW_PROP_TAG(DofMapper);
    NEW_PROP_TAG(ElementSolutionVector);
    NEW_PROP_TAG(SolutionVector);
    NEW_PROP_TAG(FVElementGeometry);
    NEW_PROP_TAG(TwoPIAIndices);
    NEW_PROP_TAG(NumEq);
    NEW_PROP_TAG(MaterialLaw);
    NEW_PROP_TAG(ElementVolumeVariables);
    NEW_PROP_TAG(AwnSurface);
    NEW_PROP_TAG(AwnSurfaceParams);
}

template<class TypeTag>
class PlotOverLine2D
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar)               Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Problem)              Problem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView)             GridView;
    typedef typename GridView::template Codim<0>::Iterator        ElementIterator;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry)    FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, DofMapper)            DofMapper;
    typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector)       SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)          MaterialLaw;
    typedef typename MaterialLaw::Params                          aterialLawParams;

    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem)          FluidSystem;

    enum {
        wPhaseIdx = FluidSystem::wPhaseIdx,
        nPhaseIdx = FluidSystem::nPhaseIdx,
        sPhaseIdx = FluidSystem::sPhaseIdx,
        wCompIdx  = FluidSystem::wCompIdx,
        nCompIdx  = FluidSystem::nCompIdx,

        // Grid and world dimension
        dim         = GridView::dimension,
        dimWorld    = GridView::dimensionworld,
    };

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

public:
    /*
     * \brief A function that writes results over a line (like paraview's plotOverline into a text file.)
     *
     *        The writer needs to be called in postTimeStep().
     *
     *        This function puts output variables (TemperaturePhase, Saturation, t, tIndex, ...) over space (1D, over a line) into a text file,
     *        so they can be read in by another program like matlab.
     *        The file can be found by the extension: dat
     */
    void write(const Problem & problem,
               const GlobalPosition & pointOne,
               const GlobalPosition & pointTwo,
               const std::string appendOutputName = "")
    {
        static_assert(dim==2, "this implements plot over Line: so far this works only for 2D");

        FVElementGeometry fvGeometry;
        ElementVolumeVariables elemVolVars;

        // so many vertices in the domain
        const unsigned int numGlobalVerts = problem.gridView().size(dim);

        // check whether a vertex was already visited by true / false vector
        bool isVisited[numGlobalVerts] ;

        // loop all vertices of the main: initialize
        for (int idx=0; idx < numGlobalVerts; idx++)
            isVisited[idx] = false ;

        // Looping over all elements of the domain
        ElementIterator eEndIt = problem.gridView().template end<0>();
        for (ElementIterator eIt = problem.gridView().template begin<0>() ; eIt not_eq eEndIt; ++eIt)
        {
            // updating the volume variables
            fvGeometry.update(problem.gridView(), *eIt);
            elemVolVars.update(problem, *eIt, fvGeometry, false);

            // number of scv
            const unsigned int numScv = fvGeometry.numScv;

            for (int scvIdx=0; scvIdx < numScv; ++scvIdx) {

                // find some global identification
                const unsigned int globalIdx = problem.vertexMapper().map(*eIt, scvIdx, dim);

                // only write out if the vertex was not already visited
                if (isVisited[globalIdx])
                    continue;

                isVisited[globalIdx] = true ;

                // Getting the spatial coordinate
                const GlobalPosition & globalPosCurrent
                    = fvGeometry.subContVol[scvIdx].global;

                std::ofstream dataFile;
                const unsigned int timeStepIndex = problem.timeManager().timeStepIndex() ;

                // filename of the output file
                std::string fileName = problem.name();

                // filename consists of problem name + some function argument + .dat
                // this way several plot overlines can be written from one simulation
                fileName += appendOutputName;
                fileName +=".dat";

                // Writing a header into the output
                if (timeStepIndex == 0)
                {
                    dataFile.open(fileName.c_str());
                    dataFile << "# This is a DuMuX output file for further processing with the preferred graphics program of your choice. \n";

                    dataFile << "# This output file was written from "<< __FILE__ << ", line " <<__LINE__ << "\n";
                    dataFile << "# This output file was generated from code compiled at " << __TIME__ <<", "<< __DATE__<< "\n";
                    dataFile << "\n";
                    dataFile << "# Header\n";
                    dataFile << "#timestep\t time\t\t \t\t x \t\t y  \t\tSw \t\t\t Tw\t\t Tn\t Ts \t xH2On \t xH2OnEquil \t xN2w \txN2wEquil " << std::endl;
                    dataFile.close();
                }

                // write output if the current location is between two specified points
                if (isBetween(globalPosCurrent, pointOne, pointTwo))
                {
                    const Scalar time         = problem.timeManager().time();
                    const Scalar saturationW  = elemVolVars[scvIdx].fluidState().saturation(wPhaseIdx);
                    const Scalar Tw           = elemVolVars[scvIdx].fluidState().temperature(wPhaseIdx);
                    const Scalar Tn           = elemVolVars[scvIdx].fluidState().temperature(nPhaseIdx);
                    const Scalar Ts           = elemVolVars[scvIdx].fluidState().temperature(sPhaseIdx);
                    const Scalar xH2On        = elemVolVars[scvIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
                    const Scalar xH2OnEquil   = elemVolVars[scvIdx].xEquil(nPhaseIdx, wCompIdx);
                    const Scalar xN2w         = elemVolVars[scvIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
                    const Scalar xN2wEquil    = elemVolVars[scvIdx].xEquil(wPhaseIdx, nCompIdx);

                    // actual output into the text file
                    // This could be done much more efficiently
                    // if writing to the text file would be done all at once.
                    dataFile.open(fileName.c_str() , std::ios::app);
                    dataFile << timeStepIndex
                            <<"\t\t\t"
                            << time
                            << "\t\t\t"
                            << globalPosCurrent[0]
                            << "\t\t\t"
                            << globalPosCurrent[1];

                    dataFile <<"\t\t\t"
                                << saturationW
                                <<"\t\t" << Tw
                                <<"\t\t" << Tn
                                <<"\t\t" << Ts
                                <<"\t\t" << xH2On
                                <<"\t\t" << xH2OnEquil
                                <<"\t\t" << xN2w
                                <<"\t\t" << xN2wEquil
                                ;
                    dataFile <<"\n";
                    dataFile.close();
                }
            }
        }
    return ;
    }

    /*!
     * \brief   Check whether the current point is on a line between two points
     */
    const bool isBetween(const GlobalPosition & globalPosCurrent,
                         const GlobalPosition & pointOne,
                         const GlobalPosition & pointTwo) const
    {
        if (    pointOne[0] - globalPosCurrent[0] <= 1e-5
            and pointOne[1] - globalPosCurrent[1] <= 1e-5
            and globalPosCurrent[0] - pointTwo[0] <= 1e-5
            and globalPosCurrent[1] - pointTwo[1] <= 1e-5 )
            return true;

        return false;
    }
};

}
#endif