BddcMlSolver.cc 12.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.

12
13
#ifdef HAVE_BDDC_ML

14
15
16
17
18
extern "C" {
#include <bddcml_interface_c.h>
}

#include "parallel/BddcMlSolver.h"
19
#include "parallel/MpiHelper.h"
20
21
22

namespace AMDiS {

23
  void BddcMlSolver::fillPetscMatrix(Matrix<DOFMatrix*> *m)
24
25
  {
    FUNCNAME("BddcMlSolver::fillPetscMatrix()");
26
27

    mat = m;
28
29
30
31
32
33
  }


  void BddcMlSolver::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("BddcMlSolver::fillPetscRhs()");
34
35

    rhsVec = vec;
36
37
38
39
40
41
42
  }


  void BddcMlSolver::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
  {
    FUNCNAME("BddcMlSolver::solvePetscMatrix()");
    
43
44
45
    TEST_EXIT(rhsVec)("Should not happen!\n");
    TEST_EXIT(mat)("Should not happen!\n");

46
    int nComponents = vec.getSize();
47
    const FiniteElemSpace *feSpace = vec.getFeSpace(0);
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    Mesh *mesh = feSpace->getMesh();
    

    // === First, create a continous leaf element index on each subdomain ===

    std::set<int> leafElIndex;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      leafElIndex.insert(elInfo->getElement()->getIndex());
      elInfo = stack.traverseNext(elInfo);
    }

    map<int, int> mapElIndex;
62
    int nLeafEls = 0;
63
64
65
66
67
68
    for (std::set<int>::iterator it = leafElIndex.begin();
	 it != leafElIndex.end(); ++it)
      mapElIndex[*it] = nLeafEls++;



69
70
71
72
73
    int nLevel = 2;
    int nSubdomains[nLevel]; 
    nSubdomains[0] = meshDistributor->getMpiSize();
    nSubdomains[1] = 1;

74
    int nSubPerProc = 1;
75
    MPI_Fint c2f = MPI_Comm_c2f(meshDistributor->getMpiComm());
76
77
78
    int verboseLevel = 2;
    int numbase = 0;

79
80
81
//     MSG("call to \"bddcml_init\" with the following arguments:\n");
//     MSG("  %d, [%d, %d], %d, %d, %d, %d, %d\n", nLevel, nSubdomains[0], nSubdomains[1], nLevel, nSubPerProc, c2f, verboseLevel, numbase);

82
    bddcml_init(&nLevel, nSubdomains, &nLevel, &nSubPerProc, 
83
84
85
86
87
88
		&c2f, &verboseLevel, &numbase);

    // global number of elements
    int nelem = mesh->getNumberOfLeaves();
    mpi::globalAdd(nelem);

89
90
    MSG("nelem = %d\n", nelem);

91
    // global number of nodes
Thomas Witkowski's avatar
Thomas Witkowski committed
92
    int nnod = dofMap[feSpace].nOverallDofs;
93
94
    
    MSG("nnod = %d\n", nnod);
95
96
97
98

    // global number of dofs
    int ndof = nnod * nComponents;

99
100
    MSG("ndof = %d\n", ndof);

101
    // space dimenstion
102
    int ndim = mesh->getDim();
103
104

    // mesh dimension
105
    int meshdim = mesh->getDim();
106
107
108
109
110
111
112

    // global indes of subdomain
    int isub = meshDistributor->getMpiRank();

    // local number of elements
    int nelems = nLeafEls;

113
114
    MSG("nelems = %d\n", nelems);

115
    // local number of nodes
116
    int nnods = feSpace->getAdmin()->getUsedSize();
117
118
119
120

    // local number of dofs
    int ndofs = nnods * nComponents;

121
122
    MSG("local nnods %d     ndofs %d\n", nnods, ndofs);

123
124
    int nVertices = mesh->getGeo(VERTEX);

125
    // Length of array inet
126
    int linet = nelems * nVertices;
127
128
129
130
131
    
    // Local array with indices of nodes on each element
    int inet[linet];
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
132
133
134
      TEST_EXIT_DBG(mapElIndex.count(elInfo->getElement()->getIndex()))
	("Should not happen!\n");

135
      int localElIndex = mapElIndex[elInfo->getElement()->getIndex()];
