main.cpp 10.9 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
#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

19
20
21
22
23
void make_vec(Storage::real v1[3], Storage::real v2[3], Storage::real out[3])
{
	out[0] = v1[0] - v2[0];
	out[1] = v1[1] - v2[1];
	out[2] = v1[2] - v2[2];
Kirill Terekhov's avatar
Kirill Terekhov committed
24
25
}

26
27
28
Storage::real dot_prod(Storage::real v1[3], Storage::real v2[3])
{
	return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
Kirill Terekhov's avatar
Kirill Terekhov committed
29
30
31
}


32
33
Storage::real func(Storage::real x[3], Storage::real tmp)
{
Kirill Terekhov's avatar
Kirill Terekhov committed
34
35
36
37
	//  	return x[0] + 2 * x[1] + 3 * x[2];
	double s0 = sin (M_PI * x[0]); 
	double s1 = sin (M_PI * x[1]);
	double s2 = sin (M_PI * x[2]);
38
39
	return s0 * s1 * s2;
	(void) tmp;
Kirill Terekhov's avatar
Kirill Terekhov committed
40
41
}

42
43
Storage::real func_rhs(Storage::real x[3], Storage::real tmp)
{
Kirill Terekhov's avatar
Kirill Terekhov committed
44
	//  	return 0;
45
	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
46
47
48
}


49
50
51
52

int main(int argc,char ** argv)
{
	Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
Kirill Terekhov's avatar
Kirill Terekhov committed
53
#if defined(USE_PARTITIONER)
54
	Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
55
#endif
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
	if( argc > 1 )
	{
		Tag phi, tensor_K, id;
		Mesh * m = new Mesh(); // Create an empty mesh
		double ttt = Timer();
		bool repartition = false;
		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
76
		BARRIER;
77
78
79
80
81
82
83
84
85
86
		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
87

Kirill Terekhov's avatar
Kirill Terekhov committed
88
89
#if defined(USE_PARTITIONER)
		if (!repartition) 
Kirill Terekhov's avatar
Kirill Terekhov committed
90
		{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
91
92
93
94
95
			ttt = Timer();
			Partitioner * p = new Partitioner(m);
			p->SetMethod(Partitioner::Inner_RCM,Partitioner::Partition); // Specify the partitioner
			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
			delete p;
Kirill Terekhov's avatar
Kirill Terekhov committed
96
			BARRIER;
97
98
99
100
101
102

			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
103
			BARRIER;
104
105
106

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

109
110
		ttt = Timer();
		m->AssignGlobalID(CELL | EDGE | FACE | NODE);
Kirill Terekhov's avatar
Kirill Terekhov committed
111
		BARRIER;
Dmitry Bagaev's avatar
Dmitry Bagaev committed
112
		if( m->GetProcessorRank() == 0 ) std::cout << "Assign id: " << Timer()-ttt << std::endl;
113
114
115
116
		id = m->GlobalIDTag(); // Get the tag of the global ID
		//m->Save("solution_check_0.vtk");
		phi = m->CreateTag("Solution",DATA_REAL,CELL,NONE,1); // Create a new tag for the solution phi
		tensor_K = m->CreateTag("K",DATA_REAL,CELL,NONE,1); // Create a new tag for K tensor
Kirill Terekhov's avatar
Kirill Terekhov committed
117
		//m->Save("solution_check_1.vtk");
118
119
120
121
122
123
124
		for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
			if( cell->GetStatus() != Element::Ghost ) // If the cell is an own one
				cell->Real(tensor_K) = 1.0; // Store the tensor K value into the tag

		ttt = Timer();
		m->ExchangeGhost(1,FACE);
		m->ExchangeData(tensor_K,CELL,0); // Exchange the tensor_K data over processors
Kirill Terekhov's avatar
Kirill Terekhov committed
125
		BARRIER;
126
127
128
129
130
131
		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
132
		S.SetParameter("absolute_tolerance", "1e-8");
Dmitry Bagaev's avatar
Dmitry Bagaev committed
133
134
135
		S.SetParameter("schwartz_overlap", "2");
    	Residual R; // Residual vector
    	Sparse::LockService Locks;
136
137
138
139
140
141
142
143
144
145
146
147
		Sparse::Vector Update; // Declare the solution and the right-hand side vectors

		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
148
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
149
			Automatizator aut;
Kirill Terekhov's avatar
Kirill Terekhov committed
150
			Automatizator::MakeCurrent(&aut);
Kirill Terekhov's avatar
Kirill Terekhov committed
151
152
			INMOST_DATA_ENUM_TYPE iphi = aut.RegisterTag(phi,CELL);
			aut.EnumerateTags();
Kirill Terekhov's avatar
Kirill Terekhov committed
153
154
155
156
157
158
159
160
161

			// Set the indeces intervals for the matrix and vectors
			R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
			R.InitLocks();
			Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
			//~ std::cout << m->GetProcessorRank() << " A,x,b interval " << idmin << ":" << idmax << " size " << idmax-idmin << std::endl;
			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
162
163
164
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
165
166
			{
				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
167
168
169
#if defined(USE_OMP)
#pragma omp for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
170
171
172
173
174
175
176
177
				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
178
					Storage::integer i1 = aut.GetIndex(r1,iphi), i2;
Kirill Terekhov's avatar
Kirill Terekhov committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
					Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1, d2, D, v[3], T;
					Storage::real f_area = face->Area(); // Get the face area
					face->UnitNormal(f_nrm); // Get the face normal
					r1->Centroid(r1_cnt);  // Get the barycenter of the cell
					face->Centroid(f_cnt); // Get the barycenter of the face
					if( !r2->isValid() ) // boundary condition
					{
						Storage::real bnd_pnt[3], dist;
						make_vec(f_cnt,r1_cnt,v);
						dist = dot_prod(f_nrm,v);
						// bnd_pnt is a projection of the cell center to the face
						bnd_pnt[0] = r1_cnt[0] + dist * f_nrm[0];
						bnd_pnt[1] = r1_cnt[1] + dist * f_nrm[1];
						bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
						T = r1->Real(tensor_K) * f_area / dist;
						//flux =  T * (func(bnd_pnt,0) - variable(aut,r1,iphi));
						R.Lock(i1);
						R[i1] -=  T * (func(bnd_pnt,0) - Phi(r1));
						R.UnLock(i1);
					}
					else
					{
Kirill Terekhov's avatar
Kirill Terekhov committed
201
						i2 = aut.GetIndex(r2,iphi);
Kirill Terekhov's avatar
Kirill Terekhov committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
						r2->Centroid(r2_cnt);
						D = dot_prod(f_nrm,f_cnt);
						d1 = fabs(dot_prod(r1_cnt,f_nrm) - D);
						d2 = fabs(dot_prod(r2_cnt,f_nrm) - D);
						T = 1.0 / (d1/r1->Real(tensor_K) + d2/r2->Real(tensor_K)) * f_area;
						//flux = T * (variable(aut,r2,iphi) - variable(aut,r1,iphi));//(unknown(aut,r2,iphi) - unknown(aut,r1,iphi));
						flux = T * (Phi(r2) - Phi(r1));
						if( s1 != Element::Ghost )
						{
							R.Lock(i1);
							R[i1] -= flux;
							R.UnLock(i1);
						}
						if( s2 != Element::Ghost )
						{
							R.Lock(i2);
							R[i2] += flux;
							R.UnLock(i2);
						}
					}
				}
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
224
225
226
#if defined(USE_OMP)
#pragma omp parallel for
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
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
			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[cell->Integer(id)] += cell->Mean(func_rhs, cell->Real(tensor_K)) * cell->Volume();
			}
			BARRIER;
			if( m->GetProcessorRank() == 0 ) std::cout << "Matrix assemble: " << Timer()-ttt << std::endl;

			m->RemoveGeometricData(table); // Clean the computed geometric data

			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
255
				std::cout << S.Residual() << " " << S.Iterations() << " " << S.ReturnReason() << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
256
257
258
259
260
261
262
263
				std::cout << "Solve system: " << Timer()-ttt << std::endl;
			}

			ttt = Timer();

			Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);

			Storage::real err_C = 0.0, err_L2 = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
