diffusion.cpp 22 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
#include "inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
2
3
4
5
6
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

Kirill Terekhov's avatar
Kirill Terekhov committed
7

Kirill Terekhov's avatar
Kirill Terekhov committed
8
9
10
11
12
13
14
15
16
17
18
19
20
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
21
typedef Storage::bulk bulk;
Kirill Terekhov's avatar
Kirill Terekhov committed
22
23
24
25
26
27
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
28
bool print_niter = false; //save file on nonlinear iterations
29
bool output_matrix = false;
30
31
32
33
34
35
36
37
38
39
40
41
bool norne_wells = true;

//data for wells
Cell cc[3] = {InvalidCell(),InvalidCell(),InvalidCell()};
real WI[3] = {50000,50000,50000};
real pbhp[3] = {265,105,110};
real ccnt[3][3] =
{
	{4.567151e+05, 7.321079e+06, 2.767665e+03},
	{4.609346e+05, 7.323503e+06, 2.597767e+03},
	{4.595400e+05, 7.326078e+06, 2.803586e+03}
};
Kirill Terekhov's avatar
Kirill Terekhov committed
42
43


Kirill Terekhov's avatar
Kirill Terekhov committed
44
//#define OPTIMIZATION
Kirill Terekhov's avatar
Kirill Terekhov committed
45