136
137
      for (int i = 0; i < nVertices; i++)
	inet[localElIndex * nVertices + i] = elInfo->getElement()->getDof(i, 0);
138
139
140
141
142
143
144
      elInfo = stack.traverseNext(elInfo);
    }


    // local array of number of nodes per element
    int nnet[nelems];
    for (int i = 0; i < nelems; i++)
145
      nnet[i] = nVertices;
146
147

    // local array with number of DOFs per node.
148
149
    int nndf[nnods];
    for (int i = 0; i < nnods; i++)
150
151
152
      nndf[i] = nComponents;

    // array of indices of subdomain nodes in global numbering
153
154
    int isngn[nnods];
    for (int i = 0; i < nnods; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
155
      isngn[i] = dofMap[feSpace][i].global; //meshDistributor->mapDofToGlobal(feSpace, i);
156
157

    // array of indices of subdomain variables in global numbering
158
    int isvgvn[ndofs];
159
160
161
    for (int j = 0; j < nnods; j++)
      for (int i = 0; i < nComponents; i++)
	isvgvn[j * nComponents + i] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
162
	  dofMap[feSpace][j].global * nComponents + i;
163
164
165
166
167
168

    // array of indices of subdomain elements in global numbering
    int isegn[nelems];
    int rStartEl, nOverallEl;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nelems, rStartEl, nOverallEl);
169
    MSG("rStartEl = %d\n", rStartEl);
170
171
172
173
174
175
    for (int i = 0; i < nelems; i++)
      isegn[i] = rStartEl + i;



    int lxyz1 = nnods;
176
    int lxyz2 = mesh->getDim();
177
178
179
180
181
182
183
184
185
186
187
188
    // local array with coordinates of nodes
    double xyz[lxyz1 * lxyz2];
    
    {
      DOFVector<WorldVector<double> > coordDofs(feSpace, "tmp");
      mesh->getDofIndexCoords(feSpace, coordDofs);

      for (int i = 0; i < lxyz2; i++)
	for (int j = 0; j < nnods; j++)
	  xyz[i * nnods + j] = coordDofs[j][i];
    }

189
190
191
    
    // === Fill for dirichlet boundary conditions. ===

192
193
194
    // local array of indices denoting dirichlet boundary data
    int ifix[ndofs];
    for (int i = 0; i < ndofs; i++)
195
      ifix[i] = 0;
196
197
198

    // local array of values for dirichlet boundary data
    double fixv[ndofs];
199
200
    for (int i = 0; i < ndofs; i++)
      fixv[i] = 0.0;
201

202
203
204
205
206
207
208
209
210
211
212
213
214
215
    for (int i = 0; i < rhsVec->getSize(); i++) {
      map<DegreeOfFreedom, double>& dcValues = 
	rhsVec->getDOFVector(i)->getDirichletValues();

      for (map<DegreeOfFreedom, double>::iterator it = dcValues.begin();
	   it != dcValues.end(); ++it) {
	int index = it->first * nComponents + i;
	TEST_EXIT_DBG(index < ndofs)("Should not happen!\n");
	ifix[index] = 1;
	fixv[index] = it->second;
      }
    }


216
217
218
219
    // local rhs data
    double rhs[ndofs];
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofvec = *(rhsVec->getDOFVector(i));
220
      for (int j = 0; j < nnods; j++)
221
222
223
224
	rhs[j * nComponents + i] = dofvec[j];
    }

    // Completenes of the rhs vector on subdomains
225
    int is_rhs_complete = 0;
226
227
228
229
230
231
232
233

    // Local array with initial solution guess
    double sol[ndofs];
    for (int i = 0; i < ndofs; i++)
      sol[i] = 0.0;

    // matrix type (set here to unsymmetric)
    int matrixtype = 0;
234
    Parameters::get("parallel->bddcml->matrix type", matrixtype);
235

236
237
238
239
240
241
242
243
244
245
246
247
248
249
    // Non zero structure of matrix
    vector<int> i_sparse;
    vector<int> j_sparse;
    vector<double> a_sparse;

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j])
	  addDofMatrix((*mat)[i][j], 
		       i_sparse, j_sparse, a_sparse, nComponents, i, j);
	

    // Number of non-zero entries in matrix
    int la = i_sparse.size();
250

251
    //    MSG("LOCAL LA = %d\n", la);
