ResidualEstimator.cc 14.5 KB
Newer Older
1
2
3
4
5
6
7
8
#include "ResidualEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "Assembler.h"
#include "Traverse.h"
#include "Parameters.h"

9
10
11
12
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "mpi.h"
#endif

13
14
namespace AMDiS {

15
  ResidualEstimator::ResidualEstimator(std::string name, int r) 
16
17
18
    : Estimator(name, r),
      C0(1.0), 
      C1(1.0), 
19
      C2(0.0), 
20
21
      C3(1.0)
  {
22
23
    FUNCNAME("ResidualEstimator::ResidualEstimator()");

24
25
26
27
28
29
30
31
32
    GET_PARAMETER(0, name + "->C0", "%f", &C0);
    GET_PARAMETER(0, name + "->C1", "%f", &C1);
    GET_PARAMETER(0, name + "->C2", "%f", &C2);
    GET_PARAMETER(0, name + "->C3", "%f", &C3);

    C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
    C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
    C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
    C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
33
34

    TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
35
36
  }

37

38
39
40
  void ResidualEstimator::init(double ts)
  {
    FUNCNAME("ResidualEstimator::init()");
41

42
43
44
45
    timestep = ts;

    mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
46
    nSystems = static_cast<int>(uh.size());
47

Thomas Witkowski's avatar
Thomas Witkowski committed
48
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
49
50

    dim = mesh->getDim();
51
52
    basFcts = new const BasisFunction*[nSystems];
    quadFast = new FastQuadrature*[nSystems];
53
54

    degree = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
55
    for (int system = 0; system < nSystems; system++) {
56
      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
57
      degree = std::max(degree, basFcts[system]->getDegree());
58
59
60
61
    }
    degree *= 2;

    quad = Quadrature::provideQuadrature(dim, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    nPoints = quad->getNumPoints();
63
64

    Flag flag = INIT_PHI | INIT_GRD_PHI;
65
66
    if (degree > 2)
      flag |= INIT_D2_PHI;    
67

68
    for (int system = 0; system < nSystems; system++)
69
70
      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
							       *quad, 
71
							       flag);    
72
  
73
74
75
    uhEl = new double*[nSystems];
    uhNeigh = new double*[nSystems];
    uhOldEl = timestep ? new double*[nSystems] : NULL;
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    for (int system = 0; system < nSystems; system++) {
78
79
      uhEl[system] = new double[basFcts[system]->getNumber()]; 
      uhNeigh[system] = new double[basFcts[system]->getNumber()];
80
      if (timestep)
81
	uhOldEl[system] = new double[basFcts[system]->getNumber()];
82
83
    }

84
85
86
    uhQP = timestep ? new double[nPoints] : NULL;
    uhOldQP = timestep ? new double[nPoints] : NULL;
    riq = new double[nPoints];
87
88
89
    grdUh_qp = NULL;
    D2uhqp = NULL;

90
    // clear error indicators and mark elements for jumpRes
91
92
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
    while (elInfo) {
      elInfo->getElement()->setEstimation(0.0, row);
      elInfo->getElement()->setMark(1);
      elInfo = stack.traverseNext(elInfo);
    }

    est_sum = 0.0;
    est_max = 0.0;
    est_t_sum = 0.0;
    est_t_max = 0.0;

    traverseFlag = 
      Mesh::FILL_NEIGH      |
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND      |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_DET        |
      Mesh::CALL_LEAF_EL;
112
113
114
    neighInfo = mesh->createNewElInfo();

    // prepare date for computing jump residual
115
116
117
118
119
120
121
122
123
124
    if (C1 > 0.0 && dim > 1) {
      surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
      nPointsSurface = surfaceQuad->getNumPoints();
      grdUhEl.resize(nPointsSurface);
      grdUhNeigh.resize(nPointsSurface);
      jump.resize(nPointsSurface);
      localJump.resize(nPointsSurface);
      nNeighbours = Global::getGeo(NEIGH, dim);
      lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
      lambda = new DimVec<double>(dim, NO_INIT);
125
    }
126
127
  }

128

129
130
131
132
  void ResidualEstimator::exit(bool output)
  {
    FUNCNAME("ResidualEstimator::exit()");

133
134
135
136
137
138
139
140
141
142
143
144
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double send_est_sum = est_sum;
    double send_est_max = est_max;
    double send_est_t_sum = est_t_sum;
    double send_est_t_max = est_t_max;

    MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX);
    MPI::COMM_WORLD.Allreduce(&send_est_t_sum, &est_t_sum, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&send_est_t_max, &est_t_max, 1, MPI_DOUBLE, MPI_MAX);
#endif

145
146
147
    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

Thomas Witkowski's avatar
Thomas Witkowski committed
148
    for (int system = 0; system < nSystems; system++) {
149
150
      delete [] uhEl[system];
      delete [] uhNeigh[system];
151
      if (timestep)
152
	delete [] uhOldEl[system];
153
154
    }

155
156
    delete [] uhEl;
    delete [] uhNeigh;
157
158