Kirill Terekhov's avatar
Kirill Terekhov committed
46
int main(int argc,char ** argv)
Kirill Terekhov's avatar
Kirill Terekhov committed
47
{
Kirill Terekhov's avatar
Kirill Terekhov committed
48
49
50
51
52
53
    Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
#if defined(USE_PARTITIONER)
    Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
#endif
    if( argc > 1 )
    {
54
55
		std::string force_name = "FORCE";
		int force_comp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
56
57
58
59
60
        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();
Kirill Terekhov's avatar
Kirill Terekhov committed
61
62
63
            //~ m->SetParallelFileStrategy(0);
			//~ m->SetParallelStrategy(1);
            //~ m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
Kirill Terekhov's avatar
Kirill Terekhov committed
64
65
66
67
68
69
70
71
72
73
            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
74

75
76
		if( argc > 2 ) force_name = std::string(argv[2]);
		if( argc > 3 ) force_comp = atoi(argv[3]);
Kirill Terekhov's avatar
Kirill Terekhov committed
77
78
79
80
81
82
83
84
85
86
87
88
		
#if defined(USE_PARTITIONER)
        if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
        {
            { // Compute mesh partitioning
                ttt = Timer();
                Partitioner p(m); //Create Partitioning object
                p.SetMethod(Partitioner::INNER_KMEANS,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;
            }
Kirill Terekhov's avatar
Kirill Terekhov committed
89

Kirill Terekhov's avatar
Kirill Terekhov committed
90
91
92
93
94
95
96
97
98
            { //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;
            }
        }
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
99

Kirill Terekhov's avatar
Kirill Terekhov committed
100
101
102
103
104
105
106
107
108
109
110
111
        //{ // 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
112

Kirill Terekhov's avatar
Kirill Terekhov committed
113
114
115
116
117
        // data tags for
        Tag tag_P;  // Pressure
        Tag tag_K;  // Diffusion tensor
        Tag tag_F;  // Forcing term
        Tag tag_W;  // Gradient matrix acting on harmonic points on faces and returning gradient on faces
118
        TagRealArray tag_BC; // Boundary conditions
119
        
Kirill Terekhov's avatar
Kirill Terekhov committed
120
		
121
122
123
		
		TagRealArray tag_PG; // Pressure gradient
		TagRealArray tag_WG; // matrix to reconstruct gradient
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
        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;
        }

        { //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.
            {
Kirill Terekhov's avatar
Kirill Terekhov committed
139
                tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,1); // create a new tag for symmetric diffusion tensor K
Kirill Terekhov's avatar
Kirill Terekhov committed
140
141
142
143
144
145
                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
Kirill Terekhov's avatar
Kirill Terekhov committed
146
147
148
149
150
                    //~ 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
Kirill Terekhov's avatar
Kirill Terekhov committed
151
152
153
154
                }

                m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
            }
155
156
157
158
159
160
161
162
163
164
165
166
            
            if( norne_wells )
            {
				for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) if( it->GetStatus() != Element::Ghost )
				{
					for(int q = 0; q < 3; ++q) if( it->Inside(ccnt[q]) )
					{
						cc[q] = it->self();
						std::cout << "proc " << m->GetProcessorRank() << " found c" << q << " " << cc[q].LocalID() << std::endl;
					}
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
167
168
169

            if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
                tag_P = m->GetTag("PRESSURE"); //Get the pressure
170
171
172
                
            tag_PG = m->CreateTag("PRESSURE_GRADIENT",DATA_REAL,CELL,NONE,3);
            tag_WG = m->CreateTag("WGRAD",DATA_REAL,CELL,NONE);
Kirill Terekhov's avatar
Kirill Terekhov committed
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197

            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
198
			
Kirill Terekhov's avatar
Kirill Terekhov committed
199
200
201
202
            ttt = Timer();
            //Assemble gradient matrix W on cells
#if defined(USE_OMP)
#pragma omp parallel
Kirill Terekhov's avatar
Kirill Terekhov committed
203
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
204
			{
205
				rMatrix NK, L, R, Areas;
Kirill Terekhov's avatar
Kirill Terekhov committed
206
				rMatrix x(1,3), xf(1,3), n(1,3);
207
				double area, dist; //area of the face
Kirill Terekhov's avatar
Kirill Terekhov committed
208
209
#if defined(USE_OMP)
#pragma omp for
Kirill Terekhov's avatar
Kirill Terekhov committed
210
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
211
212
213
214
215
216
				 for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
				 {
					 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);
217
					 cell->Barycenter(x.data());
Kirill Terekhov's avatar
Kirill Terekhov committed
218
219
220
					 //get permeability for the cell
					 rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),
													 cell->RealArrayDF(tag_K).size());
Kirill Terekhov's avatar
Kirill Terekhov committed
221
222
					 NK.Resize(NF,3); //co-normals
					 R.Resize(NF,3); //directions
223
					 L.Resize(NF,NF);
Kirill Terekhov's avatar
Kirill Terekhov committed
224
					 Areas.Resize(NF,NF); //areas
Kirill Terekhov's avatar
Kirill Terekhov committed
225
					 Areas.Zero();
226
					 L.Zero();
Kirill Terekhov's avatar
Kirill Terekhov committed
227
228
229
					 for(int k = 0; k < NF; ++k) //loop over faces
					 {
						 area = faces[k].Area();
230
						 faces[k].Barycenter(xf.data());
Kirill Terekhov's avatar
Kirill Terekhov committed
231
						 faces[k].OrientedUnitNormal(cell->self(),n.data());
232
						 dist = n.DotProduct(xf-x);
Kirill Terekhov's avatar
Kirill Terekhov committed
233
						 // assemble matrix of directions
234
						 R(k,k+1,0,3) = (xf-x);
Kirill Terekhov's avatar
Kirill Terekhov committed
235
						 // assemble matrix of co-normals
236
237
						 NK(k,k+1,0,3) = n*K*area;
						 L(k,k) =  n.DotProduct(n*K)*area/dist;
Kirill Terekhov's avatar
Kirill Terekhov committed
238
239
						 Areas(k,k) = area;
					 } //end of loop over faces
240
241
242
243
					 //~ W = NK*(NK.Transpose()*R).PseudoInvert(1.0e-12)*NK.Transpose(); //stability part
					 //~ W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())*
						//~ (2.0/(static_cast<real>(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace());
			 		 
244
245
246
			 		 
			 		 W = (NK)*((NK).Transpose()*R).Invert()*(NK).Transpose();
					 W+= L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
247
248
					 
					 //~ W = Areas*W*Areas;
Kirill Terekhov's avatar
Kirill Terekhov committed
249
250
251
252
253
254
					 //access data structure for gradient matrix in mesh
					 real_array store_W = cell->RealArrayDV(tag_W);
					 //resize the structure
					 store_W.resize(NF*NF);
					 //write down the gradient matrix
					 std::copy(W.data(),W.data()+NF*NF,store_W.data());
255
256
					 
					 tag_WG[cell].resize(3*NF);
257
					 tag_WG(cell,3,NF) = (NK.Transpose()*R).PseudoInvert(1.0e-12)*(NK).Transpose();
Kirill Terekhov's avatar
Kirill Terekhov committed
258
				 } //end of loop over cells