252
253
    
    // Matrix is assembled
254
    int is_assembled_int = 0;
255
256
257
258
259
260
261
262
263


    // Users constraints
    double user_constraints = 0.0;
   
    int luser_constraints1 = 0;

    int luser_constraints2 = 0;

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
    /*
    string tmp="";

    MSG("call to \"bddcml_upload_subdomain_data\" with the following arguments (each in one line):\n");
    MSG("  nelem = %d\n", nelem);
    MSG("  nnod = %d\n", nnod);
    MSG("  ndof = %d\n", ndof);
    MSG("  ndim = %d\n", ndim);
    MSG("  meshdim = %d\n", meshdim);
    MSG("  isub = %d\n", isub);
    MSG("  nelems = %d\n", nelems);
    MSG("  nnods = %d\n", nnods);
    MSG("  ndofs = %d\n", ndofs);
    MSG("  inet = [%d, %d, %d, %d, %d, %d]\n", inet[0], inet[1], inet[2], inet[3], inet[4], inet[5]);
    MSG("  linet = %d\n", linet);
    MSG("  nnet = [%d, %d]\n", nnet[0], nnet[1]);
    MSG("  lnnet = %d\n", nelems);
    MSG("  nndf = [%d, %d, %d, %d]\n", nndf[0], nndf[1], nndf[2], nndf[3]);
    MSG("  lnndf = %d\n", nnods);
    MSG("  isngn = [%d, %d, %d, %d]\n", isngn[0], isngn[1], isngn[2], isngn[3]);
    MSG("  lisngn = %d\n", nnods);
    MSG("  isvgvn = [%d, %d, %d, %d]\n", isvgvn[0], isvgvn[1], isvgvn[2], isvgvn[3]);
    MSG("  lisvgvn = %d\n", nnods);
    MSG("  isegn = [%d, %d]\n", isegn[0], isegn[1]);
    MSG("  lisegn = %d\n", nelems);
    MSG("  xyz = [%f, %f, %f, %f, %f, %f, %f, %f]\n", xyz[0], xyz[1], xyz[2], xyz[3], xyz[4], xyz[5], xyz[6], xyz[7]);
    MSG("  lxyz1 = %d\n", lxyz1);
    MSG("  lxyz2 = %d\n", lxyz2);
    MSG("  ifix = [%d, %d, %d, %d]\n", ifix[0], ifix[1], ifix[2], ifix[3]);
    MSG("  lifix = %d\n", ndofs);
    MSG("  fixv = [%f, %f, %f, %f]\n", fixv[0], fixv[1], fixv[2], fixv[3]);
    MSG("  lfixv = %d\n", ndofs);
    MSG("  rhs = [%f, %f, %f, %f]\n", rhs[0], rhs[1], rhs[2], rhs[3]);
    MSG("  lrhs = %d\n", ndofs);
    MSG("  is_rhs_complete = %d\n", is_rhs_complete);
    MSG("  sol = [%f, %f, %f, %f]\n", sol[0], sol[1], sol[2], sol[3]);
    MSG("  lsol = %d\n", ndofs);
    MSG("  matrixtype = %d\n", matrixtype);
    
303

304
305
306
307
308
309
    MSG("  ispare = [%d, %d, %d, %d, %d, %d, %d, %d, %d, %d]\n", i_sparse[0], i_sparse[1], i_sparse[2], i_sparse[3], i_sparse[4], i_sparse[5], i_sparse[6], i_sparse[7], i_sparse[8], i_sparse[9]);
    MSG("  jspare = [%d, %d, %d, %d, %d, %d, %d, %d, %d, %d]\n", j_sparse[0], j_sparse[1], j_sparse[2], j_sparse[3], j_sparse[4], j_sparse[5], j_sparse[6], j_sparse[7], j_sparse[8], j_sparse[9]);
    MSG("  a_spare = [%f, %f, %f, %f, %f, %f, %f, %f, %f, %f]\n", a_sparse[0], a_sparse[1], a_sparse[2], a_sparse[3], a_sparse[4], a_sparse[5], a_sparse[6], a_sparse[7], a_sparse[8], a_sparse[9]);
    MSG("  la = %d\n", la);
    MSG("  is_assembled_int = %d\n", is_assembled_int);
    */
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328

    bddcml_upload_subdomain_data(&nelem,
				 &nnod,
				 &ndof,
				 &ndim,
				 &meshdim,
				 &isub,
				 &nelems,
				 &nnods,
				 &ndofs,
				 inet,
				 &linet,
				 nnet,
				 &nelems,
				 nndf,
				 &nnods,
				 isngn,
				 &nnods,
				 isvgvn,
329
				 &ndofs,
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
				 isegn,
				 &nelems,
				 xyz,
				 &lxyz1,
				 &lxyz2,
				 ifix,
				 &ndofs,
				 fixv,
				 &ndofs,
				 rhs,
				 &ndofs,
				 &is_rhs_complete,
				 sol,
				 &ndofs,
				 &matrixtype,
				 &(i_sparse[0]),
				 &(j_sparse[0]),
				 &(a_sparse[0]),
				 &la,
349
350
351
352
				 &is_assembled_int,
 				 &user_constraints,
 				 &luser_constraints1,
 				 &luser_constraints2);