    if (timestep) {
159
160
161
      delete [] uhOldEl;
      delete [] uhQP;
      delete [] uhOldQP;
162
    } else {
163
164
      if (uhQP != NULL)
	delete [] uhQP;
165
166
167
168
    }

    if (output) {
      MSG("estimate   = %.8e\n", est_sum);
169
      if (C3)
170
171
172
	MSG("time estimate   = %.8e\n", est_t_sum);
    }

173
174
175
    delete [] riq;
    delete [] basFcts;
    delete [] quadFast;
176

177
178
179
180
    if (grdUh_qp != NULL)
      delete [] grdUh_qp;
    if (D2uhqp != NULL)
      delete [] D2uhqp;
181

182
    if (C1 && (dim > 1)) {
183
184
      delete lambdaNeigh;
      delete lambda;
185
186
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
187
    delete neighInfo;
188
189
  }

190

191
  void ResidualEstimator::estimateElement(ElInfo *elInfo)
192
  {    
193
194
    FUNCNAME("ResidualEstimator::estimateElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
195
    TEST_EXIT_DBG(nSystems > 0)("no system set\n");
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    double val = 0.0;
198
    std::vector<Operator*>::iterator it;
199
    std::vector<double*>::iterator itfac;
Thomas Witkowski's avatar
Thomas Witkowski committed
200
    Element *el = elInfo->getElement();
201
    double det = elInfo->getDet();
202
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    double est_el = el->getEstimation(row);
204
205
    double h2 = h2_from_det(det, dim);

206
    for (int iq = 0; iq < nPoints; iq++)
207
208
      riq[iq] = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
209
    for (int system = 0; system < nSystems; system++) {
210
211
212
      if (matrix[system] == NULL) 
	continue;

213
214
215
      DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
      DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);

216
217
      // === init assemblers ===

218
219
      for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
	   it != dofMat->getOperatorsEnd(); ++it, ++itfac)
220
	if (*itfac == NULL || **itfac != 0.0)
221
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
222

223
      if (C0 > 0.0)
224
	for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
225
	  (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);	
Naumann, Andreas's avatar
Naumann, Andreas committed
226
	
227
      if (timestep && uhOld[system]) {
228
229
230
	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
	uhOld[system]->getLocalVector(el, uhOldEl[system]);
  
231
232
233
	// === Compute time error. ===

	if (C0 > 0.0 || C3 > 0.0) {   
234
235
236
	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
	  
237
	  if (C3 > 0.0 && uhOldQP && system == std::max(row, 0)) {
238
	    val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
239
	    for (int iq = 0; iq < nPoints; iq++) {
240
241
242
243
244
245
246
247
248
	      double tiq = (uhQP[iq] - uhOldQP[iq]);
	      val += quad->getWeight(iq) * tiq * tiq;
	    }
	    double v = C3 * det * val;
	    est_t_sum += v;
	    est_t_max = max(est_t_max, v);
	  }
	}
      }
249
           
250
251
      // === Compute element residual. ===
      if (C0 > 0.0) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
252
	for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
253
	     it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
254
	  if (*itfac == NULL || **itfac != 0.0) {
255
	    if (uhQP == NULL && (*it)->zeroOrderTerms()) {
256
	      uhQP = new double[nPoints];
257
258
	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
	    }
259
260
	    if (grdUh_qp == NULL && 
		((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
261
262
263
	      grdUh_qp = new WorldVector<double>[nPoints];
	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
	    }
264
	    if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) { 
265
266
267
	      D2uhqp = new WorldMatrix<double>[nPoints];
	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
	    }
268
269
270
	  }
	}
	
271
272
	// === Compute the element residual and store it in irq. ===

273
	r(elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
274
	  nPoints, 
275
276
277
278
279
280
	  uhQP,
	  grdUh_qp,
	  D2uhqp,
	  uhOldQP,
	  NULL,  // grdUhOldQP 
	  NULL,  // D2UhOldQP
Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
	  dofMat, 
	  dofVec,
283
	  quad,
284
	  riq);
285
286
287
288
      }     
    }

    // add integral over r square
