ResidualParallelEstimator.cc 2.8 KB
Newer Older
1
2
3
4
5
6
#include "ResidualParallelEstimator.h"
#include "ResidualEstimator.h"
#include "TraverseParallel.h"
#include "Parameters.h"
#include "OpenMP.h"

7
8
#ifdef _OPENMP

9
namespace AMDiS {
10
  ResidualParallelEstimator::ResidualParallelEstimator(std::string name, int r)
11
12
13
14
15
16
    : Estimator(name, r),
      C3(1.0)
  {
    GET_PARAMETER(0, name + "->C3", "%f", &C3);
    C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;

17
    seqEstimators_.resize(omp_get_overall_max_threads());
18

Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
19
    for (int i = 0; i < omp_get_overall_max_threads(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
20
      seqEstimators_[i] = new ResidualEstimator(name, r);
21
22
23
24
  }

  ResidualParallelEstimator::~ResidualParallelEstimator()
  {
Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
25
26
    for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++)
      delete seqEstimators_[i];    
27
28
29
30
31
32
33
34
35
  }

  void ResidualParallelEstimator::addSystem(DOFMatrix *matrix_,
					    DOFVector<double> *uh_,
					    DOFVector<double> *fh_,
					    DOFVector<double> *uhOld_)
  {
    Estimator::addSystem(matrix_, uh_, fh_, uhOld_);

Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
36
    for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++)
37
38
39
40
41
42
43
44
      seqEstimators_[i]->addSystem(matrix_, uh_, fh_, uhOld_);
  }

  void ResidualParallelEstimator::addUhOldToSystem(int system, 
						   DOFVector<double> *uhOld_)
  {
    Estimator::addUhOldToSystem(system, uhOld_);

Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
45
    for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++)
46
47
48
49
50
      seqEstimators_[i]->addUhOldToSystem(system, uhOld_);
  }

  double ResidualParallelEstimator::estimate(double ts)
  {
51
52
    FUNCNAME("ResidualParallelEstimator::estimate()");

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

Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
55
    for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++)
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
      seqEstimators_[i]->init(ts);

    TraverseParallelStack stack;

#pragma omp parallel
    {
      ResidualEstimator *estimator = seqEstimators_[omp_get_thread_num()];     
   
      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;
   
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);

      while (elInfo) {
	estimator->estimateElement(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }
    }

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

    for (int i = 0; i < static_cast<int>(seqEstimators_.size()); i++) {
      seqEstimators_[i]->exit(false);
      est_sum += pow(seqEstimators_[i]->getErrorSum(), 2);
      est_max = max(est_max, seqEstimators_[i]->getErrorMax());
      est_t_sum += pow(seqEstimators_[i]->getTimeEst(), 2);
      est_t_max = max(est_t_max, seqEstimators_[i]->getTimeEstMax());
    }

    est_sum = sqrt(est_sum);
    est_t_sum = sqrt(est_t_sum);

    MSG("estimate   = %.8e\n", est_sum);
Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
98
    if (C3)
99
100
101
102
103
104
105
106
      MSG("time estimate   = %.8e\n", est_t_sum);

    return est_sum;
  }


}

107
#endif // _OPENMP