Kirill Terekhov's avatar
Kirill Terekhov committed
259
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
260
261
            std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
			 
262
            if( m->HaveTag(force_name) ) //Is there force on the mesh?
Kirill Terekhov's avatar
Kirill Terekhov committed
263
            {
264
                tag_F = m->GetTag(force_name); //initial force
Kirill Terekhov's avatar
Kirill Terekhov committed
265
266
267
268
269
270
271
272
273
274
275
276
277
278
                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; // declare class to help manage unknowns
            Automatizator::MakeCurrent(&aut);
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
            
            MarkerType unk = m->CreateMarker();
            for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
				m->CellByLocalID(q).SetMarker(unk);
			for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
			{
				Face face = m->FaceByLocalID(q);
				face.SetMarker(unk);
				if( face.Boundary() )
				{
					double alpha = 0, beta = 1, gamma = 0;
					if( tag_BC.isValid() && face.HaveData(tag_BC) )
					{
						alpha = tag_BC[face][0];
						beta  = tag_BC[face][1];
						gamma = tag_BC[face][2];
					}
					if( alpha != 0 && beta == 0 )
					{
						face.RemMarker(unk);
						face.Real(tag_P) = gamma/alpha;
					}
				}
			}
            
            dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE,unk)); //register pressure as primary unknown
Kirill Terekhov's avatar
Kirill Terekhov committed
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
            aut.EnumerateEntries(); //enumerate all primary variables
            std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << 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);
322
                    if( face.GetStatus() != Element::Ghost && face.GetMarker(unk) )
Kirill Terekhov's avatar
Kirill Terekhov committed
323
                    {
324
325
326
327
328
329
330
331
332
333
334
335
336
                        if( face.Boundary() )
                        {
							double alpha = 0, beta = 1, gamma = 0;
							if( tag_BC.isValid() && face.HaveData(tag_BC) )
							{
								alpha = tag_BC[face][0];
								beta  = tag_BC[face][1];
								gamma = tag_BC[face][2];
							}
							std::stringstream str;
							str << "Pressure guided by boundary condition " << alpha << " " << beta << " " << gamma;
                            Text.SetAnnotation(P.Index(face),str.str());
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
337
338
339
340
341
342
343
344
345
                        else
                            Text.SetAnnotation(P.Index(face),"Interface pressure");
                    }
                }
            }

            std::cout << "Matrix was annotated" << std::endl;
			 
            do
Kirill Terekhov's avatar
Kirill Terekhov committed
346
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
347
348
                R.Clear(); //clean up the residual
                double tttt = Timer();
349
                //~ int total = 0, dmp = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
350
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
351
#pragma omp parallel
Kirill Terekhov's avatar
Kirill Terekhov committed
352
#endif
353
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
354
355
					vMatrix pF; //vector of pressure differences on faces
					vMatrix FLUX; //computed flux on faces