289
    val = 0.0;
290
    for (int iq = 0; iq < nPoints; iq++)
291
      val += quad->getWeight(iq) * riq[iq] * riq[iq];
292
   
293
    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
294
      val = C0 * h2 * h2 * det * val;
295
    else
296
      val = C0 * h2 * det * val;
297
    
298
299
    est_el += val;

300
301
302

    // === Compute jump residuals. ===

303
    if (C1 && dim > 1) {
304
      int dow = Global::getGeo(WORLD);
305

306
      for (int face = 0; face < nNeighbours; face++) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
307
	Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
308
	if (neigh && neigh->getMark()) {     
309
310
	  int oppV = elInfo->getOppVertex(face);
	      
311
312
	  el->sortFaceIndices(face, &faceIndEl);
	  neigh->sortFaceIndices(oppV, &faceIndNeigh);
313
314
315
	    
	  neighInfo->setElement(const_cast<Element*>(neigh));
	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
316
	      	
317
318
319
320
321
322
323
324
325
	  for (int i = 0; i < dow; i++)
	    neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
		
	  // periodic leaf data ?
	  ElementData *ldp = el->getElementData()->getElementData(PERIODIC);

	  bool periodicCoords = false;

	  if (ldp) {
326
327
	    std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	    std::list<LeafDataPeriodic::PeriodicInfo>& infoList = 
328
		dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList();
329

330
	    for (it = infoList.begin(); it != infoList.end(); ++it) {
331
332
	      if (it->elementSide == face) {
		for (int i = 0; i < dim; i++) {
333
334
		  int i1 = faceIndEl[i];
		  int i2 = faceIndNeigh[i];
335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
		  int j = 0;
		  for (; j < dim; j++) {
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
		    if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, 
								   dim),
						      face,
						      j)) {
		      break;
		    }
		  }

		  TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
		      
		  neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
		}
		periodicCoords = true;
		break;
	      }
	    }
	  }
355

