BddcMlSolver.cc 12.1 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 93
    ParallelDofMapping &dofMap = meshDistributor->getDofMap();
    int nnod = dofMap[feSpace].nOverallDofs;
94 95
    
    MSG("nnod = %d\n", nnod);
96 97 98 99

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

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

102 103 104 105 106 107 108 109 110 111 112 113
    // space dimenstion
    int ndim = 2;

    // mesh dimension
    int meshdim = 2;

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

    // local number of elements
    int nelems = nLeafEls;

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

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

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

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

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

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


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

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

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

    // array of indices of subdomain variables in global numbering
157
    int isvgvn[ndofs];
158 159 160
    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
161
	  dofMap[feSpace][j].global * nComponents + i;
162 163 164 165 166 167

    // array of indices of subdomain elements in global numbering
    int isegn[nelems];
    int rStartEl, nOverallEl;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nelems, rStartEl, nOverallEl);
168
    MSG("rStartEl = %d\n", rStartEl);
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
    for (int i = 0; i < nelems; i++)
      isegn[i] = rStartEl + i;



    int lxyz1 = nnods;
    int lxyz2 = 2;
    // 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];
    }

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

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

    // local rhs data
    double rhs[ndofs];
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofvec = *(rhsVec->getDOFVector(i));
202
      for (int j = 0; j < nnods; j++)
203 204 205 206
	rhs[j * nComponents + i] = dofvec[j];
    }

    // Completenes of the rhs vector on subdomains
207
    int is_rhs_complete = 0;
208 209 210 211 212 213 214 215

    // 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;
216

217 218 219 220 221 222 223 224 225 226 227 228 229 230
    // 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();
231

232
    //    MSG("LOCAL LA = %d\n", la);
233 234
    
    // Matrix is assembled
235
    int is_assembled_int = 0;
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274
    /*
    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);
    
275

276 277 278 279 280 281
    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);
    */
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300

    bddcml_upload_subdomain_data(&nelem,
				 &nnod,
				 &ndof,
				 &ndim,
				 &meshdim,
				 &isub,
				 &nelems,
				 &nnods,
				 &ndofs,
				 inet,
				 &linet,
				 nnet,
				 &nelems,
				 nndf,
				 &nnods,
				 isngn,
				 &nnods,
				 isvgvn,
301
				 &ndofs,
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
				 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,
				 &is_assembled_int);


    int use_defaults_int = 1;
    int parallel_division_int = 1;
    int use_arithmetic_int = 1;
327
    int use_adaptive_int = 0;
328 329 330 331 332 333 334
    // 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);

335 336 337 338 339 340 341 342
    bddcml_setup_preconditioner(&matrixtype, 
				&use_defaults_int,
				&parallel_division_int,
				&use_arithmetic_int,
				&use_adaptive_int);

    int method = 1;
    double tol = 1.e-6;
343
    int maxit = 1000;
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
    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);
360

361 362 363 364 365 366 367 368
    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];
    }

369 370 371 372 373 374 375 376 377
    bddcml_finalize();
  }


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

378

Thomas Witkowski's avatar
Thomas Witkowski committed
379 380 381 382 383 384
  void BddcMlSolver::destroyVectorData()
  {
    FUNCNAME("BddcMlSolver::destroyVectorData()");
  }


385
  void BddcMlSolver::addDofMatrix(DOFMatrix* dmat, 
386 387 388
				  vector<int>& i_sparse, 
				  vector<int>& j_sparse,
				  vector<double>& a_sparse,
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
				  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;

Thomas Witkowski's avatar
Thomas Witkowski committed
409 410
    ParallelDofMapping &dofMap = meshDistributor->getDofMap();

411 412 413
    for (cursor_type cursor = begin<row>(dmat->getBaseMatrix()), 
	   cend = end<row>(dmat->getBaseMatrix()); cursor != cend; ++cursor) {
      int rowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
414
	dofMap[feSpace][*cursor].global * nComponents + ithRowComponent;
415 416 417 418

      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {	
	int colIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
419
	  dofMap[feSpace][col(*icursor)].global * nComponents + ithColComponent;
420 421 422

	double val = value(*icursor);
	
Thomas Witkowski's avatar
Thomas Witkowski committed
423 424
	// 	i_sparse.push_back(rowIndex);
	// 	j_sparse.push_back(colIndex);
425 426
	i_sparse.push_back(*cursor * nComponents + ithRowComponent);
	j_sparse.push_back(col(*icursor) * nComponents + ithColComponent);
427 428 429 430 431 432
	a_sparse.push_back(val);
      }
    }

  }

433 434
}

435 436
#endif