main.cpp 28.6 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
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
typedef Storage::real_array real_array;
typedef Storage::var_array var_array;

Dmitry Bagaev's avatar
Dmitry Bagaev committed
27
bool rt0 = true;
Kirill Terekhov's avatar
Kirill Terekhov committed
28

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

31
32
33
int main(int argc,char ** argv)
{
    Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
34
#if defined(USE_PARTITIONER)
35
    Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
36
#endif
37
38
    if( argc > 1 )
    {
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;
         */
57
58
59
60
61
62
63
		if( argc > 2 )
		{
			if( std::string(argv[2]) == "MFD" )
				rt0 = false;
			else if ( std::string(argv[2]) == "MHFE" )
				rt0 = true;
		}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
64

Kirill Terekhov's avatar
Kirill Terekhov committed
65
66
        double ttt; // Variable used to measure timing
        bool repartition = false; // Is it required to redistribute the mesh?
67
        Mesh * m = new Mesh(); // Create an empty mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
68
69
70
        { // Load the mesh
            ttt = Timer();
            m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
71
72
            if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
            if( m->isParallelFileFormat(argv[1]) ) //The format is
Kirill Terekhov's avatar
Kirill Terekhov committed
73
74
75
            {
                m->Load(argv[1]); // Load mesh from the parallel file format
                repartition = true; // Ask to repartition the mesh
76
77
            }
            else if( m->GetProcessorRank() == 0 ) m->Load(argv[1]); // Load mesh from the serial file format
Kirill Terekhov's avatar
Kirill Terekhov committed
78
            BARRIER
79
            if( m->GetProcessorRank() == 0 ) std::cout << "Load the mesh: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
80
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
81

82
83
84
85
		if( rt0 )
			std::cout << "Running MHFE RT0" << std::endl;
		else
			std::cout << "Running MFD" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
86
87


Kirill Terekhov's avatar
Kirill Terekhov committed
88
#if defined(USE_PARTITIONER)
89
        if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
Kirill Terekhov's avatar
Kirill Terekhov committed
90
        {
Kirill Terekhov's avatar
Kirill Terekhov committed
91
92
93
            { // Compute mesh partitioning
                ttt = Timer();
                Partitioner p(m); //Create Partitioning object
94
                p.SetMethod(Partitioner::Inner_RCM,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
Kirill Terekhov's avatar
Kirill Terekhov committed
95
96
                p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
                BARRIER
97
                if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
98
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
99

Kirill Terekhov's avatar
Kirill Terekhov committed
100
101
102
            { //Distribute the mesh
                ttt = Timer();
                m->Redistribute(); // Redistribute the mesh data
103
                m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
Kirill Terekhov's avatar
Kirill Terekhov committed
104
                BARRIER
105
                if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
106
107
            }
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
108
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
109

Kirill Terekhov's avatar
Kirill Terekhov committed
110
111
112
        { // prepare geometrical data on the mesh
            ttt = Timer();
            Mesh::GeomParam table;
113
114
            table[CENTROID]    = CELL | FACE; //Compute averaged center of mass
            table[NORMAL]      = FACE;        //Compute normals
Kirill Terekhov's avatar
Kirill Terekhov committed
115
            table[ORIENTATION] = FACE;        //Check and fix normal orientation
116
            table[MEASURE]     = CELL | FACE; //Compute volumes and areas
Kirill Terekhov's avatar
Kirill Terekhov committed
117
118
119
            //table[BARYCENTER]  = CELL | FACE; //Compute volumetric center of mass
            m->PrepareGeometricData(table); //Ask to precompute the data
            BARRIER
120
            if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
121
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
122

Kirill Terekhov's avatar
Kirill Terekhov committed
123
124
125
126
127
128
        // 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
129
        Tag tag_DMP; // Indicates weather local W matrix satisfy DMP condition
Dmitry Bagaev's avatar
Dmitry Bagaev committed
130

131
        if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
Kirill Terekhov's avatar
Kirill Terekhov committed
132
133
        { // Exchange ghost cells
            ttt = Timer();
134
            m->ExchangeGhost(1,FACE); // Produce layer of ghost cells
Kirill Terekhov's avatar
Kirill Terekhov committed
135
            BARRIER
136
            if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
137
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
138

Kirill Terekhov's avatar
Kirill Terekhov committed
139
        { //initialize data
140
            if( m->HaveTag("PERM") ) // is diffusion tensor already defined on the mesh? (PERM from permeability)
Kirill Terekhov's avatar
Kirill Terekhov committed
141
                tag_K = m->GetTag("PERM"); // get the diffusion tensor
Dmitry Bagaev's avatar
Dmitry Bagaev committed
142

143
            if( !tag_K.isValid() || !tag_K.isDefined(CELL) ) // diffusion tensor was not initialized or was not defined on cells.
Kirill Terekhov's avatar
Kirill Terekhov committed
144
            {
145
146
147
148
149
150
151
152
153
154
155
156
157
                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
                }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
158

159
                m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
Kirill Terekhov's avatar
Kirill Terekhov committed
160
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
161

162
            if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
Kirill Terekhov's avatar
Kirill Terekhov committed
163
                tag_P = m->GetTag("PRESSURE"); //Get the pressure
Dmitry Bagaev's avatar
Dmitry Bagaev committed
164

165
            if( !tag_P.isValid() || !tag_P.isDefined(CELL) ) // Pressure was not initialized or was not defined on nodes
Kirill Terekhov's avatar
Kirill Terekhov committed
166
167
            {
                srand(1); // Randomization
168
169
                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
Kirill Terekhov's avatar
Kirill Terekhov committed
170
171
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
172

173
174
175
176
            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
Kirill Terekhov's avatar
Kirill Terekhov committed
177
178
                    e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
179
180


181
182

            if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh?
Kirill Terekhov's avatar
Kirill Terekhov committed
183
184
185
186
            {
                tag_BC = m->GetTag("BOUNDARY_CONDITION");
                //initialize unknowns at boundary
            }
187
188
189
            m->ExchangeData(tag_P,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
            tag_W = m->CreateTag("nKGRAD",DATA_REAL,CELL,NONE);
            tag_DMP = m->CreateTag("isDMP",DATA_BULK,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
190
191
            ttt = Timer();
            //Assemble gradient matrix W on cells
Kirill Terekhov's avatar
Kirill Terekhov committed
192
            int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
193
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
194
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
195
#endif
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
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
            {
				Mesh * mesh = m;
                Cell cell = m->CellByLocalID(q);
				ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
				int NF = (int)faces.size(); //number of faces;
				rMatrix W(NF,NF);

				if( cell->GetGeometricType() == Element::Tet && rt0 ) // RT0 consturction of W matrix
				{
					double V = cell.Volume();
					double dN[12];
					double J[9];
					double J_T[9];
					double K_inv[9];
					double K_inv_ref1[9];
					double W_RT0[12];
					double B_RT0[16];
					double gauss_pt_xyz[12];
					double gauss_wei[4];
					double tetra[12];
					double xyz_i[3];
					double xyz_j[3];
					int i, j,k,Z;
					double w,l,m, J_det,x,y,z;


					rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
					rMatrix K_inverse (3,3); // inverse of permeability
					K_inverse= K.Invert(true).first;
					// gauss points for intgration
					gauss_pt_xyz[0]=(5.0+3.0*(sqrt(5.0)))/20.0;
					gauss_pt_xyz[1]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[2]=(5.0-(sqrt(5.0)))/20.0;

					gauss_pt_xyz[3]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[4]=(5.0+3.0*(sqrt(5.0)))/20.0;
					gauss_pt_xyz[5]=(5.0-(sqrt(5.0)))/20.0;

					gauss_pt_xyz[6]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[7]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[8]=(5.0+3.0*(sqrt(5.0)))/20.0;

					gauss_pt_xyz[9]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[10]=(5.0-(sqrt(5.0)))/20.0;
					gauss_pt_xyz[11]=(5.0-(sqrt(5.0)))/20.0;

					// gauss weights
					gauss_wei[0]=  1.0/24.0;
					gauss_wei[1]=  1.0/24.0;
					gauss_wei[2]=  1.0/24.0;
					gauss_wei[3]= 1.0/24.0;

					// derivative of the finite element shape functions for tetraherdon (local)
					dN[0]=-1;
					dN[1]=-1;
					dN[2]=-1;

					dN[3]=1;
					dN[4]=0;
					dN[5]=0;

					dN[6]=0;
					dN[7]=1;
					dN[8]=0;

					dN[9]=0;
					dN[10]=0;
					dN[11]=1;


					// local to reference tranformation
					ElementArray<Node> node(mesh);
					for(j = 0; j < 4; ++j)
					{
						ElementArray<Node> cell_nodes = cell->getNodes();
						cell_nodes.Subtract(faces[j].getNodes());
						node.Unite(cell_nodes);
					}

					for ( j = 0; j < 4; j++ )
					{

						for ( i = 0; i < 3; i++ )
							tetra[i+j*3]= node[j].Coords()[i]; // <-- X
					}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
282

283
284
285
					J[0]=dN[0]*tetra[0]+dN[3]*tetra[3]+dN[6]*tetra[6]+dN[9]*tetra[9];
					J[1]=dN[0]*tetra[1]+dN[3]*tetra[4]+dN[6]*tetra[7]+dN[9]*tetra[10];
					J[2]=dN[0]*tetra[2]+dN[3]*tetra[5]+dN[6]*tetra[8]+dN[9]*tetra[11];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
286

287
288
289
					J[3]=dN[1]*tetra[0]+dN[4]*tetra[3]+dN[7]*tetra[6]+dN[10]*tetra[9];
					J[4]=dN[1]*tetra[1]+dN[4]*tetra[4]+dN[7]*tetra[7]+dN[10]*tetra[10];
					J[5]=dN[1]*tetra[2]+dN[4]*tetra[5]+dN[7]*tetra[8]+dN[10]*tetra[11];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
290

291
292
293
					J[6]=dN[2]*tetra[0]+dN[5]*tetra[3]+dN[8]*tetra[6]+dN[11]*tetra[9];
					J[7]=dN[2]*tetra[1]+dN[5]*tetra[4]+dN[8]*tetra[7]+dN[11]*tetra[10];
					J[8]=dN[2]*tetra[2]+dN[5]*tetra[5]+dN[8]*tetra[8]+dN[11]*tetra[11];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
294

295
296
297
298
299
300
301
302
303
					w = J[0] * (J[4] * J[8] -J [5] * J[7]);
					l = J[1] * (J[3] * J[8] -J [5] * J[6]);
					m = J[2] * (J[3] * J[7] -J [4] * J[6]);
					J_det = w - l + m;
					// cout << " Determinent equals " << J_det << endl;
					if (J_det == 0.0)
					{
						std::cout << "As Determinent=0 so it is singular matrix and its inverse cannot exist for element "
							<< std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
304

305
306
307
308
309
310
						exit(0);
					}
					// transpose of jacobian
					J_T[0]=J[0]; J_T[1]=J[3]; J_T[2]=J[6];
					J_T[3]=J[1]; J_T[4]=J[4]; J_T[5]=J[7];
					J_T[6]=J[2]; J_T[7]=J[5]; J_T[8]=J[8];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
311

312
313
314
					// inverse of K tensor
					for ( j = 0; j < 3; j++ )
					{
Dmitry Bagaev's avatar
Dmitry Bagaev committed
315

316
317
318
319
						for ( i = 0; i < 3; i++ )
						{
							K_inv[i+j*3]=K_inverse (i,j);
						}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
320

321
					}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
322

323
324
325
326
327
328
329
330
					for ( i = 0; i < 3; i++) { // row number of output
						for ( j = 0; j < 3; j++) { // column number of output
							K_inv_ref1[3*i+j] = 0.0;
							for ( k = 0; k < 3; k++) { // three elements are added for this output
								K_inv_ref1[3*i+j] += J[3*i+k] * K_inv[3*k+j];
							}
						}
					}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
331

332
333
334
335
					for ( i = 0; i < 16; i++)
					{
						B_RT0[i]=0.0;
					}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
336

337
					for ( i = 0; i < 4; i++) { // row number of output
Dmitry Bagaev's avatar
Dmitry Bagaev committed
338

339
						for ( j = 0; j < 4; j++) { // column number of output
Dmitry Bagaev's avatar
Dmitry Bagaev committed
340
341


342
343
344
345
							for ( Z = 0; Z < 4; Z++) { // GAUSS POINT IN Z DIRECTION
								x=gauss_pt_xyz[Z*3];
								y=gauss_pt_xyz[Z*3+1];
								z=gauss_pt_xyz[Z*3+2];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
346
347
348



349
350
351
								W_RT0[0]=  2*x;///(sqrt(3.0));
								W_RT0[1]=  2*y;///(sqrt(3.0));
								W_RT0[2]=  2*z;///(sqrt(3.0));
Dmitry Bagaev's avatar
Dmitry Bagaev committed
352

353
354
355
								W_RT0[3]=  2*(-1+x);
								W_RT0[4]=  2*y;
								W_RT0[5]=  2*z;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
356

357
358
359
360
361
362
363
								W_RT0[6]=  (2*x);///(sqrt(3.0));
								W_RT0[7]=  2*(-1+y);;///(sqrt(3.0));
								W_RT0[8]=  (2*z);///(sqrt(3.0));

								W_RT0[9]=  2*x;
								W_RT0[10]= 2*y;
								W_RT0[11]= 2*(-1+z);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
364
365
366
367
368





369
370
371
								xyz_i[0]=W_RT0[i*3]*J[0]+W_RT0[i*3+1]*J[3]+W_RT0[i*3+2]*J[6];
								xyz_i[1]=W_RT0[i*3]*J[1]+W_RT0[i*3+1]*J[4]+W_RT0[i*3+2]*J[7];
								xyz_i[2]=W_RT0[i*3]*J[2]+W_RT0[i*3+1]*J[5]+W_RT0[i*3+2]*J[8];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
372
373
374



375
376
377
								xyz_j[0]=W_RT0[j*3]*K_inv_ref1[0]+W_RT0[j*3+1]*K_inv_ref1[3]+W_RT0[j*3+2]*K_inv_ref1[6];
								xyz_j[1]=W_RT0[j*3]*K_inv_ref1[1]+W_RT0[j*3+1]*K_inv_ref1[4]+W_RT0[j*3+2]*K_inv_ref1[7];
								xyz_j[2]=W_RT0[j*3]*K_inv_ref1[2]+W_RT0[j*3+1]*K_inv_ref1[5]+W_RT0[j*3+2]*K_inv_ref1[8];
Dmitry Bagaev's avatar
Dmitry Bagaev committed
378

379
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
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
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
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532

								B_RT0[i+j*4]+=gauss_wei[Z]*(xyz_i[0]*xyz_j[0]+ xyz_i[1]*xyz_j[1]+xyz_i[2]*xyz_j[2])/std::abs(J_det);


							}

						}
					}
					for ( j = 0; j < 4; j++) { // row number of output
						double sum = 0;
						for ( i = 0; i < 4; i++) { // column number of output

							sum += W(i,j);
							W(i,j)= B_RT0[i+j*4];
						}
						//std::cout << "row " << j << " sum " << sum << std::endl;
					}


					//std::cout << "W" << std::endl;
					//W.Print();

					W = W.Invert(true).first ;


					for ( j = 0; j < 4; j++) { // row number of output
						double sum = 0;
						for ( i = 0; i < 4; i++) { // column number of output

							sum += W(i,j);
						}
						//std::cout << "row inverted " << j << " sum " << sum << std::endl;
					}

					//std::cout << "W inverted" << std::endl;

					//W.Print();



					//scanf("%*c");

				}
				else
				{
					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);
					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 NK(NF,3), R(NF,3), D(NF,NF), U(NF,NF), Areas(NF,NF); //big gradient matrix, co-normals, directions
					Areas.Zero();
					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);

						Areas(k,k) = aF;
					} //end of loop over faces
					rMatrix SU,SS,SV;
					W = 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());

					//std::cout << "UDtR" << std::endl;
					//(U*D.Transpose()*R).Print();

					W += D*U*D.Transpose();


					W = Areas*W*Areas;
				}
					//std::cout << "W: " << std::endl;
					//W.Print();
				bulk & isDMP = cell->Bulk(tag_DMP);
				isDMP = 1;
				for(int k = 0; k < NF; ++k)
				{
					real row_sum = 0;
					if( W(k,k) < 0.0 ) isDMP = 0;
					for(int j = 0; j < NF; ++j)
						row_sum += W(k,j);
					if( row_sum < 0.0 ) isDMP = 0;
					for(int j = k+1; j < NF; ++j)
						if( W(k,j) > 0.0 )
							isDMP = 0;
				}
				++total;
				if( isDMP ) ++dmp;
				real_array store_W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
				store_W.resize(NF*NF); //resize the structure
				std::copy(W.data(),W.data()+NF*NF,store_W.data()); //write down the gradient matrix
            } //end of loop over cells
Kirill Terekhov's avatar
Kirill Terekhov committed
533
            std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
534
            std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
535

536
            if( m->HaveTag("FORCE") ) //Is there force on the mesh?
Kirill Terekhov's avatar
Kirill Terekhov committed
537
538
539
540
541
            {
                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
Dmitry Bagaev's avatar
Dmitry Bagaev committed
542

Kirill Terekhov's avatar
Kirill Terekhov committed
543
        std::cout << "Initialization done" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
544
545


Kirill Terekhov's avatar
Kirill Terekhov committed
546
547
        integer nit = 0;
        ttt = Timer();
Dmitry Bagaev's avatar
Dmitry Bagaev committed
548

Kirill Terekhov's avatar
Kirill Terekhov committed
549
550
551
        { //Main loop for problem solution
            Automatizator aut(m); // declare class to help manage unknowns
            Automatizator::MakeCurrent(&aut);
552
            dynamic_variable P(aut,aut.RegisterDynamicTag(tag_P,CELL|FACE)); //register pressure as primary unknown
Kirill Terekhov's avatar
Kirill Terekhov committed
553
            aut.EnumerateDynamicTags(); //enumerate all primary variables
Dmitry Bagaev's avatar
Dmitry Bagaev committed
554
555
556
557

            std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;


558
559
560
561
            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
Kirill Terekhov's avatar
Kirill Terekhov committed
562
            {//Annotate matrix
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
                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
578
                    }
579
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
580
            }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
581

Kirill Terekhov's avatar
Kirill Terekhov committed
582
            std::cout << "Matrix was annotated" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
583

584
585
            do
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
586
587
                R.Clear(); //clean up the residual
                double tttt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
588
                int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
589
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
590
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
591
#endif
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
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
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
                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
                    bulk & isDMP = cell->Bulk(tag_DMP);
                    vMatrix pF(NF,1); //vector of pressure differences on faces
                    vMatrix FLUX(NF,1); //computed flux on faces
                    vMatrix FLUX_MASS(NF,1); //computed flux on faces for mass balance on cell
					vMatrix POT (2,1); //vector of pressure differences on faces

					for(int k = 0; k < NF; ++k)
					{


					if (faces[k]->FrontCell().isValid())
					{
						Cell cell_n = cell.Neighbour(faces[k]);
						ElementArray<Face> faces_n = cell_n->getFaces(); //obtain faces of the cell neighoubr
					    rMatrix nKGRAD_n(cell_n->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
						double B_1, B_2, L_1, L_2;
						int face_n_k;

                        L_1 =0.0;
						B_1= nKGRAD(k,k) ;
						for(int j = 0; j < NF; ++j)
						L_1 +=nKGRAD(k,j) ;
						for(int j = 0; j < NF; ++j) if (faces_n[j] == faces[k])
						{
							face_n_k= j ;
						}
						L_2 =0.0;
						for(int j = 0; j < NF; ++j)
						L_2 +=nKGRAD_n(face_n_k,j) ;

						B_2 = nKGRAD_n(face_n_k,face_n_k) ;

						POT(0,0)=L_1*P(cell) ;
						POT(1,0)=L_2*P(cell_n) ;

						for(int j = 0; j < NF; ++j)if (j != k)
						POT(0,0) -= nKGRAD(k,j)* P(faces[j]);

						for(int j = 0; j < NF; ++j) if (j != face_n_k)
						POT(1,0) -= nKGRAD_n(face_n_k,j)* P(faces_n[j]);

                        FLUX_MASS(k,0) = -(B_2*POT(0,0)-B_1 *POT(1,0))/(B_1+B_2) ;


					}
					else {
						FLUX_MASS(k,0)=0.0;

					for(int j = 0; j < NF; ++j)
						FLUX_MASS(k,0)+= nKGRAD(k,j)* (P(faces[j]) - P(cell)) ;
					}
					}


				//	for(int k = 0; k < NF; ++k)
				//	{
				//	FLUX_MASS(k,0)=0.0;
				////	Cell cell_l =faces[k].BackCell();
				//	Cell cell_l= cell.Neighbour(faces[k]);
				//	Cell cell_n= cell_l.Neighbour(faces[k]);
				//	for(int j = 0; j < NF; ++j)
				//	FLUX_MASS(k,0)+= nKGRAD(k,j)* (P(faces[j]) - P(cell_n)) ;
				//	}

	                    for(int k = 0; k < NF; ++k)
                    pF(k,0) = (P(faces[k]) - P(cell));//*faces[k].Area();
                    FLUX = nKGRAD*pF; //fluxes on faces
                    //Change matrix by nonlinear correction for DMP
                    ++total;
                    if( cell.GetStatus() != Element::Ghost )
Kirill Terekhov's avatar
Kirill Terekhov committed
668
                    {
669
670
671
672
673
674
675
676
677
678
679
                        for(int k = 0; k < NF; ++k) //loop over faces of current cell
                        { //R[P.Index(cell)] += FLUX(k,0);//faces[k].Area();
						   R[P.Index(cell)] +=  FLUX_MASS(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) )
Kirill Terekhov's avatar
Kirill Terekhov committed
680
                        {
681
682
                            real_array BC = faces[k].RealArray(tag_BC);
                            R[index] -= BC[0]*P(faces[k]) + BC[1]*FLUX(k,0) - BC[2];
Kirill Terekhov's avatar
Kirill Terekhov committed
683
                        }
684
685
686
687
688
                        else
                            R[index] -= FLUX(k,0);
                        Locks.UnLock(index);
                    }
                } //end of loop over cells
Dmitry Bagaev's avatar
Dmitry Bagaev committed
689

Kirill Terekhov's avatar
Kirill Terekhov committed
690
                std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
691
692


693
694
                if( tag_F.isValid() )
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
695
696
697
#if defined(USE_OMP)
#pragma omp parallel for
#endif
698
699
700
701
702
703
                    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();
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
704
                }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
705

Kirill Terekhov's avatar
Kirill Terekhov committed
706
                std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
707
708


Kirill Terekhov's avatar
Kirill Terekhov committed
709
                R.Rescale();
Kirill Terekhov's avatar
Kirill Terekhov committed
710
711
                //R.GetJacobian().Save("jacobian.mtx",&Text);
                //R.GetResidual().Save("residual.mtx");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
712
713


Kirill Terekhov's avatar
Kirill Terekhov committed
714
                std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
715

716
                if( R.Norm() < 1.0e-4 ) break;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
717

718
				//Solver S(Solver::INNER_ILU2);
719
                //Solver S(Solver::INNER_MPTILUC);
720
				Solver S("superlu");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
721
722
723
724
                S.SetParameter("relative_tolerance", "1.0e-14");
                S.SetParameter("absolute_tolerance", "1.0e-12");
                S.SetParameter("drop_tolerance", "1.0e-1");
                S.SetParameter("reuse_tolerance", "1.0e-2");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
725

Kirill Terekhov's avatar
Kirill Terekhov committed
726
                S.SetMatrix(R.GetJacobian());
Kirill Terekhov's avatar
Kirill Terekhov committed
727
                //std::fill(Update.Begin(),Update.End(),0.0);
728
729
                if( S.Solve(R.GetResidual(),Update) )
                {
Kirill Terekhov's avatar
Kirill Terekhov committed
730
731
732
#if defined(USE_OMP)
#pragma omp parallel for
#endif
733
734
735
736
737
                    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
738
739
740
#if defined(USE_OMP)
#pragma omp parallel for
#endif
741
742
743
744
745
746
                    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);
Kirill Terekhov's avatar
Kirill Terekhov committed
747
748
749
                    {
                        std::stringstream str;
                        str << "iter" << nit;
750
                        if( m->GetProcessorsNumber() == 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
751
752
753
754
755
                            str << ".vtk";
                        else
                            str << ".pvtk";
                        m->Save(str.str());
                    }
756
757
758
                }
                else
                {
Dmitry Bagaev's avatar
Dmitry Bagaev committed
759
                    std::cout << "Unable to solve: " << S.ReturnReason() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
760
761
762
                    break;
                }
                ++nit;
763
            } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm
Kirill Terekhov's avatar
Kirill Terekhov committed
764
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
765
        std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
766

767
768
769
        if( m->HaveTag("REFERENCE_SOLUTION") )
        {
            Tag tag_E = m->CreateTag("ERRROR",DATA_REAL,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
770
771
772
            Tag tag_R = m->GetTag("REFERENCE_SOLUTION");
            real C, L2, volume;
            C = L2 = volume = 0.0;
773
774
775
776
777
778
779
780
781
782
783
            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);
Dmitry Bagaev's avatar
Dmitry Bagaev committed
784
785
            std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
            C = L2 = volume = 0.0;
786
787
788
789
790
791
792
793
794
795
796
797
798
799
            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);
Kirill Terekhov's avatar
Kirill Terekhov committed
800
                std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
801
802
            }
            else std::cout << "Reference solution was not defined on faces" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
803
        }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
804

805
        if( m->GetProcessorsNumber() == 1 )
Kirill Terekhov's avatar
Kirill Terekhov committed
806
807
808
            m->Save("out.vtk");
        else
            m->Save("out.pvtk");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
809

Kirill Terekhov's avatar
Kirill Terekhov committed
810
        delete m; //clean up the mesh
811
812
813
    }
    else
    {
Kirill Terekhov's avatar
Kirill Terekhov committed
814
        std::cout << argv[0] << " mesh_file" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
815
    }
Dmitry Bagaev's avatar
Dmitry Bagaev committed
816

Kirill Terekhov's avatar
Kirill Terekhov committed
817
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
818
    Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
819
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
820
821
    Solver::Finalize(); // Finalize solver and close MPI activity
    return 0;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
822
}