BddcMlSolver.cc 8.66 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
    bddcml_init(&nLevel, nSubdomains, &nLevel, &nSubPerProc, 
80 81 82 83 84 85
		&c2f, &verboseLevel, &numbase);

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

86 87
    MSG("nelem = %d\n", nelem);

88 89
    // global number of nodes
    int nnod = meshDistributor->getNumberOverallDofs(feSpace);
90 91
    
    MSG("nnod = %d\n", nnod);
92 93 94 95

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

96 97
    MSG("ndof = %d\n", ndof);

98 99 100 101 102 103 104 105 106 107 108 109
    // 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;

110 111
    MSG("nelems = %d\n", nelems);

112
    // local number of nodes
113
    int nnods = feSpace->getAdmin()->getUsedSize();
114 115 116 117

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

118 119
    MSG("local nnods %d     ndofs %d\n", nnods, ndofs);

120 121 122 123 124 125 126
    // 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) {
127 128 129
      TEST_EXIT_DBG(mapElIndex.count(elInfo->getElement()->getIndex()))
	("Should not happen!\n");

130 131 132 133 134 135 136 137 138 139 140 141 142
      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.
143 144
    int nndf[nnods];
    for (int i = 0; i < nnods; i++)
145 146 147
      nndf[i] = nComponents;

    // array of indices of subdomain nodes in global numbering
148 149
    int isngn[nnods];
    for (int i = 0; i < nnods; i++)
150
      isngn[i] = meshDistributor->mapDofToGlobal(feSpace, i);
151 152

    // array of indices of subdomain variables in global numbering
153
    int isvgvn[ndofs];
154 155 156
    for (int j = 0; j < nnods; j++)
      for (int i = 0; i < nComponents; i++)
	isvgvn[j * nComponents + i] = 
157
	  meshDistributor->mapDofToGlobal(feSpace, j) * nComponents + i;
158 159 160 161 162 163

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

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

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

    // Completenes of the rhs vector on subdomains
203
    int is_rhs_complete = 0;
204 205 206 207 208 209 210 211

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

213 214 215 216 217 218 219 220 221 222 223 224 225 226
    // 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();
227 228

    MSG("LOCAL LA = %d\n", la);
229 230
    
    // Matrix is assembled
231
    int is_assembled_int = 0;
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251


    bddcml_upload_subdomain_data(&nelem,
				 &nnod,
				 &ndof,
				 &ndim,
				 &meshdim,
				 &isub,
				 &nelems,
				 &nnods,
				 &ndofs,
				 inet,
				 &linet,
				 nnet,
				 &nelems,
				 nndf,
				 &nnods,
				 isngn,
				 &nnods,
				 isvgvn,
252
				 &ndofs,
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
				 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;
278 279
    int use_adaptive_int = 0;
    MSG("BDDC POINT A\n");
280 281 282 283 284
    bddcml_setup_preconditioner(&matrixtype, 
				&use_defaults_int,
				&parallel_division_int,
				&use_arithmetic_int,
				&use_adaptive_int);
285
    MSG("BDDC POINT B\n");
286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303

    int method = 1;
    double tol = 1.e-6;
    int maxit = 1000;
    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);

304 305
    MSG("BDDC POINT C\n");

306 307
    MSG("BDDCML converged reason: %d within %d iterations \n", 
	converged_reason, num_iter);
308 309 310 311 312 313 314 315 316 317

    bddcml_finalize();
  }


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

318 319

  void BddcMlSolver::addDofMatrix(DOFMatrix* dmat, 
320 321 322
				  vector<int>& i_sparse, 
				  vector<int>& j_sparse,
				  vector<double>& a_sparse,
323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
				  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) {
      int rowIndex = 
346
	meshDistributor->mapDofToGlobal(feSpace, *cursor) * nComponents +
347 348 349 350 351
	ithRowComponent;

      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {	
	int colIndex = 
352
	  meshDistributor->mapDofToGlobal(feSpace, col(*icursor)) * nComponents +
353 354 355 356 357 358 359 360 361 362 363 364
	  ithColComponent;

	double val = value(*icursor);
	
	i_sparse.push_back(rowIndex);
	j_sparse.push_back(colIndex);
	a_sparse.push_back(val);
      }
    }

  }

365 366
}

367 368
#endif