HeartSolver.hpp 13.7 KB
Newer Older
thomaskummer's avatar
thomaskummer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
//
//  HeartSolver.cpp
//  LifeV
//
//  Created by Thomas Kummer on 03.05.16.
//  Copyright © 2016 Thomas Kummer. All rights reserved.
//

#ifndef _HEARTDATA_H_
#define _HEARTDATA_H_


#include <stdio.h>
thomaskummer's avatar
thomaskummer committed
14
#include <lifev/em/solver/EMSolver.hpp>
thomaskummer's avatar
thomaskummer committed
15
#include <lifev/em/solver/circulation/Circulation.hpp>
thomaskummer's avatar
thomaskummer committed
16
17


thomaskummer's avatar
thomaskummer committed
18
19
20
21
22

namespace LifeV
{

    
thomaskummer's avatar
thomaskummer committed
23
24
25
class HeartData
{
public:
thomaskummer's avatar
thomaskummer committed
26
    
thomaskummer's avatar
thomaskummer committed
27
28
    HeartData() {}
    
thomaskummer's avatar
thomaskummer committed
29
    virtual ~HeartData() {};
thomaskummer's avatar
thomaskummer committed
30
    
thomaskummer's avatar
thomaskummer committed
31
32
33
34
35
36
    void setup(const GetPot& datafile)
    {
        setGetPotFile(datafile);
        readFile();
    }
    
thomaskummer's avatar
thomaskummer committed
37
    const Real& dt_activation () const { return M_dt_activation; }
thomaskummer's avatar
thomaskummer committed
38
    const Real& dt_loadstep () const { return M_dt_loadstep; }
thomaskummer's avatar
thomaskummer committed
39
    const Real& activationLimit_loadstep () const { return M_activationLimit_loadstep; }
thomaskummer's avatar
thomaskummer committed
40
41
    const Real& dt_mechanics () const { return M_dt_mechanics; }
    const Real& dt_save () const { return M_dt_save; }
thomaskummer's avatar
thomaskummer committed
42
    const Real& endtime () const { return M_endtime; }
thomaskummer's avatar
thomaskummer committed
43
44
45
46
    const UInt& mechanicsLoadstepIter () const { return M_mechanicsLoadstepIter; }
    const UInt& mechanicsCouplingIter () const { return M_mechanicsCouplingIter; }
    const UInt& maxiter () const { return M_maxiter; }
    const UInt& preloadSteps () const { return M_preloadSteps; }
thomaskummer's avatar
thomaskummer committed
47
    const bool& safePreload () const { return M_safePreload; }
thomaskummer's avatar
thomaskummer committed
48

thomaskummer's avatar
thomaskummer committed
49
50
51
52
53
54
55
    const UInt& pPerturbationFe () const { return M_pPerturbationFe; }
    const UInt& pPerturbationCirc () const { return M_pPerturbationCirc; }
    const UInt& couplingError () const { return M_couplingError; }
    const UInt& couplingJFeSubIter () const { return M_couplingJFeSubIter; }
    const UInt& couplingJFeSubStart () const { return M_couplingJFeSubStart; }
    const UInt& couplingJFeIter () const { return M_couplingJFeIter; }

thomaskummer's avatar
thomaskummer committed
56
57
    const GetPot& datafile () { return M_datafile; }
    
thomaskummer's avatar
thomaskummer committed
58
protected:
thomaskummer's avatar
thomaskummer committed
59
60
    
    
thomaskummer's avatar
thomaskummer committed
61
    void setGetPotFile(const GetPot& datafile)
thomaskummer's avatar
thomaskummer committed
62
    {
thomaskummer's avatar
thomaskummer committed
63
        M_datafile = datafile;
thomaskummer's avatar
thomaskummer committed
64
65
    }
    
thomaskummer's avatar
thomaskummer committed
66
    void readFile()
thomaskummer's avatar
thomaskummer committed
67
    {
thomaskummer's avatar
thomaskummer committed
68
        M_dt_activation = M_datafile ("activation/time_discretization/timestep", 0.05 );
thomaskummer's avatar
thomaskummer committed
69
70
        M_dt_loadstep =  M_datafile ( "solid/time_discretization/dt_loadstep", 1.0 );
        M_activationLimit_loadstep =  M_datafile ( "solid/time_discretization/activation_limit_loadstep", 0.0 );
thomaskummer's avatar
thomaskummer committed
71
        M_dt_mechanics = M_datafile ("solid/time_discretization/timestep", 1.0 );
thomaskummer's avatar
thomaskummer committed
72
        M_dt_save = M_datafile ( "exporter/save", 10. );
thomaskummer's avatar
thomaskummer committed
73
        M_endtime = M_datafile ("solid/time_discretization/endtime", 100000);
thomaskummer's avatar
thomaskummer committed
74
75
76
        M_mechanicsLoadstepIter = static_cast<UInt>( M_dt_loadstep / M_dt_activation );
        M_mechanicsCouplingIter = static_cast<UInt>( M_dt_mechanics / M_dt_activation );
        M_maxiter = static_cast<UInt>( M_endtime / M_dt_activation ) ;
thomaskummer's avatar
thomaskummer committed
77
        M_preloadSteps = M_datafile ( "solid/boundary_conditions/numPreloadSteps", 0);
thomaskummer's avatar
thomaskummer committed
78
        M_safePreload = M_datafile ( "exporter/savePreload", false );
thomaskummer's avatar
thomaskummer committed
79

thomaskummer's avatar
thomaskummer committed
80
81
82
83
84
85
        M_pPerturbationFe = M_datafile ( "solid/coupling/pPerturbationFe", 1e-2 );
        M_pPerturbationCirc = M_datafile ( "solid/coupling/pPerturbationCirc", 1e-3 );
        M_couplingError = M_datafile ( "solid/coupling/couplingError", 1e-6 );
        M_couplingJFeSubIter = M_datafile ( "solid/coupling/couplingJFeSubIter", 1 );
        M_couplingJFeSubStart = M_datafile ( "solid/coupling/couplingJFeSubStart", 1 );
        M_couplingJFeIter = M_datafile ( "solid/coupling/couplingJFeIter", 1 );
thomaskummer's avatar
thomaskummer committed
86
        
thomaskummer's avatar
thomaskummer committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
    //        const Real dpMax = dataFile ( "solid/coupling/dpMax", 0.1 );
    //
    //        std::vector<std::vector<std::string> > bcNames { { "lv" , "p" } , { "rv" , "p" } };
    //        std::vector<double> bcValues { p ( "lv" ) , p ( "rv") };
    //        std::vector<double> bcValuesPre ( bcValues );
    //
    //        VectorSmall<4> ABdplv, ABdprv, ABcoef;
    //        ABcoef (0) = 55/24; ABcoef (1) = -59/24; ABcoef (2) = 37/24; ABcoef (3) = -3/8;
    //
    //        VectorSmall<2> VCirc, VCircNew, VCircPert, VFe, VFeNew, VFePert, R, dp;
    //        MatrixSmall<2,2> JFe, JCirc, JR;
    //        
    //        UInt iter (0);
    //        Real t (0);
    
thomaskummer's avatar
thomaskummer committed
102
103
104
        
    }
    
thomaskummer's avatar
thomaskummer committed
105
106
107
108
109
110
111
112
113
    Real M_dt_activation;
    Real M_dt_loadstep;
    Real M_activationLimit_loadstep;
    Real M_dt_mechanics;
    Real M_dt_save;
    Real M_endtime;
    UInt M_mechanicsLoadstepIter;
    UInt M_mechanicsCouplingIter;
    UInt M_maxiter;
thomaskummer's avatar
thomaskummer committed
114
    UInt M_preloadSteps;
thomaskummer's avatar
thomaskummer committed
115
116
    bool M_safePreload;

thomaskummer's avatar
thomaskummer committed
117
118
119
120
121
122
    Real M_pPerturbationFe;
    Real M_pPerturbationCirc;
    Real M_couplingError;
    UInt M_couplingJFeSubIter;
    UInt M_couplingJFeSubStart;
    UInt M_couplingJFeIter;
thomaskummer's avatar
thomaskummer committed
123
    
thomaskummer's avatar
thomaskummer committed
124
    GetPot M_datafile;
thomaskummer's avatar
thomaskummer committed
125
    
thomaskummer's avatar
thomaskummer committed
126
};
thomaskummer's avatar
thomaskummer committed
127
    
thomaskummer's avatar
thomaskummer committed
128
    
thomaskummer's avatar
thomaskummer committed
129
130
131
132
template <class EmSolver>
class HeartSolver {
   
public:
thomaskummer's avatar
thomaskummer committed
133
    
thomaskummer's avatar
thomaskummer committed
134
135
    HeartSolver(EmSolver& emSolver,  Circulation& circulationSolver) :
        M_emSolver          (emSolver),
thomaskummer's avatar
thomaskummer committed
136
137
        M_circulationSolver (circulationSolver),
        M_heartData         (HeartData())
thomaskummer's avatar
thomaskummer committed
138
    {}
thomaskummer's avatar
thomaskummer committed
139
    
thomaskummer's avatar
thomaskummer committed
140
    virtual ~HeartSolver() {}
thomaskummer's avatar
thomaskummer committed
141
142
143
144
145
146
147
148
149
150