356
					rMatrix n(3,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
357
358
359
#if defined(USE_OMP)
#pragma omp for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
360
					for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
Dmitry Bagaev's avatar
Dmitry Bagaev committed
361
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
362
363
364
						Cell cell = m->CellByLocalID(q);
						ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
						int NF = (int)faces.size();
Kirill Terekhov's avatar
Kirill Terekhov committed
365
						
Kirill Terekhov's avatar
Kirill Terekhov committed
366
						raMatrix W = raMatrixMake(cell->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
Kirill Terekhov's avatar
Kirill Terekhov committed
367
						
Kirill Terekhov's avatar
Kirill Terekhov committed
368
369
						pF.Resize(NF,1);
						FLUX.Resize(NF,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
370
371
						
						for(int k = 0; k < NF; ++k)
372
							pF(k,0) = (P[faces[k]] - P[cell]);
Kirill Terekhov's avatar
Kirill Terekhov committed
373
						FLUX = W*pF; //fluxes on faces
374
						tag_PG(cell,3,1) = tag_WG(cell,3,NF)*pF;
Kirill Terekhov's avatar
Kirill Terekhov committed
375
						if( cell.GetStatus() != Element::Ghost )
376
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
377
378
							for(int k = 0; k < NF; ++k) //loop over faces of current cell
								R[P.Index(cell)] += FLUX(k,0);//faces[k].Area();
379
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
380
						for(int k = 0; k < NF; ++k) //loop over faces of current cell
381
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
382
							if( faces[k].GetStatus() == Element::Ghost ) continue;
383
							if( !faces[k].GetMarker(unk) ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
384
385
							int index = P.Index(faces[k]);
							Locks.Lock(index);
386
387
							
							if( faces[k].Boundary() )
Kirill Terekhov's avatar
Kirill Terekhov committed
388
							{
389
390
391
392
393
394
395
396
397
								double a = 0, b = 1, c = 0;
								//~ faces[k].UnitNormal(n.data());
								//~ double a = 1, b = 0, c = 0;
								//~ if( fabs(n(2,0)-1) < 1.0e-4 )
								//~ {
									//~ a = 0;
									//~ b = 1;
								//~ }
								
398
399
400
401
402
403
								if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
								{
									real_array BC = faces[k].RealArray(tag_BC);
									a = BC[0], b = BC[1], c = BC[2];
								}
								R[index] -= a*P(faces[k]) + b*FLUX(k,0) - c;
Kirill Terekhov's avatar
Kirill Terekhov committed
404
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
405
406
407
							else
								R[index] -= FLUX(k,0);
							Locks.UnLock(index);
408
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
409
410
411
412
413
					} //end of loop over cells


					 if( tag_F.isValid() )
					 {
Kirill Terekhov's avatar
Kirill Terekhov committed
414
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
415
#pragma omp for
Kirill Terekhov's avatar
Kirill Terekhov committed
416
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
417
418
419
420
						 for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
						 {
							 Cell cell = m->CellByLocalID(q);
							 if( cell.GetStatus() == Element::Ghost ) continue;
421
							 if( cell->HaveData(tag_F) ) R[P.Index(cell)] += cell->RealArray(tag_F)[force_comp]*cell->Volume();
Kirill Terekhov's avatar
Kirill Terekhov committed
422
423
						 }
					 }
424
					 
425
426
427
428
429
430
					 if( norne_wells )
					 {
						 for(int q = 0; q < 3; ++q) if( cc[q].isValid() )
							R[P.Index(cc[q])] += WI[q]*(pbhp[q] - P[cc[q]])*cc[q].Volume();
					 }
					 
431
					 //R[P.Index(m->BeginCell()->self())] = P[m->BeginCell()->self()];
Kirill Terekhov's avatar
Kirill Terekhov committed
432
433
				}
				std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
434
435
436
437
438

                std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;

                if( R.Norm() < 1.0e-4 ) break;

439
440
				//~ R.GetJacobian().Save("A.mtx");
				//~ R.GetResidual().Save("b.mtx");
441
				//Solver S(Solver::INNER_ILU2);
442
443
444
                Solver S(Solver::INNER_MPTILUC);
                //~ Solver S(Solver::INNER_MLMPTILUC);
                //~ Solver S(Solver::K3BIILU2);
Kirill Terekhov's avatar
Kirill Terekhov committed
445
				//Solver S("superlu");
446
447
				S.SetParameter("verbosity","3");
				S.SetParameter("rescale_iterations", "10");
Kirill Terekhov's avatar
Kirill Terekhov committed
448
449
                S.SetParameter("relative_tolerance", "1.0e-14");
                S.SetParameter("absolute_tolerance", "1.0e-12");
450
451
452
                //~ S.SetParameter("drop_tolerance", "5.0e-2");
                //~ S.SetParameter("reuse_tolerance", "2.5e-3");
                S.SetParameter("drop_tolerance", "1.0e-2");
453
                S.SetParameter("reuse_tolerance", "1.0e-4");
454
455
                S.SetParameter("pivot_condition", "20");
                double tset = Timer(), titr;
Kirill Terekhov's avatar
Kirill Terekhov committed
456
457

                S.SetMatrix(R.GetJacobian());
458
459
460
461
462
463
464
465
466
467
468
469
470
                
                tset = Timer() - tset;
                
				if( output_matrix )
				{
					std::cout << "write A.mtx" << std::endl;
					R.GetJacobian().Save("A.mtx");//,&Text);
					std::cout << "write b.mtx" << std::endl;
					R.GetResidual().Save("b.mtx");
					std::cout << "write done, solve" << std::endl;
				}
				
				titr = Timer();
Kirill Terekhov's avatar
Kirill Terekhov committed
471
				
472
473
474
475
476
				bool success =  S.Solve(R.GetResidual(),Update);
				
				titr = Timer() - titr;
				
                if( success )
Kirill Terekhov's avatar
Kirill Terekhov committed
477
                {
478
479
480
481
482
483
484
485
					std::cout << "Solved system in " << S.Iterations() << " iterations, time for setup " << tset << "s iterations " << titr << "s" << std::endl;
					if( output_matrix )
					{
						std::cout << "write x.mtx" << std::endl;
						Update.Save("x.mtx");
						//std::cout << "exit" << std::endl;
						//exit(-1);
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
486
487
488
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
489
490
491
492
                    for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
                    {
                        Cell cell = m->CellByLocalID(q);
						if( cell->GetStatus() == Element::Ghost ) continue;
493
						if( !cell->GetMarker(unk) ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
494
495
                        cell->Real(tag_P) -= Update[P.Index(cell)];
                    }
Kirill Terekhov's avatar
Kirill Terekhov committed
496
497
498
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
499
500
501
502
                    for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
                    {
                        Face face = m->FaceByLocalID(q);
						if( face->GetStatus() == Element::Ghost ) continue;
503
						if( !face->GetMarker(unk) ) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
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
                        face->Real(tag_P) -= Update[P.Index(face)];
                    }
                    m->ExchangeData(tag_P, CELL|FACE, 0);
					
					if( print_niter )
                    {
                        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.ReturnReason() << std::endl;
                    break;
                }
                ++nit;
            } while( R.Norm() > 1.0e-4 && nit < 10); //check the residual norm
        }
        std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;

        if( m->HaveTag("REFERENCE_SOLUTION") )
        {
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
531
            Tag tag_E = m->CreateTag("ERROR",DATA_REAL,CELL,NONE,1);
Kirill Terekhov's avatar
Kirill Terekhov committed
532
533
534
535
536
537
538
539
540
541
542
543
544
            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;
            }
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
545
546
547
			C = m->AggregateMax(C);
			L2 = m->Integrate(L2);
			volume = m->Integrate(volume);
Kirill Terekhov's avatar
Kirill Terekhov committed
548
            L2 = sqrt(L2/volume);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
549
            if( m->GetProcessorRank() == 0 ) std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
550
551
552
553
554
555
556
557
558
559
560
561
562
563
            C = L2 = volume = 0.0;
            if( tag_R.isDefined(FACE) )
            {
                tag_E = m->CreateTag("ERROR",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;
                }
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
564
565
566
				C = m->AggregateMax(C);
				L2 = m->Integrate(L2);
				volume = m->Integrate(volume);
Kirill Terekhov's avatar
Kirill Terekhov committed
567
                L2 = sqrt(L2/volume);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
568
                if( m->GetProcessorRank() == 0 ) std::cout << "Error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
569
570
571
            }
            else std::cout << "Reference solution was not defined on faces" << std::endl;
        }
Kirill Terekhov's avatar
Kirill Terekhov committed
572
573
574
575
576
577
        
        tag_W = m->DeleteTag(tag_W);
        tag_WG = m->DeleteTag(tag_WG);
        tag_K = m->DeleteTag(tag_K);
        tag_BC = m->DeleteTag(tag_BC);
        m->RemoveGeometricData(table);
Kirill Terekhov's avatar
Kirill Terekhov committed
578
579
580
581
582

        if( m->GetProcessorsNumber() == 1 )
            m->Save("out.vtk");
        else
            m->Save("out.pvtk");
583
584
            
        m->Save("solution.pmf");
Kirill Terekhov's avatar
Kirill Terekhov committed
585
586
587
588
589
590
591
592

        delete m; //clean up the mesh
    }
    else
    {
        std::cout << argv[0] << " mesh_file" << std::endl;
    }

Kirill Terekhov's avatar
Kirill Terekhov committed
593
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
Kirill Terekhov committed
594
    Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
595
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
596
597
    Solver::Finalize(); // Finalize solver and close MPI activity
    return 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
598
}