main.cpp 12.1 KB
Newer Older
1
#include "inmost.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.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

Kirill Terekhov's avatar
Kirill Terekhov committed
19
double func(double x[3], double tmp)
20
{
Kirill Terekhov's avatar
Kirill Terekhov committed
21
	//  	return x[0] + 2 * x[1] + 3 * x[2];
Kirill Terekhov's avatar
Kirill Terekhov committed
22
	return sin (M_PI * x[0]) * sin (M_PI * x[1]) * sin (M_PI * x[2]);
23
	(void) tmp;
Kirill Terekhov's avatar
Kirill Terekhov committed
24
25
}

Kirill Terekhov's avatar
Kirill Terekhov committed
26
27

double func_rhs(double x[3], double tmp)
28
{
Kirill Terekhov's avatar
Kirill Terekhov committed
29
	//  	return 0;
30
	return -3 * tmp * M_PI * M_PI * sin (M_PI * x[0]) * sin (M_PI * x[1]) * sin (M_PI * x[2]);
Kirill Terekhov's avatar
Kirill Terekhov committed
31
32
33
}


34
35
36
37

int main(int argc,char ** argv)
{
	Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
38
#if defined(USE_PARTITIONER)
39
	Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
40
#endif
41
42
	if( argc > 1 )
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
43
44
45
46
47
		TagReal phi;
		TagReal tag_F;
		TagRealArray tag_K;
		TagRealArray tag_BC;
		TagReal phi_ref;
48
49
50
		Mesh * m = new Mesh(); // Create an empty mesh
		double ttt = Timer();
		bool repartition = false;
51
		(void)repartition;
52
53
54
55
56
57
58
59
60
61
62
63
64
65
		m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
		if( m->GetProcessorRank() == 0 ) // If the current process is the master one
			std::cout << argv[0] << std::endl;

		if( m->isParallelFileFormat(argv[1]) )
		{
			m->Load(argv[1]); // Load mesh from the parallel file format
			repartition = true;
		}
		else
		{
			if( m->GetProcessorRank() == 0 )
				m->Load(argv[1]); // Load mesh from the serial file format
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
66
		BARRIER;
67
68
69
70
71
72
73
74
75
76
		if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
		if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_File): " << Timer()-ttt << std::endl;

		//~ double ttt2 = Timer();
		//~ Mesh t;
		//~ t.SetCommunicator(INMOST_MPI_COMM_WORLD);
		//~ t.SetParallelFileStrategy(0);
		//~ t.Load(argv[1]);
		//~ BARRIER
		//~ if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_Scatter): " << Timer()-ttt2 << std::endl;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
77

Kirill Terekhov's avatar
Kirill Terekhov committed
78
#if defined(USE_PARTITIONER)
Kirill Terekhov's avatar
update    
Kirill Terekhov committed
79
		if (m->GetProcessorsNumber() > 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
80
		{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
81
82
			ttt = Timer();
			Partitioner * p = new Partitioner(m);
83
			p->SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition); // Specify the partitioner
84
85
			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
			delete p;
Kirill Terekhov's avatar
Kirill Terekhov committed
86
			BARRIER;
87
88
89
90
91
92

			if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;

			ttt = Timer();
			m->Redistribute(); // Redistribute the mesh data
			m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
Kirill Terekhov's avatar
Kirill Terekhov committed
93
			BARRIER;
94
95
96

			if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
97
#endif
Dmitry Bagaev's avatar
Dmitry Bagaev committed
98

99
100
		ttt = Timer();
		phi = m->CreateTag("Solution",DATA_REAL,CELL,NONE,1); // Create a new tag for the solution phi
Kirill Terekhov's avatar
Kirill Terekhov committed
101
102
103
104
105
106
107
		
		bool makerefsol = true;
		
		if( m->HaveTag("PERM" ) )
		{
			tag_K = m->GetTag("PERM");
			makerefsol = false;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
108
			std::cout << "Permeability from grid" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
109
110
111
		}
		else
		{
Kirill Terekhov's avatar
update    
Kirill Terekhov committed
112
			std::cout << "Set perm" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
113
114
115
116
117
			tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,1); // Create a new tag for K tensor
			for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
				tag_K[*cell][0] = 1.0; // Store the tensor K value into the tag
		}
		
Kirill Terekhov's avatar
Kirill Terekhov committed
118
		
Kirill Terekhov's avatar
Kirill Terekhov committed
119
120
121
122
123
		
		if( m->HaveTag("BOUNDARY_CONDITION") )
		{
			tag_BC = m->GetTag("BOUNDARY_CONDITION");
			makerefsol = false;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
124
			std::cout << "Boundary conditions from grid" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
125
126
127
		}
		else
		{
Kirill Terekhov's avatar
update    
Kirill Terekhov committed
128
			std::cout << "Set boundary conditions" << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
129
			double x[3];
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
130
			tag_BC = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
Kirill Terekhov's avatar
Kirill Terekhov committed
131
			for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
Kirill Terekhov's avatar
update    
Kirill Terekhov committed
132
				if( face->Boundary() && !(face->GetStatus() == Element::Ghost) )
Kirill Terekhov's avatar
Kirill Terekhov committed
133
134
135
136
137
138
				{
					face->Centroid(x);
					tag_BC[*face][0] = 1; //dirichlet
					tag_BC[*face][1] = 0; //neumann
					tag_BC[*face][2] = func(x,0);//face->Mean(func, 0);
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
139
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
140
		
Kirill Terekhov's avatar
Kirill Terekhov committed
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
		if( m->HaveTag("FORCE") )
		{
			tag_F = m->GetTag("FORCE");
			makerefsol = false;
			std::cout << "Force from grid" << std::endl;
		}
		else if( makerefsol )
		{
			std::cout << "Set rhs" << std::endl;
			tag_F = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,1); // Create a new tag for external force
			double x[3];
			for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
			{
				cell->Centroid(x);
				tag_F[*cell] = -func_rhs(x,1);
				//tag_F[*cell] = -cell->Mean(func_rhs,1);
			}
		}
		
Kirill Terekhov's avatar
Kirill Terekhov committed
160
161
162
163
164
165
166
167
168
169
170
171
		if(m->HaveTag("REFERENCE_SOLUTION") )
			phi_ref = m->GetTag("REFERENCE_SOLUTION");
		else if( makerefsol )
		{
			phi_ref = m->CreateTag("REFRENCE_SOLUTION",DATA_REAL,CELL,NONE,1);
			double x[3];
			for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell )
			{
				cell->Centroid(x);
				phi_ref[*cell] = func(x,0);//cell->Mean(func, 0);
			}
		}
172
173
174

		ttt = Timer();
		m->ExchangeGhost(1,FACE);
Kirill Terekhov's avatar
Kirill Terekhov committed
175
		m->ExchangeData(tag_K,CELL,0); // Exchange the tag_K data over processors
Kirill Terekhov's avatar
Kirill Terekhov committed
176
		BARRIER;
177
178
179
180
181
182
		if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;



		ttt = Timer();
		Solver S("inner_ilu2"); // Specify the linear solver to ASM+ILU2+BiCGStab one
Dmitry Bagaev's avatar
Dmitry Bagaev committed
183
		S.SetParameter("absolute_tolerance", "1e-8");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
184
185
186
		S.SetParameter("schwartz_overlap", "2");
    	Residual R; // Residual vector
    	Sparse::LockService Locks;
187
188
		Sparse::Vector Update; // Declare the solution and the right-hand side vectors

Kirill Terekhov's avatar
Kirill Terekhov committed
189
190
191
192
193
194
195
196
197
198
199
200
		{
			Mesh::GeomParam table;
			table[CENTROID] = CELL | FACE;
			table[NORMAL] = FACE;
			table[ORIENTATION] = FACE;
			table[MEASURE] = CELL | FACE;
			table[BARYCENTER] = CELL | FACE;
			m->PrepareGeometricData(table);
		}
		BARRIER
		if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
		
Kirill Terekhov's avatar
Kirill Terekhov committed
201
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
202
			Automatizator aut;
Kirill Terekhov's avatar
Kirill Terekhov committed
203
			Automatizator::MakeCurrent(&aut);
Kirill Terekhov's avatar
Kirill Terekhov committed
204
			INMOST_DATA_ENUM_TYPE iphi = aut.RegisterTag(phi,CELL);
205
			aut.EnumerateEntries();
Kirill Terekhov's avatar
Kirill Terekhov committed
206
207
208
209
210

			// Set the indeces intervals for the matrix and vectors
			R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
			R.InitLocks();
			Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
Kirill Terekhov's avatar
Kirill Terekhov committed
211
			
Kirill Terekhov's avatar
Kirill Terekhov committed
212
213
214
			dynamic_variable Phi(aut,iphi);
			// Solve \nabla \cdot \nabla phi = f equation
			//for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
Kirill Terekhov's avatar
Kirill Terekhov committed
215
216
217
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
218
219
			{
				variable flux; //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
Kirill Terekhov's avatar
Kirill Terekhov committed
220
221
				rMatrix x1(3,1), x2(3,1), xf(3,1), n(3,1);
				double d1, d2, k1, k2, area, T, a, b, c;
Kirill Terekhov's avatar
Kirill Terekhov committed
222
223
224
#if defined(USE_OMP)
#pragma omp for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
225
226
227
228
229
230
231
232
				for(Storage::integer iface = 0; iface < m->FaceLastLocalID(); ++iface ) if( m->isValidFace(iface) )
				{
					Face face = Face(m,ComposeFaceHandle(iface));
					Element::Status s1,s2;
					Cell r1 = face->BackCell();
					Cell r2 = face->FrontCell();
					if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
						((!r2->isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
Kirill Terekhov's avatar
Kirill Terekhov committed
233
234
235
236
237
238
239
240
241
					
					
					area = face->Area(); // Get the face area
					face->UnitNormal(n.data()); // Get the face normal
					face->Centroid(xf.data()); // Get the barycenter of the face
					r1->Centroid(x1.data());  // Get the barycenter of the cell
					k1 = n.DotProduct(rMatrix::FromTensor(tag_K[r1].data(),
														  tag_K[r1].size(),3)*n);
					d1 = fabs(n.DotProduct(xf-x1));
Kirill Terekhov's avatar
Kirill Terekhov committed
242
243
244
					if( !r2->isValid() ) // boundary condition
					{
						// bnd_pnt is a projection of the cell center to the face
Kirill Terekhov's avatar
Kirill Terekhov committed
245
246
247
248
249
						// a*pb + bT(pb-p1) = c
						// F = T(pb-p1)
						// pb = (c + bTp1)/(a+bT)
						// F = T/(a+bT)(c - ap1)
						T = k1/d1;
Kirill Terekhov's avatar
minifix    
Kirill Terekhov committed
250
251
252
						a = 0;
						b = 1;
						c = 0;
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
253
						if( tag_BC.isValid() && face.HaveData(tag_BC) )
Kirill Terekhov's avatar
minifix    
Kirill Terekhov committed
254
255
256
257
						{
							a = tag_BC[face][0];
							b = tag_BC[face][1];
							c = tag_BC[face][2];
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
258
							//std::cout << "a " << a << " b " << b << " c " << c << std::endl;
Kirill Terekhov's avatar
minifix    
Kirill Terekhov committed
259
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
260
261
262
						R.Lock(Phi.Index(r1));
						R[Phi.Index(r1)] -=  T/(a + b*T) * area * (c - a*Phi(r1));
						R.UnLock(Phi.Index(r1));
Kirill Terekhov's avatar
Kirill Terekhov committed
263
264
265
					}
					else
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
266
267
268
269
270
271
						r2->Centroid(x2.data());
						k2 = n.DotProduct(rMatrix::FromTensor(tag_K[r2].data(),
															  tag_K[r2].size(),3)*n);
						d2 = fabs(n.DotProduct(x2-xf));
						T = 1.0/(d1/k1 + d2/k2);
						flux = T* area * (Phi(r2) - Phi(r1));
Kirill Terekhov's avatar
Kirill Terekhov committed
272
273
						if( s1 != Element::Ghost )
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
274
275
276
							R.Lock(Phi.Index(r1));
							R[Phi.Index(r1)] -= flux;
							R.UnLock(Phi.Index(r1));
Kirill Terekhov's avatar
Kirill Terekhov committed
277
278
279
						}
						if( s2 != Element::Ghost )
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
280
281
282
							R.Lock(Phi.Index(r2));
							R[Phi.Index(r2)] += flux;
							R.UnLock(Phi.Index(r2));
Kirill Terekhov's avatar
Kirill Terekhov committed
283
284
285
286
						}
					}
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
287
288
			if( tag_F.isValid() )
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
289
290
291
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
292
293
294
295
296
297
				for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell ) if( m->isValidCell(icell) )
				{
					Cell cell = Cell(m,ComposeCellHandle(icell));
					if( cell->GetStatus() != Element::Ghost )
						R[Phi.Index(cell)] -= tag_F[cell] * cell->Volume();
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
298
299
300
301
			}
			BARRIER;
			if( m->GetProcessorRank() == 0 ) std::cout << "Matrix assemble: " << Timer()-ttt << std::endl;

Kirill Terekhov's avatar
Kirill Terekhov committed
302
			//m->RemoveGeometricData(table); // Clean the computed geometric data
Kirill Terekhov's avatar
Kirill Terekhov committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320

			if( argc > 3 ) // Save the matrix and RHS if required
			{
				ttt = Timer();
				R.GetJacobian().Save(std::string(argv[2])); // "A.mtx"
				R.GetResidual().Save(std::string(argv[3])); // "b.rhs"
				BARRIER;
				if( m->GetProcessorRank() == 0 ) std::cout << "Save matrix \"" << argv[2] << "\" and RHS \"" << argv[3] << "\": " << Timer()-ttt << std::endl;
			}

			ttt = Timer();

			S.SetMatrix(R.GetJacobian()); // Compute the preconditioner for the original matrix
			S.Solve(R.GetResidual(),Update);   // Solve the linear system with the previously computted preconditioner

			BARRIER;
			if( m->GetProcessorRank() == 0 ) 
			{
Dmitry Bagaev's avatar
Dmitry Bagaev committed
321
				std::cout << S.Residual() << " " << S.Iterations() << " " << S.ReturnReason() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
322
323
324
325
326
				std::cout << "Solve system: " << Timer()-ttt << std::endl;
			}

			ttt = Timer();

Kirill Terekhov's avatar
Kirill Terekhov committed
327
328
329
330
			
			if( phi_ref.isValid() )
			{
				Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
331
				double err_C = 0.0, err_L2 = 0.0, vol = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
332
333
334
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
335
336
				{
					double local_err_C = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
337
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
338
#pragma omp for reduction(+:err_L2) reduction(+:vol)
Kirill Terekhov's avatar
Kirill Terekhov committed
339
#endif
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
340
					for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell ) if( m->isValidCell(icell) )
Kirill Terekhov's avatar
Kirill Terekhov committed
341
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
342
343
344
345
346
347
348
349
350
351
						Cell cell = Cell(m,ComposeCellHandle(icell));
						if( cell->GetStatus() != Element::Ghost )
						{
							double old = phi[cell];
							double exact = phi_ref[cell];
							double res = Update[Phi.Index(cell)];
							double sol = old-res;
							double err = fabs (sol - exact);
							if (err > local_err_C) local_err_C = err;
							err_L2 += err * err * cell->Volume();
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
352
							vol += cell->Volume();
Kirill Terekhov's avatar
Kirill Terekhov committed
353
354
355
							cell->Real(error) = err;
							phi[cell] = sol;
						}
Kirill Terekhov's avatar
Kirill Terekhov committed
356
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
357
358
359
#if defined(USE_OMP)
#pragma omp critical
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
360
361
362
					{
						if( local_err_C > err_C ) err_C = local_err_C;
					}
Kirill Terekhov's avatar
Kirill Terekhov committed
363
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
364
				err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
365
				err_L2 = sqrt(m->Integrate(err_L2)/m->Integrate(vol)); // Compute the global L2 norm for the error
Kirill Terekhov's avatar
Kirill Terekhov committed
366
367
				if( m->GetProcessorRank() == 0 ) std::cout << "err_C  = " << err_C << std::endl;
				if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
368
369
370
			}
		}
		BARRIER;
371
372
373
374
		if( m->GetProcessorRank() == 0 ) std::cout << "Compute true residual: " << Timer()-ttt << std::endl;

		ttt = Timer();
		m->ExchangeData(phi,CELL,0); // Data exchange over processors
Kirill Terekhov's avatar
Kirill Terekhov committed
375
		BARRIER;
376
377
378
379
380
381
382
383
384
		if( m->GetProcessorRank() == 0 ) std::cout << "Exchange phi: " << Timer()-ttt << std::endl;

		std::string filename = "result";
		if( m->GetProcessorsNumber() == 1 )
			filename += ".vtk";
		else
			filename += ".pvtk";
		ttt = Timer();
		m->Save(filename);
Kirill Terekhov's avatar
Kirill Terekhov committed
385
386
		m->Save("result.pmf");
		BARRIER;
387
388
389
390
391
392
393
394
395
		if( m->GetProcessorRank() == 0 ) std::cout << "Save \"" << filename << "\": " << Timer()-ttt << std::endl;


		delete m;
	}
	else
	{
		std::cout << argv[0] << " mesh_file [A.mtx b.rhs]" << std::endl;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
396
397

#if defined(USE_PARTITIONER)
398
	Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
399
#endif
400
401
	Solver::Finalize(); // Finalize solver and close MPI activity
	return 0;
402
}