    EmSolver& emSolver()
    {
        return M_emSolver;
    }
    
    Circulation& circulation()
    {
        return M_heartData;
    }
thomaskummer's avatar
thomaskummer committed
151
    
thomaskummer's avatar
thomaskummer committed
152
    const HeartData& data() const
thomaskummer's avatar
thomaskummer committed
153
154
155
156
    {
        return M_heartData;
    }
    
thomaskummer's avatar
thomaskummer committed
157
    void setup(const GetPot& datafile)
thomaskummer's avatar
thomaskummer committed
158
    {
thomaskummer's avatar
thomaskummer committed
159
        M_heartData.setup(datafile);
thomaskummer's avatar
thomaskummer committed
160
    }
thomaskummer's avatar
thomaskummer committed
161
    
thomaskummer's avatar
thomaskummer committed
162
    template <class lambda>
thomaskummer's avatar
thomaskummer committed
163
    void preload(const lambda& modifyFeBC, const std::vector<Real>& bcValues)
thomaskummer's avatar
thomaskummer committed
164
    {
thomaskummer's avatar
thomaskummer committed
165
        M_emSolver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(0.0);
thomaskummer's avatar
thomaskummer committed
166
167
168
169
170
171
172
173
174
175
        
        auto preloadPressure = [] (std::vector<double> p, const int& step, const int& steps)
        {
            for (auto& i : p) {i *= double(step) / double(steps);}
            return p;
        };
        
        LifeChrono chronoSave;
        chronoSave.start();
        
thomaskummer's avatar
thomaskummer committed
176
        M_emSolver.saveSolution (-1.0);
thomaskummer's avatar
thomaskummer committed
177
        
thomaskummer's avatar
thomaskummer committed
178
        if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
179
180
181
182
183
184
185
186
187
188
189
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nData stored in " << chronoSave.diff() << " s";
            std::cout << "\n*****************************************************************\n";
        }
        