353
354


355
    int use_defaults_int = 0;
356
    int parallel_division_int = 1;
357
    int use_arithmetic_int = 1;
358
    int use_adaptive_int = 0;
359
360
    int use_user_constraint_int = 0;

361
362
    Parameters::get("parallel->bddcml->arithmetic constraints", use_arithmetic_int);

363
364
365
366
367
368
369
    // MSG("call to \"bddcml_setup_preconditioner\" with the following arguments (each in one line):\n");
//     MSG("  %d\n", matrixtype);
//     MSG("  %d\n", use_defaults_int);
//     MSG("  %d\n", parallel_division_int);
//     MSG("  %d\n", use_arithmetic_int);
//     MSG("  %d\n", use_adaptive_int);

370
371
372
373
    bddcml_setup_preconditioner(&matrixtype, 
				&use_defaults_int,
				&parallel_division_int,
				&use_arithmetic_int,
374
375
				&use_adaptive_int,
				&use_user_constraint_int);
376
377
378

    int method = 1;
    double tol = 1.e-6;
379
    int maxit = 1000;
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
    int ndecrmax = 30;
    int num_iter = 0;
    int converged_reason = 0;
    double condition_number = 0.0;

    bddcml_solve(&c2f, 
		 &method, 
		 &tol,
		 &maxit,
		 &ndecrmax,
		 &num_iter,
		 &converged_reason,
		 &condition_number);

    MSG("BDDCML converged reason: %d within %d iterations \n", 
	converged_reason, num_iter);
396

397
398
399
400
401
402
403
404
    bddcml_download_local_solution(&isub, rhs, &ndofs);

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofvec = *(vec.getDOFVector(i));
      for (int j = 0; j < nnods; j++)
	dofvec[j] = rhs[j * nComponents + i];
    }

405
406
407
408
409
410
411
412
413
    bddcml_finalize();
  }


  void BddcMlSolver::destroyMatrixData()
  {
    FUNCNAME("BddcMlSolver::destroyMatrixData()");
  }

414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
  void BddcMlSolver::destroyVectorData()
  {
    FUNCNAME("BddcMlSolver::destroyVectorData()");
  }


421
  void BddcMlSolver::addDofMatrix(DOFMatrix* dmat, 
422
423
424
				  vector<int>& i_sparse, 
				  vector<int>& j_sparse,
				  vector<double>& a_sparse,
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
				  int nComponents,
				  int ithRowComponent,
				  int ithColComponent)
  {
    FUNCNAME("BddcMlSolver::addDofMatrix()");

    TEST_EXIT(dmat)("Should not happen!\n");

    const FiniteElemSpace *feSpace = dmat->getFeSpace();

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(dmat->getBaseMatrix());
    traits::const_value<Matrix>::type value(dmat->getBaseMatrix());

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    for (cursor_type cursor = begin<row>(dmat->getBaseMatrix()), 
	   cend = end<row>(dmat->getBaseMatrix()); cursor != cend; ++cursor) {
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {	
449
450
	i_sparse.push_back(*cursor * nComponents + ithRowComponent);
	j_sparse.push_back(col(*icursor) * nComponents + ithColComponent);
451
	a_sparse.push_back(value(*icursor));
452
453
454
455
456
      }
    }

  }

457
458
}

459
460
#endif