264
265
266
#if defined(USE_OMP)
#pragma omp parallel
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
267
268
			{
				Storage::real local_err_C = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
269
270
271
#if defined(USE_OMP)
#pragma omp for reduction(+:err_L2)
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
272
273
274
275
276
277
278
				for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell )
				{
					Cell cell = Cell(m,ComposeCellHandle(icell));
					if( cell->GetStatus() != Element::Ghost )
					{
						Storage::real old = cell->Real(phi);
						Storage::real exact = cell->Mean(func, 0); // Compute the mean value of the function over the cell
Kirill Terekhov's avatar
Kirill Terekhov committed
279
						Storage::real res = Update[aut.GetIndex(cell->self(),iphi)];
Kirill Terekhov's avatar
Kirill Terekhov committed
280
281
282
283
284
285
286
287
						Storage::real sol = old-res;
						Storage::real err = fabs (sol - exact);
						if (err > local_err_C) local_err_C = err;
						err_L2 += err * err * cell->Volume();
						cell->Real(error) = err;
						cell->Real(phi) = sol;
					}
				}
Kirill Terekhov's avatar
Kirill Terekhov committed
288
289
290
#if defined(USE_OMP)
#pragma omp critical
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
291
292
293
294
295
296
297
298
299
300
				{
					if( local_err_C > err_C ) err_C = local_err_C;
				}
			}
			err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
			err_L2 = sqrt(m->Integrate(err_L2)); // Compute the global L2 norm for the error
			if( m->GetProcessorRank() == 0 ) std::cout << "err_C  = " << err_C << std::endl;
			if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
		}
		BARRIER;
301
302
303
304
		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
305
		BARRIER;
306
307
308
309
310
311
312
313
314
		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
315
316
		m->Save("result.pmf");
		BARRIER;
317
318
319
320
321
322
323
324
325
		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
326
327

#if defined(USE_PARTITIONER)
328
	Partitioner::Finalize(); // Finalize the partitioner activity
Kirill Terekhov's avatar
Kirill Terekhov committed
329
#endif
330
331
332
	Solver::Finalize(); // Finalize solver and close MPI activity
	return 0;
}