        LifeChrono chronoPreload;
        chronoPreload.start();
        
        for (int i (1); i <= data().preloadSteps(); i++)
        {
thomaskummer's avatar
thomaskummer committed
190
            if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
191
192
193
194
195
196
197
            {
                std::cout << "\n*****************************************************************";
                std::cout << "\nPreload step: " << i << " / " << data().preloadSteps();
                std::cout << "\n*****************************************************************\n";
            }
            
            // Update pressure b.c.
thomaskummer's avatar
thomaskummer committed
198
            modifyFeBC(preloadPressure(bcValues, i, data().preloadSteps() ));
thomaskummer's avatar
thomaskummer committed
199
200
            
            // Solve mechanics
thomaskummer's avatar
thomaskummer committed
201
202
            M_emSolver.bcInterfacePtr() -> updatePhysicalSolverVariables();
            M_emSolver.solveMechanics();
thomaskummer's avatar
thomaskummer committed
203
204
            
            // Safe preload steps
thomaskummer's avatar
thomaskummer committed
205
            if ( data().safePreload() ) M_emSolver.saveSolution (i-1);
thomaskummer's avatar
thomaskummer committed
206
207
        }
        
thomaskummer's avatar
thomaskummer committed
208
        if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
209
210
211
212
213
214
215
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nPreload done in: " << chronoPreload.diff();
            std::cout << "\n*****************************************************************\n";
        }

    }
thomaskummer's avatar
thomaskummer committed
216
217
    
    
thomaskummer's avatar
thomaskummer committed
218
    void restart(std::string& restartInput, const GetPot& command_line)
