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

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

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

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

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

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

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

142
            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
143
            {
144
145
146
147
148
149
150
151
152
153
154
155
156
                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
157

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

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

164
            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
165
166
            {
                srand(1); // Randomization
167
168
                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
169
170
                    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
171

172
173
174
175
            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
176
177
                    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
178
179


180
181

            if( m->HaveTag("BOUNDARY_CONDITION") ) //Is there boundary condition on the mesh?
Kirill Terekhov's avatar
Kirill Terekhov committed
182
183
184
185
            {
                tag_BC = m->GetTag("BOUNDARY_CONDITION");
                //initialize unknowns at boundary
            }
186
187
188
            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
189
190
            ttt = Timer();
            //Assemble gradient matrix W on cells
Kirill Terekhov's avatar
Kirill Terekhov committed
191
            int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
192
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
193
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
194
#endif
195
196
197
198
199
200
            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

ahmadabushaikha's avatar
ahmadabushaikha committed
201
202
203
				int NF = (int)faces.size(); //number of faces;	
				rMatrix W(NF,NF);
				
Dmitry Bagaev's avatar
Dmitry Bagaev committed
204
				if( cell->GetGeometricType() == Element::Tet && rt0 ) // RT0 consturction of W matrix
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
				{
					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;

ahmadabushaikha's avatar
ahmadabushaikha committed
222
223
					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 
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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
					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;

ahmadabushaikha's avatar
ahmadabushaikha committed
249
250
251
252
253
254
255
256
257
258
259
260
261
262
					// 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; 
Dmitry Bagaev's avatar
Dmitry Bagaev committed
263

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
					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
281

282
283
284
					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
285

286
287
288
					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
289

290
291
292
					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
293

294
295
296
297
298
299
300
301
302
					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
303

304
305
306
307
308
309
						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
310

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

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

320
					}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
321

322
323
324
325
326
327
328
329
					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
330

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

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

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

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



347
348
349
								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
350

351
352
								W_RT0[3]=  2*(-1+x);
								W_RT0[4]=  2*y;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
353

354
								W_RT0[5]=  2*z;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
355

356
357
358
359
360
								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;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
361

362
								W_RT0[10]= 2*y;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
363

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





370
371
372
								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
373
374
375



376
377
378
								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
379

380
381
382
383
384
385
386

								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);


							}

						}
ahmadabushaikha's avatar
ahmadabushaikha committed
387
					}	
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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
					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;
					 */
Kirill Terekhov's avatar
Kirill Terekhov committed
468
                
ahmadabushaikha's avatar
ahmadabushaikha committed
469
					int rank = 0; //size of matrix U
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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
					{ //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
        { //Main loop for problem solution
Kirill Terekhov's avatar
Kirill Terekhov committed
550
            Automatizator aut; // declare class to help manage unknowns
Kirill Terekhov's avatar
Kirill Terekhov committed
551
            Automatizator::MakeCurrent(&aut);
Kirill Terekhov's avatar
Kirill Terekhov committed
552
553
            dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE)); //register pressure as primary unknown
            aut.EnumerateTags(); //enumerate all primary variables
Dmitry Bagaev's avatar
Dmitry Bagaev committed
554
555
            std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;

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

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

582
583
            do
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
584
585
                R.Clear(); //clean up the residual
                double tttt = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
586
                int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
587
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
588
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
Kirill Terekhov's avatar
Kirill Terekhov committed
589
#endif
590
591
592
593
594
595
596
597
598
599
600
601
602
603
                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)
					{
ahmadabushaikha's avatar
ahmadabushaikha committed
604
605
606
607
					
					if (faces[k]->FrontCell().isValid())
					{
						Cell cell_n = cell.Neighbour(faces[k]); 
Dmitry Bagaev's avatar
Dmitry Bagaev committed
608

609
610
611
612
613
						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;

ahmadabushaikha's avatar
ahmadabushaikha committed
614
                        L_1 =0.0; 
Dmitry Bagaev's avatar
Dmitry Bagaev committed
615

616
617
618
619
620
621
622
						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 ;
						}
Dmitry Bagaev's avatar
Dmitry Bagaev committed
623

ahmadabushaikha's avatar
ahmadabushaikha committed
624
						L_2 =0.0; 
Dmitry Bagaev's avatar
Dmitry Bagaev committed
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 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
}