main.cpp 29.4 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "inmost.h"
using namespace INMOST;

#ifndef M_PI
#define M_PI 3.141592653589
#endif

#if defined(USE_MPI)
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
#else
#define BARRIER
#endif

//shortcuts
Kirill Terekhov's avatar
Kirill Terekhov committed
20
typedef Storage::bulk bulk;
Kirill Terekhov's avatar
Kirill Terekhov committed
21
22
23
24
25
26
27
28
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
typedef Storage::real_array real_array;
typedef Storage::var_array var_array;



Kirill Terekhov's avatar
Kirill Terekhov committed
29
//#define OPTIMIZATION
Kirill Terekhov's avatar
Kirill Terekhov committed
30
31
32

int main(int argc,char ** argv)
{
Kirill Terekhov's avatar
Kirill Terekhov committed
33
    Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
34
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
35
    Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
36
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
37
    if( argc > 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
38
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
39
        /*
Kirill Terekhov's avatar
Kirill Terekhov committed
40
41
42
43
44
        std::cout << "Row: " << sizeof(Sparse::Row) << std::endl;
        std::cout << "HessianRow: " << sizeof(Sparse::HessianRow) << std::endl;
        std::cout << "Matrix: " << sizeof(Sparse::Matrix) << std::endl;
        std::cout << "HessianMatrix: " << sizeof(Sparse::HessianMatrix) << std::endl;
        std::cout << "variable: " << sizeof(multivar_expression) << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
45
46
47
48
49
50
51
52
53
54
55
56
        Sparse::Row J;
        Sparse::HessianRow H;
        unknown a = unknown(0.5,1);
        unknown b = unknown(0.1,2);
        (a*a*b*(a+b)).GetHessian(1.0,J,1.0,H);
        std::cout << "J: "; J.Print();
        std::cout << "H: "; H.Print();
        J.Clear();
        (a*a*b*(a+b)).GetJacobian(1.0,J);
        std::cout << "J: "; J.Print();
        return 0;
         */
Kirill Terekhov's avatar
Kirill Terekhov committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
        double ttt; // Variable used to measure timing
        bool repartition = false; // Is it required to redistribute the mesh?
        Mesh * m = new Mesh(); // Create an empty mesh
        { // Load the mesh
            ttt = Timer();
            m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
            if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
            if( m->isParallelFileFormat(argv[1]) ) //The format is
            {
                m->Load(argv[1]); // Load mesh from the parallel file format
                repartition = true; // Ask to repartition the mesh
            }
            else if( m->GetProcessorRank() == 0 ) m->Load(argv[1]); // Load mesh from the serial file format
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Load the mesh: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
72
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
73
74
75
76
        
        
#if defined(USE_PARTITIONER)
        if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
Kirill Terekhov's avatar
Kirill Terekhov committed
77
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
            { // Compute mesh partitioning
                ttt = Timer();
                Partitioner p(m); //Create Partitioning object
                p.SetMethod(Partitioner::Inner_RCM,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
                p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;
            }
            
            { //Distribute the mesh
                ttt = Timer();
                m->Redistribute(); // Redistribute the mesh data
                m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
                BARRIER
                if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
            }
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
95
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
96
97
98
99
100
101
102
103
104
105
106
107
        
        { // prepare geometrical data on the mesh
            ttt = Timer();
            Mesh::GeomParam table;
            table[CENTROID]    = CELL | FACE; //Compute averaged center of mass
            table[NORMAL]      = FACE;        //Compute normals
            table[ORIENTATION] = FACE;        //Check and fix normal orientation
            table[MEASURE]     = CELL | FACE; //Compute volumes and areas
            //table[BARYCENTER]  = CELL | FACE; //Compute volumetric center of mass
            m->PrepareGeometricData(table); //Ask to precompute the data
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
108
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
109
110
111
112
113
114
115
        
        // data tags for
        Tag tag_P;  // Pressure
        Tag tag_K;  // Diffusion tensor
        Tag tag_F;  // Forcing term
        Tag tag_BC; // Boundary conditions
        Tag tag_W;  // Gradient matrix acting on harmonic points on faces and returning gradient on faces
Kirill Terekhov's avatar
Kirill Terekhov committed
116
        Tag tag_DMP; // Indicates weather local W matrix satisfy DMP condition
Kirill Terekhov's avatar
Kirill Terekhov committed
117
118
119
120
121
122
123
        
        if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
        { // Exchange ghost cells
            ttt = Timer();
            m->ExchangeGhost(1,FACE); // Produce layer of ghost cells
            BARRIER
            if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
124
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
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
        
        { //initialize data
            if( m->HaveTag("PERM") ) // is diffusion tensor already defined on the mesh? (PERM from permeability)
                tag_K = m->GetTag("PERM"); // get the diffusion tensor
            
            if( !tag_K.isValid() || !tag_K.isDefined(CELL) ) // diffusion tensor was not initialized or was not defined on cells.
            {
                tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,6); // create a new tag for symmetric diffusion tensor K
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) // loop over mesh cells
                {
                    Cell cell = m->CellByLocalID(q);
                    real_array K = cell->RealArray(tag_K);
                    // assign a symmetric positive definite tensor K
                    K[0] = 1.0; //XX
                    K[1] = 0.0; //XY
                    K[2] = 0.0; //XZ
                    K[3] = 1.0; //YY
                    K[4] = 0.0; //YZ
                    K[5] = 1.0; //ZZ
                }
                
                m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
            }
            
            if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
                tag_P = m->GetTag("PRESSURE"); //Get the pressure
            
            if( !tag_P.isValid() || !tag_P.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes
            {
                srand(1); // Randomization
                tag_P = m->CreateTag("PRESSURE",DATA_REAL,CELL|FACE,NONE,1); // Create a new tag for the pressure
                for(Mesh::iteratorElement e = m->BeginElement(CELL|FACE); e != m->EndElement(); ++e) //Loop over mesh cells
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
            
            if( !tag_P.isDefined(FACE) )
            {
                tag_P = m->CreateTag("PRESSURE",DATA_REAL,FACE,NONE,1);
                for(Mesh::iteratorElement e = m->BeginElement(FACE); e != m->EndElement(); ++e) //Loop over mesh cells
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
            
            
            
            if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh?
            {
                tag_BC = m->GetTag("BOUNDARY_CONDITION");
                //initialize unknowns at boundary
            }
            m->ExchangeData(tag_P,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
            tag_W = m->CreateTag("nKGRAD",DATA_REAL,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
176
            tag_DMP = m->CreateTag("isDMP",DATA_BULK,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
177
178
            ttt = Timer();
            //Assemble gradient matrix W on cells
Kirill Terekhov's avatar
Kirill Terekhov committed
179
            int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
180
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
181
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
182
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
183
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
Kirill Terekhov's avatar
Kirill Terekhov committed
184
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
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
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
290
291
292
293
294
295
296
297
298
299
300
301
302
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
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
                Cell cell = m->CellByLocalID(q);
                real xP[3]; //center of the cell
                real yF[3]; //center of the face
                real nF[3]; //normal to the face
                real aF; //area of the face
                real vP = cell->Volume(); //volume of the cell
                cell->Centroid(xP);
                ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
                int NF = (int)faces.size(); //number of faces;
                rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
                //rMatrix U,S,V;
                //K0.SVD(U,S,V);
                //for(int k = 0; k < 3; ++k) S(k,k) = sqrt(S(k,k));
                //rMatrix K = U*S*V;
                rMatrix nKGRAD(NF,NF), NK(NF,3), R(NF,3), D(NF,NF), U(NF,NF), Areas(NF,1); //big gradient matrix, co-normals, directions
                for(int k = 0; k < NF; ++k) //loop over faces
                {
                    aF = faces[k].Area();
                    faces[k].Centroid(yF);
                    faces[k].OrientedUnitNormal(cell->self(),nF);
                    // assemble matrix of directions
                    R(k,0) = (yF[0]-xP[0])*aF;
                    R(k,1) = (yF[1]-xP[1])*aF;
                    R(k,2) = (yF[2]-xP[2])*aF;
                    // assemble matrix of co-normals
                    rMatrix nK = rMatrix::FromVector(nF,3).Transpose()*K;
                    NK(k,0) = nK(0,0);
                    NK(k,1) = nK(0,1);
                    NK(k,2) = nK(0,2);
                } //end of loop over faces
                rMatrix SU,SS,SV;
                nKGRAD = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part
                /*
                 std::cout << "W" << std::endl;
                 nKGRAD.Print();
                 nKGRAD.SVD(SU,SS,SV);
                 std::cout << "U" << std::endl;
                 SU.Print();
                 std::cout << "S" << std::endl;
                 SS.Print();
                 std::cout << "V" << std::endl;
                 SV.Print();
                 std::cout << "Check " << (nKGRAD - SU*SS*SV.Transpose()).FrobeniusNorm() << std::endl;
                 */
                
                int rank = 0; //size of matrix U
                
                { //Retrive orthogonal to R matrix D
                    //Symmetric orthogonal matrix
                    rMatrix DUD = (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose());
                    //perfrom singular value decomposition
                    //S should be unity matrix with rank NF-3
                    rMatrix DUD_U,DUD_S,DUD_V;
                    DUD.SVD(DUD_U,DUD_S,DUD_V);
                    //compute the rank
                    for(int q = 0; q < NF; ++q)
                        if( DUD_S(q,q) > 1.0e-2 )
                            rank++;
                    rank = NF-3;
                    if( rank != NF-3)
                    {
                        std::cout << "rank: " << rank << " expected " << NF-3 << std::endl;
                        DUD_S.Print();
                    }
                    //chop matrix to the full rank
                    DUD_S.RemoveSubset(rank,NF,rank,NF);
                    DUD_V.RemoveColumns(rank,NF);
                    //assign the matrix
                    D = DUD_V;
                    U = DUD_S;
                }
                //std::cout << "D" << std::endl;
                //D.Print();
                //std::cout << "U" << std::endl;
                //U.Print();
                //std::cout << "DtR" << std::endl;
                //(D.Transpose()*R).Print();
                
                U *=(2.0/(static_cast<real>(NF)*vP)*(NK*K.Invert(true).first*NK.Transpose()).Trace());
#if defined(OPTIMIZATION)
                { //Make W a Z-matrix
                    vMatrix vL(rank,rank), vD(rank,rank), vW(NF,NF);
                    int unk = 0;
                    // U = L*D*L^T
                    vD.Zero();
                    //diagonal D
                    for(int i = 0; i < rank; ++i)
                    {
                        vD(i,i) = unknown(U(i,i),unk);
                        unk++;
                    }
                    //off-diagonal
                    vL.Zero();
                    for(int i = 0; i < rank; ++i)
                        vL(i,i) = 1.0;
                    for(int i = 1; i < rank; ++i)
                        for(int j = 0; j < i; ++j)
                        {
                            vL(i,j) = unknown(0.0,unk);
                            unk++;
                        }
                    //std::cout << "unknowns: " << unk << std::endl;
                    variable phi,s;
                    //std::cout << "vD" << std::endl;
                    //vD.Print();
                    //std::cout << "vL" << std::endl;
                    //vL.Print();
                    
                    int iter = 0;
                    do
                    { //Optimize U matrix
                        vW = nKGRAD + D*vL*vD*vL.Transpose()*D.Transpose();
                        //construct minimization functional phi(W)
                        phi = 0.0;
                        for(int i = 0; i < NF; ++i)
                        {
                            //phi += 1.0 / (vW(i,i)*vW(i,i));
                            s = vW(i,i)*faces[i].Area();
                            for(int j = 0; j < NF; ++j) if( i != j )
                            {
                                s += vW(i,j)*faces[j].Area();
                                phi += (vW(i,j)+fabs(vW(i,j)))*(vW(i,j)+fabs(vW(i,j)));
                            }
                            phi += (s - fabs(s))*(s - fabs(s));
                        }
                        Sparse::Row & der = phi.GetRow(); //row of derivatives
                        //std::sort(der.Begin(),der.End());
                        //for(int i = 0; i < der.Size(); ++i)
                        //    std::cout << "(" << der.GetIndex(i) << "," << der.GetValue(i) << ") ";
                        //std::cout<<std::endl;
                        int q = 0;
                        real a = 0.00005;
                        real minvD = 1.0e20;
                        //diagonal
                        for(int i = 0; i < rank; ++i)
                        {
                            real d = a*der[q++];
                            //if( vD(i,i)-d > 0.0 )
                            vD(i,i) -= d;
                            if( vD(i,i) < minvD ) minvD = get_value(vD(i,i));
                        }
                        std::cout << "[" << iter << "] phi: " << get_value(phi) << " minD " << minvD << std::endl;
                        //off-diagonal
                        for(int i = 1; i < rank; ++i)
                            for(int j = 0; j < i; ++j)
                            {
                                real d = a*der[q++];
                                vL(i,j) -= d;
                            }
                        iter++;
                        //std::cout << "vD" << std::endl;
                        //vD.Print();
                        //std::cout << "vL" << std::endl;
                        //vL.Print();
                    } while(iter < 100 && phi > 1.0e-3);
                    {
                        vMatrix vU = vL*vD*vL.Transpose();
                        for(int i = 0; i < rank; ++i)
                            for(int j = 0; j < rank; ++j)
                                U(i,j) = get_value(vU(i,j));
                    }
                    //std::cout << "U: " << std::endl;
                    //U.Print();
                    
                    
                }
#endif
                //std::cout << "UDtR" << std::endl;
                //(U*D.Transpose()*R).Print();
                
                nKGRAD += D*U*D.Transpose();
                
                //std::cout << "W: " << std::endl;
                //nKGRAD.Print();
Kirill Terekhov's avatar
Kirill Terekhov committed
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
                bulk & isDMP = cell->Bulk(tag_DMP);
                isDMP = 1;
                for(int k = 0; k < NF; ++k)
                {
                    real row_sum = 0;
                    if( nKGRAD(k,k) < 0.0 ) isDMP = 0;
                    for(int j = 0; j < NF; ++j)
                        row_sum += nKGRAD(k,j);
                    if( row_sum < 0.0 ) isDMP = 0;
                    for(int j = k+1; j < NF; ++j)
                        if( nKGRAD(k,j) > 0.0 )
                            isDMP = 0;
                }
                ++total;
                if( isDMP ) ++dmp;
Kirill Terekhov's avatar
Kirill Terekhov committed
374
375
376
377
378
                real_array W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
                W.resize(NF*NF); //resize the structure
                std::copy(nKGRAD.data(),nKGRAD.data()+NF*NF,W.data()); //write down the gradient matrix
            } //end of loop over cells
            std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
379
            std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
            
            if( m->HaveTag("FORCE") ) //Is there force on the mesh?
            {
                tag_F = m->GetTag("FORCE"); //initial force
                assert(tag_F.isDefined(CELL)); //assuming it was defined on cells
            } // end of force
        } //end of initialize data
        
        std::cout << "Initialization done" << std::endl;
        
        
        integer nit = 0;
        ttt = Timer();
        
        { //Main loop for problem solution
            Automatizator aut(m); // declare class to help manage unknowns
            Automatizator::MakeCurrent(&aut);
            dynamic_variable P(aut,aut.RegisterDynamicTag(tag_P,CELL|FACE)); //register pressure as primary unknown
            aut.EnumerateDynamicTags(); //enumerate all primary variables
            
            std::cout << "Enumeration done" << std::endl;
            
            
            Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex());
            Sparse::Vector Update  ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
            {//Annotate matrix
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                {
                    Cell cell = m->CellByLocalID(q);
                    if( cell.GetStatus() != Element::Ghost )
                        Text.SetAnnotation(P.Index(cell),"Cell-centered pressure value");
                }
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    if( face.GetStatus() != Element::Ghost )
                    {
                        if( tag_BC.isValid() && face.HaveData(tag_BC) )
                            Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
                        else
                            Text.SetAnnotation(P.Index(face),"Interface pressure");
                    }
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
425
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
426
427
428
429
430
431
432
            
            std::cout << "Matrix was annotated" << std::endl;
            
            do
            {
                R.Clear(); //clean up the residual
                double tttt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
433
                int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
434
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
435
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
436
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
437
438
439
440
441
442
                for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
                {
                    Cell cell = m->CellByLocalID(q);
                    ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
                    int NF = (int)faces.size();
                    rMatrix nKGRAD(cell->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
Kirill Terekhov's avatar
Kirill Terekhov committed
443
                    bulk & isDMP = cell->Bulk(tag_DMP);
Kirill Terekhov's avatar
Kirill Terekhov committed
444
445
446
447
448
                    vMatrix pF(NF,1); //vector of pressure differences on faces
                    vMatrix FLUX(NF,1); //computed flux on faces
                    for(int k = 0; k < NF; ++k)
                        pF(k,0) = (P(faces[k]) - P(cell))*faces[k].Area();
                    FLUX = nKGRAD*pF; //fluxes on faces
Kirill Terekhov's avatar
Kirill Terekhov committed
449
450
451
452
453
454
455
456
457
458
                    //Change matrix by nonlinear correction for DMP
                    ++total;
                    
                    if( !isDMP )
                    //if( false )
                    {
                        const real var = 0;//0.25;
                        vMatrix vnKGRAD = nKGRAD;
                        vMatrix card(NF,1);
                        vMatrix dpF(NF,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
459
                        variable div = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
460
461
462
463
                        //denominators of corrections
                        //this is sum of absolute values of differences
                        for(int k = 0; k < NF; ++k)
                        {
Kirill Terekhov's avatar
Kirill Terekhov committed
464
                            div += FLUX(k,0)*faces[k].Area();;
Kirill Terekhov's avatar
Kirill Terekhov committed
465
466
467
468
469
470
471
472
473
474
475
476
477
                            for(int j = 0; j < NF; ++j)
                            {
                                dpF(k,0) += soft_fabs(pF(j,0) - pF(k,0),1.0e-9);
                                card(k,0) += (soft_fabs(pF(j,0) - pF(k,0),1.0e-9)) / (soft_fabs(pF(j,0) - pF(k,0),1.0e-9)+1.0e-9);
                            }
                        }
                        //actual corrections
                        variable beta;
                        for(int k = 0; k < NF; ++k)
                        {
                            for(int j = k+1; j < NF; ++j) //if( nKGRAD(k,j) > 0.0 )
                            {
                                //beta = soft_fabs(FLUX(k,0),1.0e-9)/(dpF(k,0)) + soft_fabs(FLUX(j,0),1.0e-9)/(dpF(j,0));
Kirill Terekhov's avatar
Kirill Terekhov committed
478
479
480
                                beta = soft_max(soft_fabs(FLUX(k,0)-FLUX(j,0),1.0e-9)/card(k,0),soft_fabs(FLUX(j,0)-FLUX(k,0),1.0e-9)/card(j,0),1.0e-9)/soft_fabs(pF(k,0)-pF(j,0),1.0e-9);
                                //beta = soft_max(soft_fabs(FLUX(k,0),1.0e-9)/card(k,0),soft_fabs(FLUX(j,0),1.0e-9)/card(j,0),1.0e-9)/soft_fabs(pF(k,0)-pF(j,0),1.0e-9);
                                //beta = soft_max(soft_fabs(div,1.0e-9)/card(k,0),soft_fabs(div,1.0e-9)/card(j,0),1.0e-9)/soft_fabs(pF(k,0)-pF(j,0),1.0e-9);
Kirill Terekhov's avatar
Kirill Terekhov committed
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
                                beta = variation(beta,var);
                                vnKGRAD(k,j) -= (beta);
                                vnKGRAD(j,k) -= (beta);
                                vnKGRAD(k,k) += (beta);
                                vnKGRAD(j,j) += (beta);
                            }
                        }
                        FLUX = vnKGRAD*pF;
                        bool haveDMP = true;
                        for(int k = 0; k < NF; ++k)
                        {
                            real row_sum = 0;
                            if( vnKGRAD(k,k) < 0.0 ) haveDMP = false;
                            for(int j = 0; j < NF; ++j)
                                row_sum += get_value(vnKGRAD(k,j));
                            if( row_sum < 0.0 ) haveDMP = false;
                            for(int j = k+1; j < NF; ++j)
                                if( vnKGRAD(k,j) > 0.0 )
                                    haveDMP = false;
                        }
                        if( !haveDMP )
                        {
                            //std::cout << "Failed to correct, matrix:" << std::endl;
                            //nKGRAD.Print();
                        }
                        else ++dmp;
                    }
                    else ++dmp;
                    
Kirill Terekhov's avatar
Kirill Terekhov committed
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
                    if( cell.GetStatus() != Element::Ghost )
                    {
                        for(int k = 0; k < NF; ++k) //loop over faces of current cell
                            R[P.Index(cell)] += FLUX(k,0)*faces[k].Area();
                    }
                    for(int k = 0; k < NF; ++k) //loop over faces of current cell
                    {
                        if( faces[k].GetStatus() == Element::Ghost ) continue;
                        int index = P.Index(faces[k]);
                        Locks.Lock(index);
                        if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
                        {
                            real_array BC = faces[k].RealArray(tag_BC);
                            R[index] -= BC[0]*P(faces[k]) + BC[1]*FLUX(k,0) - BC[2];
                        }
                        else
                            R[index] -= FLUX(k,0);
                        Locks.UnLock(index);
                    }
                } //end of loop over cells
                
Kirill Terekhov's avatar
Kirill Terekhov committed
531
532
                std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
                
Kirill Terekhov's avatar
Kirill Terekhov committed
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
                
                if( tag_F.isValid() )
                {
#if defined(USE_OMP)
#pragma omp parallel for
#endif
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
                        if( cell.GetStatus() == Element::Ghost ) continue;
                        if( cell->HaveData(tag_F) ) R[P.Index(cell)] += cell->Real(tag_F)*cell->Volume();
                    }
                }
                
                std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
                
                
                R.Rescale();
Kirill Terekhov's avatar
Kirill Terekhov committed
551
552
                //R.GetJacobian().Save("jacobian.mtx",&Text);
                //R.GetResidual().Save("residual.mtx");
Kirill Terekhov's avatar
Kirill Terekhov committed
553
554
555
556
557
558
559
560
561
562
                
                
                std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;
                
                if( R.Norm() < 1.0e-4 ) break;
                
                Solver S(Solver::INNER_MPTILUC);
                S.SetMatrix(R.GetJacobian());
                S.SetParameterReal("relative_tolerance", 1.0e-14);
                S.SetParameterReal("absolute_tolerance", 1.0e-12);
Kirill Terekhov's avatar
Kirill Terekhov committed
563
                S.SetParameterReal("drop_tolerance", 1.0e-2);
Kirill Terekhov's avatar
Kirill Terekhov committed
564
565
566
567
                S.SetParameterReal("reuse_tolerance", 1.0e-4);
                //std::fill(Update.Begin(),Update.End(),0.0);
                if( S.Solve(R.GetResidual(),Update) )
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
568
569
570
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
571
572
573
574
575
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
                        cell->Real(tag_P) -= Update[P.Index(cell)];
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
576
577
578
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
                    for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                    {
                        Face face = m->FaceByLocalID(q);
                        face->Real(tag_P) -= Update[P.Index(face)];
                    }
                    m->ExchangeData(tag_P, CELL|FACE, 0);
                    {
                        std::stringstream str;
                        str << "iter" << nit;
                        if( m->GetProcessorsNumber() == 1 )
                            str << ".vtk";
                        else
                            str << ".pvtk";
                        m->Save(str.str());
                    }
                }
                else
                {
                    std::cout << "Unable to solve: " << S.GetReason() << std::endl;
                    break;
                }
                ++nit;
            } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm
Kirill Terekhov's avatar
Kirill Terekhov committed
602
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
603
604
605
        std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;
        
        if( m->HaveTag("REFERENCE_SOLUTION") )
Kirill Terekhov's avatar
Kirill Terekhov committed
606
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
            Tag tag_E = m->CreateTag("ERRROR",DATA_REAL,CELL,NONE,1);
            Tag tag_R = m->GetTag("REFERENCE_SOLUTION");
            real C, L2, volume;
            C = L2 = volume = 0.0;
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
            {
                Cell cell = m->CellByLocalID(q);
                real err = cell->Real(tag_P) - cell->Real(tag_R);
                real vol = cell->Volume();
                if( C < fabs(err) ) C = fabs(err);
                L2 += err*err*vol;
                volume += vol;
                cell->Real(tag_E) = err;
            }
            L2 = sqrt(L2/volume);
            std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
            C = L2 = volume = 0.0;
            if( tag_R.isDefined(FACE) )
            {
                tag_E = m->CreateTag("ERRROR",DATA_REAL,FACE,NONE,1);
                for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                {
                    Face face = m->FaceByLocalID(q);
                    real err = face->Real(tag_P) - face->Real(tag_R);
                    real vol = (face->BackCell()->Volume() + (face->FrontCell().isValid() ? face->FrontCell()->Volume() : 0))*0.5;
                    if( C < fabs(err) ) C = fabs(err);
                    L2 += err*err*vol;
                    volume += vol;
                    face->Real(tag_E) = err;
                }
                L2 = sqrt(L2/volume);
                std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
            }
            else std::cout << "Reference solution was not defined on faces" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
641
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
642
643
644
645
646
647
648
        
        if( m->GetProcessorsNumber() == 1 )
            m->Save("out.vtk");
        else
            m->Save("out.pvtk");
        
        delete m; //clean up the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
649
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
650
    else
Kirill Terekhov's avatar
Kirill Terekhov committed
651
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
652
        std::cout << argv[0] << " mesh_file" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
653
    }
Kirill Terekhov's avatar
Kirill Terekhov committed
654
    
Kirill Terekhov's avatar
Kirill Terekhov committed
655
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
656
    Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
657
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
658
659
    Solver::Finalize(); // Finalize solver and close MPI activity
    return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
660
}