thomaskummer's avatar
thomaskummer committed
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
    {
        const std::string restartDir = command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
        
        Real dtExport = 10.;
        
        // Get most recent restart index
        if ( restartInput == "." )
        {
            restartInput = pipeToString( ("tail -n 1 " + restartDir + "solution.dat | awk -F '[. ]' '{print $1 \".\" $2}' | awk '{printf \"%05g\", int($1*1000/" + std::to_string(dtExport) + ") + 1}'").c_str() );
        }
        
        // Set time variable
        const unsigned int restartInputStr = std::stoi(restartInput);
        const unsigned int nIter = (restartInputStr - 1) * dtExport / data().dt_mechanics();
        t = nIter * data().dt_mechanics();
        
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        
        // Set time exporter time index
        M_emSolver.setTimeIndex(restartInputStr + 1);
        //solver.importHdf5();
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        // Load restart solutions from output files
        std::string polynomialDegree = data().datafile() ( "solid/space_discretization/order", "P1");
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        ElectrophysiologyUtility::importVectorField ( M_emSolver.structuralOperatorPtr() -> displacementPtr(), "MechanicalSolution" , "displacement", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        for ( unsigned int i = 0; i < M_emSolver.electroSolverPtr()->globalSolution().size() ; ++i )
        {
            ElectrophysiologyUtility::importScalarField (M_emSolver.electroSolverPtr()->globalSolution().at(i), "ElectroSolution" , ("Variable" + std::to_string(i)), M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        }
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        ElectrophysiologyUtility::importScalarField (M_emSolver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        //ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "ActivationTimeSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
        circulationSolver.restartFromFile ( restartDir + "solution.dat" , nIter );
    }
    
    
thomaskummer's avatar
thomaskummer committed
276
    
thomaskummer's avatar
thomaskummer committed
277
protected:
thomaskummer's avatar
thomaskummer committed
278
    
thomaskummer's avatar
thomaskummer committed
279
    
thomaskummer's avatar
thomaskummer committed
280
281
    EmSolver& M_emSolver;
    Circulation& M_circulationSolver;
thomaskummer's avatar
thomaskummer committed
282
    
thomaskummer's avatar
thomaskummer committed
283
    HeartData M_heartData;
thomaskummer's avatar
thomaskummer committed
284
    
thomaskummer's avatar
thomaskummer committed
285
    
thomaskummer's avatar
thomaskummer committed
286
287
288
289
    VectorSmall<2> M_pressure;
    VectorSmall<2> M_volume;

    
thomaskummer's avatar
thomaskummer committed
290

thomaskummer's avatar
thomaskummer committed
291
    
thomaskummer's avatar
thomaskummer committed
292
293
294
295
296
297
298
299
300
301
302
    std::string pipeToString ( const char* command )
    {
        FILE* file = popen( command, "r" ) ;
        std::ostringstream stm ;
        char line[6] ;
        fgets( line, 6, file );
        stm << line;
        pclose(file) ;
        return stm.str() ;
    };
    
thomaskummer's avatar
thomaskummer committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
    Real patchForce (const Real& t, const Real& Tmax, const Real& tmax, const Real& tduration) const
    {
        bool time ( fmod(t-tmax+0.5*tduration, 700.) < tduration && fmod(t-tmax+0.5*tduration, 700.) > 0);
        Real force = std::pow( std::sin(fmod(t-tmax+0.5*tduration, 700.)*3.14159265359/tduration) , 2 ) * Tmax;
        return ( time ? force : 0 );
    }
    
    Real patchFunction (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/) const
    {
        Real disp = std::pow( std::sin(fmod(t, 700.) * 3.14159265359/300) , 2 )*15;
        return disp;
    }
    
    Real Iapp (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/) const
    {
        bool coords ( Y < -7. );
        //bool coords ( Y > 4. ); //( Y > 1.5 && Y < 3 );
        bool time ( fmod(t, 700.) < 4 && fmod(t, 700.) > 2);
        return ( coords && time ? 30 : 0 );
    }
    
    Real potentialMultiplyerFcn (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/) const
    {
        bool time ( fmod(t, 700.) < 4 && fmod(t, 700.) > 2);
        return 1.4 * time; // ( Y < 2.5 && Y > 0.5 ? 1.0 : 0.0 );
    }
    
thomaskummer's avatar
thomaskummer committed
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    Real
    externalPower ( const VectorEpetra& dispCurrent,
                    const VectorEpetra& dispPrevious,
                    const boost::shared_ptr<ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > > dispETFESpace,
                    Real pressure,
                    Real dt,
                    const unsigned int bdFlag) const
    {
        VectorEpetra traction ( dispCurrent.map() );
        VectorEpetra velocity ( (dispCurrent - dispPrevious) / dt );
        
        MatrixSmall<3,3> Id;
        Id(0,0) = 1.; Id(0,1) = 0., Id(0,2) = 0.;
        Id(1,0) = 0.; Id(1,1) = 1., Id(1,2) = 0.;
        Id(2,0) = 0.; Id(2,1) = 0., Id(2,2) = 1.;
        
        {
            using namespace ExpressionAssembly;
            
            auto I = value(Id);
            auto Grad_u = grad( dispETFESpace, dispCurrent, 0);
            auto F =  Grad_u + I;
            auto FmT = minusT(F);
            auto J = det(F);
            auto p = value(pressure);
            
            QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
            
            integrate ( boundary ( dispETFESpace->mesh(), bdFlag),
                       myBDQR,
                       dispETFESpace,
                       p * J * dot( FmT * Nface,  phi_i)
                       //p * J * dot( FmT * Nface,  phi_i)
                       //value(-1.0) * J * dot (vE1, FmT * Nface) * phi_i) >> intergral
                       ) >> traction;
            
            traction.globalAssemble();
        }
        
        return traction.dot(velocity);
    }
    
    
thomaskummer's avatar
thomaskummer committed
373
374
375
376
};

}

thomaskummer's avatar
thomaskummer committed
377
#endif