1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * functions for the top level MLI data structure
11 *
12 *****************************************************************************/
13
14 #include <string.h>
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include "mli.h"
18 #include "mli_utils.h"
19
20 /*****************************************************************************
21 * constructor
22 *---------------------------------------------------------------------------*/
23
MLI(MPI_Comm comm)24 MLI::MLI( MPI_Comm comm )
25 {
26 #ifdef MLI_DEBUG_DETAILED
27 printf("MLI::MLI\n");
28 #endif
29 mpiComm_ = comm;
30 maxLevels_ = 40;
31 numLevels_ = 40;
32 coarsestLevel_ = 0;
33 outputLevel_ = 0;
34 assembled_ = MLI_FALSE;
35 tolerance_ = 1.0e-6;
36 maxIterations_ = 20;
37 currIter_ = 0;
38 oneLevels_ = new MLI_OneLevel*[maxLevels_];
39 for (int j = 0; j < maxLevels_; j++) oneLevels_[j] = new MLI_OneLevel(this);
40 for (int i = 0; i < maxLevels_; i++)
41 {
42 oneLevels_[i]->setLevelNum(i);
43 if ( i < (maxLevels_-1) ) oneLevels_[i]->setNextLevel(oneLevels_[i+1]);
44 if ( i > 0 ) oneLevels_[i]->setPrevLevel(oneLevels_[i-1]);
45 }
46 coarseSolver_ = NULL;
47 methodPtr_ = NULL;
48 solveTime_ = 0.0;
49 buildTime_ = 0.0;
50 }
51
52 /*****************************************************************************
53 * destructor
54 *---------------------------------------------------------------------------*/
55
~MLI()56 MLI::~MLI()
57 {
58 #ifdef MLI_DEBUG_DETAILED
59 printf("MLI::~MLI\n");
60 #endif
61 for ( int i = 0; i < maxLevels_; i++ ) delete oneLevels_[i];
62 delete [] oneLevels_;
63 if ( coarseSolver_ != NULL ) delete coarseSolver_;
64 if ( methodPtr_ != NULL ) delete methodPtr_;
65 }
66
67 /*****************************************************************************
68 * set discretization matrix
69 *---------------------------------------------------------------------------*/
70
setSystemMatrix(int level,MLI_Matrix * A)71 int MLI::setSystemMatrix( int level, MLI_Matrix *A )
72 {
73 #ifdef MLI_DEBUG_DETAILED
74 printf("MLI::setSystemMatrix, level = %d\n", level);
75 #endif
76 if ( level >= 0 && level < maxLevels_ ) oneLevels_[level]->setAmat( A );
77 else
78 {
79 printf("MLI::setSystemMatrix ERROR : wrong level = %d\n", level);
80 exit(1);
81 }
82 return 0;
83 }
84
85 /*****************************************************************************
86 * set restriction operator
87 *---------------------------------------------------------------------------*/
88
setRestriction(int level,MLI_Matrix * R)89 int MLI::setRestriction( int level, MLI_Matrix *R )
90 {
91 #ifdef MLI_DEBUG_DETAILED
92 printf("MLI::setRestriction, level = %d\n", level);
93 #endif
94 if ( level >= 0 && level < maxLevels_ ) oneLevels_[level]->setRmat( R );
95 else
96 {
97 printf("MLI::setRestriction ERROR : wrong level = %d\n", level);
98 exit(1);
99 }
100 return 0;
101 }
102
103 /*****************************************************************************
104 * set prolongation operator
105 *---------------------------------------------------------------------------*/
106
setProlongation(int level,MLI_Matrix * P)107 int MLI::setProlongation( int level, MLI_Matrix *P )
108 {
109 #ifdef MLI_DEBUG_DETAILED
110 printf("MLI::setProlongation, level = %d\n", level);
111 #endif
112 if ( level >= 0 && level < maxLevels_ ) oneLevels_[level]->setPmat( P );
113 else
114 {
115 printf("MLI::setProlongation ERROR : wrong level = %d\n", level);
116 exit(1);
117 }
118 return 0;
119 }
120
121 /*****************************************************************************
122 * set smoother
123 *---------------------------------------------------------------------------*/
124
setSmoother(int level,int pre_post,MLI_Solver * smoother)125 int MLI::setSmoother( int level, int pre_post, MLI_Solver *smoother )
126 {
127 #ifdef MLI_DEBUG_DETAILED
128 printf("MLI::setSmoother, level = %d\n", level);
129 #endif
130 if ( level >= 0 && level < maxLevels_ )
131 {
132 oneLevels_[level]->setSmoother( pre_post, smoother );
133 }
134 else
135 {
136 printf("MLI::setSmoother ERROR : wrong level = %d\n", level);
137 exit(1);
138 }
139 return 0;
140 }
141
142 /*****************************************************************************
143 * set coarse solver
144 *---------------------------------------------------------------------------*/
145
setCoarseSolve(MLI_Solver * solver)146 int MLI::setCoarseSolve( MLI_Solver *solver )
147 {
148 #ifdef MLI_DEBUG_DETAILED
149 printf("MLI::setCoarseSolve\n");
150 #endif
151 if ( ! assembled_ ) coarseSolver_ = solver;
152 else oneLevels_[coarsestLevel_]->setCoarseSolve(solver);
153 return 0;
154 }
155
156 /*****************************************************************************
157 * set finite element data information
158 *---------------------------------------------------------------------------*/
159
setFEData(int level,MLI_FEData * fedata,MLI_Mapper * map)160 int MLI::setFEData( int level, MLI_FEData *fedata, MLI_Mapper *map )
161 {
162 #ifdef MLI_DEBUG_DETAILED
163 printf("MLI::setFEData\n");
164 #endif
165 if ( level >= 0 && level < maxLevels_ )
166 {
167 oneLevels_[level]->setFEData( fedata, map );
168 }
169 else
170 {
171 printf("MLI::setFEData ERROR : wrong level = %d\n", level);
172 exit(1);
173 }
174 return 0;
175 }
176
177 /*****************************************************************************
178 * set finite element data information
179 *---------------------------------------------------------------------------*/
180
setSFEI(int level,MLI_SFEI * sfei)181 int MLI::setSFEI( int level, MLI_SFEI *sfei )
182 {
183 #ifdef MLI_DEBUG_DETAILED
184 printf("MLI::setSFEI\n");
185 #endif
186 if ( level >= 0 && level < maxLevels_ )
187 {
188 oneLevels_[level]->setSFEI( sfei );
189 }
190 else
191 {
192 printf("MLI::setSFEI ERROR : wrong level = %d\n", level);
193 exit(1);
194 }
195 return 0;
196 }
197
198 /*****************************************************************************
199 * set cycle type at various levels
200 *---------------------------------------------------------------------------*/
201
setCyclesAtLevel(int level,int cycles)202 int MLI::setCyclesAtLevel( int level, int cycles )
203 {
204 #ifdef MLI_DEBUG_DETAILED
205 printf("MLI::setCyclesAtLevel at level %d, cycles = %d\n",level,cycles);
206 #endif
207 if ( level >= 0 && level < maxLevels_ )
208 {
209 oneLevels_[level]->setCycles( cycles );
210 }
211 else if ( level == -1 )
212 {
213 for (int i = 0; i < maxLevels_; i++) oneLevels_[i]->setCycles( cycles );
214 }
215 else
216 {
217 printf("MLI::setCyclesAtLevel ERROR : wrong level = %d\n",level);
218 exit(1);
219 }
220 return 0;
221 }
222
223 /*****************************************************************************
224 * set ML method
225 *---------------------------------------------------------------------------*/
226
setMethod(MLI_Method * object)227 int MLI::setMethod( MLI_Method *object )
228 {
229 #ifdef MLI_DEBUG_DETAILED
230 printf("MLI::setMethod = %s\n", object->getName());
231 #endif
232 methodPtr_ = object;
233 return 0;
234 }
235
236 /*****************************************************************************
237 * set up the grid hierarchy
238 *---------------------------------------------------------------------------*/
239
setup()240 int MLI::setup()
241 {
242 int nlevels, status=0;
243 char paramString[100];
244
245 currIter_ = 0;
246 buildTime_ = MLI_Utils_WTime();
247 sprintf( paramString, "setOutputLevel %d", outputLevel_ );
248 methodPtr_->setParams(paramString, 0, NULL);
249 nlevels = methodPtr_->setup(this);
250 coarsestLevel_ = nlevels - 1;
251 buildTime_ = MLI_Utils_WTime() - buildTime_;
252 for (int i = 0; i < nlevels; i++) status += oneLevels_[i]->setup();
253 if ( coarseSolver_ != NULL )
254 {
255 oneLevels_[coarsestLevel_]->setCoarseSolve(coarseSolver_);
256 coarseSolver_ = NULL;
257 }
258 assembled_ = 1;
259 return status;
260 }
261
262 /*****************************************************************************
263 * perform one cycle
264 *---------------------------------------------------------------------------*/
265
cycle(MLI_Vector * sol,MLI_Vector * rhs)266 int MLI::cycle( MLI_Vector *sol, MLI_Vector *rhs )
267 {
268 oneLevels_[0]->setSolutionVector( sol );
269 oneLevels_[0]->setRHSVector( rhs );
270 int status = oneLevels_[0]->solve1Cycle();
271 return status;
272 }
273
274 /*****************************************************************************
275 * perform solve
276 *---------------------------------------------------------------------------*/
277
solve(MLI_Vector * sol,MLI_Vector * rhs)278 int MLI::solve( MLI_Vector *sol, MLI_Vector *rhs )
279 {
280 int iter=0, mypid;
281 double norm2, relTol, oldNorm2, zero=0.0;
282 MLI_Matrix *Amat;
283 MLI_Vector *res;
284 #if 0
285 char paramString[30];
286 MLI_Solver *preSmoother;
287 #endif
288
289 /*-------------------------------------------------------------------*/
290 /* check for error */
291 /*-------------------------------------------------------------------*/
292
293 if ( ! assembled_ )
294 {
295 printf("MLI::solve ERROR - setup not called yet.\n");
296 exit(1);
297 }
298
299 /*-------------------------------------------------------------------*/
300 /* if coarse solver was set before setup, put it in the coarse level */
301 /*-------------------------------------------------------------------*/
302
303 if ( coarseSolver_ != NULL )
304 {
305 oneLevels_[coarsestLevel_]->setCoarseSolve(coarseSolver_);
306 coarseSolver_ = NULL;
307 }
308
309 /*-------------------------------------------------------------------*/
310 /* compute initial residual norm and convergence tolerance */
311 /*-------------------------------------------------------------------*/
312
313 MPI_Comm_rank(mpiComm_, &mypid);
314 res = oneLevels_[0]->getResidualVector();
315 Amat = oneLevels_[0]->getAmat();
316 solveTime_ = MLI_Utils_WTime();
317 if ( maxIterations_ == 1 )
318 {
319 norm2 = 1.0;
320 relTol = 0.1;
321 sol->setConstantValue(zero);
322 #if 0
323 strcpy( paramString, "zeroInitialGuess" );
324 preSmoother = oneLevels_[0]->getPreSmoother();
325 if (preSmoother != NULL) preSmoother->setParams(paramString, 0, NULL);
326 #endif
327 }
328 else
329 {
330 Amat->apply( -1.0, sol, 1.0, rhs, res );
331 norm2 = res->norm2();
332 relTol = tolerance_ * norm2;
333 if ( outputLevel_ > 0 && currIter_ == 0 )
334 printf("\tMLI Initial norm = %16.8e (%16.8e)\n", norm2, relTol);
335 }
336
337 while ( norm2 > relTol && iter < maxIterations_ )
338 {
339 iter++;
340 currIter_++;
341 cycle( sol, rhs );
342 if ( maxIterations_ > 1 )
343 {
344 Amat->apply( -1.0, sol, 1.0, rhs, res );
345 oldNorm2 = norm2;
346 norm2 = res->norm2();
347 if ( outputLevel_ > 0 && mypid == 0 && maxIterations_ > 1 )
348 printf("\tMLI iteration = %5d, rnorm = %14.6e (%14.6e)\n",
349 currIter_, norm2, norm2/oldNorm2);
350 }
351 if ( iter < maxIterations_ )
352 {
353 oneLevels_[0]->resetSolutionVector();
354 oneLevels_[0]->resetRHSVector();
355 }
356 }
357 solveTime_ = MLI_Utils_WTime() - solveTime_;
358 if ( norm2 > tolerance_ || iter >= maxIterations_ ) return 1;
359 else return 0;
360 }
361
362 /*****************************************************************************
363 * print
364 *---------------------------------------------------------------------------*/
365
print()366 int MLI::print()
367 {
368 int mypid;
369 MPI_Comm_rank(mpiComm_, &mypid);
370 if ( mypid == 0 )
371 {
372 printf("\t***************** MLI Information *********************\n");
373 printf("\t*** maxLevels = %d\n", maxLevels_);
374 printf("\t*** output level = %d\n", outputLevel_);
375 printf("\t*** max iterations = %d\n", maxIterations_);
376 printf("\t*** tolerance = %e\n", tolerance_);
377 printf("\t*******************************************************\n");
378 }
379 return 0;
380 }
381
382 /*****************************************************************************
383 * output timing information
384 *---------------------------------------------------------------------------*/
385
printTiming()386 int MLI::printTiming()
387 {
388 int mypid;
389
390 MPI_Comm_rank( mpiComm_, &mypid );
391 if ( mypid == 0 )
392 {
393 printf("\t***************** MLI Timing Information **************\n");
394 printf("\t*** MLI Build time = %e seconds\n", buildTime_);
395 printf("\t*** MLI Solve time = %e seconds\n", solveTime_);
396 printf("\t*******************************************************\n");
397 }
398 return 0;
399 }
400
401 /*****************************************************************************
402 * get oneLevel object
403 *---------------------------------------------------------------------------*/
404
getOneLevelObject(int level)405 MLI_OneLevel* MLI::getOneLevelObject( int level )
406 {
407 if ( level >= 0 && level < maxLevels_ ) return oneLevels_[level];
408 else
409 {
410 printf("MLI::getOneLevelObject ERROR : wrong level = %d\n", level);
411 return NULL;
412 }
413 }
414
415 /*****************************************************************************
416 * get system matrix
417 *---------------------------------------------------------------------------*/
418
getSystemMatrix(int level)419 MLI_Matrix* MLI::getSystemMatrix( int level )
420 {
421 if ( level >= 0 && level < maxLevels_ )
422 return oneLevels_[level]->getAmat();
423 else
424 {
425 printf("MLI::getSystemMatrix ERROR : wrong level = %d\n", level);
426 return NULL;
427 }
428 }
429
430 /*****************************************************************************
431 * get prolongation operator
432 *---------------------------------------------------------------------------*/
433
getProlongation(int level)434 MLI_Matrix* MLI::getProlongation( int level )
435 {
436 if ( level >= 0 && level < maxLevels_ )
437 return oneLevels_[level]->getPmat();
438 else
439 {
440 printf("MLI::getProlongation ERROR : wrong level = %d\n", level);
441 return NULL;
442 }
443 }
444
445 /*****************************************************************************
446 * get restriction operator
447 *---------------------------------------------------------------------------*/
448
getRestriction(int level)449 MLI_Matrix* MLI::getRestriction( int level )
450 {
451 if ( level >= 0 && level < maxLevels_ )
452 return oneLevels_[level]->getRmat();
453 else
454 {
455 printf("MLI::getRestriction ERROR : wrong level = %d\n", level);
456 return NULL;
457 }
458 }
459
460 /*****************************************************************************
461 * get smoother
462 *---------------------------------------------------------------------------*/
463
getSmoother(int level,int pre_post)464 MLI_Solver* MLI::getSmoother( int level, int pre_post )
465 {
466 if ( level >= 0 && level < maxLevels_ )
467 {
468 if ( pre_post == MLI_SMOOTHER_PRE )
469 return oneLevels_[level]->getPreSmoother();
470 else if ( pre_post == MLI_SMOOTHER_POST )
471 return oneLevels_[level]->getPostSmoother();
472 else
473 {
474 printf("MLI::getSmoother ERROR : pre or post ? \n");
475 return ((MLI_Solver *) NULL);
476 }
477 }
478 else
479 {
480 printf("MLI::getRestriction ERROR : wrong level = %d\n", level);
481 return ((MLI_Solver *) NULL);
482 }
483 }
484
485 /*****************************************************************************
486 * get fedata
487 *---------------------------------------------------------------------------*/
488
getFEData(int level)489 MLI_FEData* MLI::getFEData( int level )
490 {
491 if ( level >= 0 && level < maxLevels_ )
492 return oneLevels_[level]->getFEData();
493 else
494 {
495 printf("MLI::getFEData ERROR : wrong level = %d\n", level);
496 return ((MLI_FEData *) NULL);
497 }
498 }
499
500 /*****************************************************************************
501 * get sfei
502 *---------------------------------------------------------------------------*/
503
getSFEI(int level)504 MLI_SFEI* MLI::getSFEI( int level )
505 {
506 if ( level >= 0 && level < maxLevels_ )
507 return oneLevels_[level]->getSFEI();
508 else
509 {
510 printf("MLI::getSFEI ERROR : wrong level = %d\n", level);
511 return ((MLI_SFEI *) NULL);
512 }
513 }
514
515 /*****************************************************************************
516 * get node to equation map
517 *---------------------------------------------------------------------------*/
518
getNodeEqnMap(int level)519 MLI_Mapper* MLI::getNodeEqnMap( int level )
520 {
521 if ( level >= 0 && level < maxLevels_ )
522 return oneLevels_[level]->getNodeEqnMap();
523 else
524 {
525 printf("MLI::getNodeEqnMap ERROR : wrong level = %d\n", level);
526 return ((MLI_Mapper *) NULL);
527 }
528 }
529
530 /*****************************************************************************
531 * reset discretization matrix
532 *---------------------------------------------------------------------------*/
533
resetSystemMatrix(int level)534 int MLI::resetSystemMatrix( int level )
535 {
536 #ifdef MLI_DEBUG_DETAILED
537 printf("MLI::resetSystemMatrix, level = %d\n", level);
538 #endif
539 if ( level >= 0 && level < maxLevels_ ) oneLevels_[level]->resetAmat();
540 else
541 {
542 printf("MLI::resetSystemMatrix ERROR : wrong level = %d\n", level);
543 exit(1);
544 }
545 return 0;
546 }
547
548