356
357
	  if (!periodicCoords) {
	    for (int i = 0; i < dim; i++) {
358
359
	      int i1 = faceIndEl[i];
	      int i2 = faceIndNeigh[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
360
	      for (int j = 0; j < dow; j++)
361
362
363
		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
	    }
	  }
364

365
	  Parametric *parametric = mesh->getParametric();
366
367
	  if (parametric)
	    neighInfo = parametric->addParametricInfo(neighInfo);	  
368
	      
369
	  double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
370
	      
371
372
	  for (int iq = 0; iq < nPointsSurface; iq++)
	    jump[iq].set(0.0);     
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
	  for (int system = 0; system < nSystems; system++) {	
375
376
377
	    if (matrix[system] == NULL) 
	      continue;
	      
378
379
	    uh[system]->getLocalVector(el, uhEl[system]);	
	    uh[system]->getLocalVector(neigh, uhNeigh[system]);
380
			
381
382
383
384
	    for (int iq = 0; iq < nPointsSurface; iq++) {
	      (*lambda)[face] = 0.0;
	      for (int i = 0; i < dim; i++)
		(*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
385
		  
386
	      basFcts[system]->evalGrdUh(*lambda, 
387
					 grdLambda, 
388
					 uhEl[system], 
389
					 &grdUhEl[iq]);
390
		  
391
392
393
	      (*lambda)[oppV] = 0.0;
	      for (int i = 0; i < dim; i++)
		(*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
394
		  
395
396
	      basFcts[system]->evalGrdUh(*lambda, 
					 *lambdaNeigh, 
397
					 uhNeigh[system], 
398
					 &grdUhNeigh[iq]);
399
		  
400
	      grdUhEl[iq] -= grdUhNeigh[iq];
401
402
	    }				

403
	    std::vector<double*>::iterator fac;
404
405
406
407
408

	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
		 ++it, ++fac) {
409

410
	      if (*fac == NULL || **fac != 0.0) {
411
412
		for (int iq = 0; iq < nPointsSurface; iq++)
		  localJump[iq].set(0.0);
413
		
414
415
416
		(*it)->weakEvalSecondOrder(nPointsSurface,
					   grdUhEl.getValArray(),
					   localJump.getValArray());
417
		double factor = *fac ? **fac : 1.0;
418
419
420
		if (factor != 1.0)
		  for (int i = 0; i < nPointsSurface; i++)
		    localJump[i] *= factor;
421
		
422
423
		for (int i = 0; i < nPointsSurface; i++)
		  jump[i] += localJump[i];
424
425
	      }		
	    }
426
	  }
427

428
	  val = 0.0;
429
430
	  for (int iq = 0; iq < nPointsSurface; iq++)
	    val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
431
432
433
434
435
436
437
438
	      
	  double d = 0.5 * (det + detNeigh);

	  if (norm == NO_NORM || norm == L2_NORM)
	    val *= C1 * h2_from_det(d, dim) * d;
	  else
	    val *= C1 * d;
	      
439
	  if (parametric)
440
441
442
443
444
445
446
	    neighInfo = parametric->removeParametricInfo(neighInfo);

	  neigh->setEstimation(neigh->getEstimation(row) + val, row);
	  est_el += val;
	} 
      } 
       
Thomas Witkowski's avatar
Thomas Witkowski committed
447
      val = fh[std::max(row, 0)]->getBoundaryManager()->
448
	boundResidual(elInfo, matrix[std::max(row, 0)], uh[std::max(row, 0)]);
Thomas Witkowski's avatar
Thomas Witkowski committed
449

450
451
452
453
454
455
456
457
458
459
460
461
462
463
      if (norm == NO_NORM || norm == L2_NORM)
	val *= C1 * h2;
      else
	val *= C1;
	
      est_el += val;
    } 
  

    el->setEstimation(est_el, row);

    est_sum += est_el;
    est_max = max(est_max, est_el);

464
    elInfo->getElement()->setMark(0);  
465
466
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
467

468
  void r(const ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
469
	 int nPoints,
470
	 const double *uhIq,
471
472
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
473
	 const double *uhOldIq,
474
475
476
477
478
479
480
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
	 double *result)
  {
481
482
    std::vector<Operator*>::iterator it;
    std::vector<double*>::iterator fac;
483
484

    // lhs
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
	 it != A->getOperatorsEnd(); ++it, ++fac) {
487
488
489
     
      double factor = *fac ? **fac : 1.0;

490
      if (factor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
	if (D2UhIq)
	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	
493
494

	if (grdUhIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
495
496
	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
497
498
	}
	
Thomas Witkowski's avatar
Thomas Witkowski committed
499
500
	if (uhIq)
	  (*it)->evalZeroOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);	
501
502
503
504
      }
    }
    
    // rhs
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
	 it != fh->getOperatorsEnd(); ++it, ++fac) {
507

508
509
      double factor = *fac ? **fac : 1.0;

510
511
      if (factor) {
	if ((*it)->getUhOld()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
	  if (D2UhOldIq)
	    (*it)->evalSecondOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, factor);
	  
515
	  if (grdUhOldIq) {
Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
	    (*it)->evalFirstOrderGrdPsi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
	    (*it)->evalFirstOrderGrdPhi(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
518
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
519
520
521

	  if (uhOldIq)
	    (*it)->evalZeroOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);	  
522
	} else {
523
	  std::vector<double> fx(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
524
	  (*it)->getC(elInfo, nPoints, fx);
525

526
	  for (int iq = 0; iq < nPoints; iq++)
527
528
529
530
531
532
533
534
	    result[iq] -= factor * fx[iq];
	}
      }
    }    
  }


}