LORENE
scalar_test_poisson.C
1 /*
2  * Method of class Scalar to check if a Poisson equation has been correctly
3  * solved.
4  *
5  * (see file cmp.h for documentation)
6  *
7  */
8 
9 /*
10  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
11  *
12  * Copyright (c) 2000-2001 Eric Gourgoulhon (for a preceding Cmp version)
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation; either version 2 of the License, or
19  * (at your option) any later version.
20  *
21  * LORENE is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with LORENE; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  *
30  */
31 
32 
33 
34 
35 /*
36  * $Id: scalar_test_poisson.C,v 1.5 2016/12/05 16:18:19 j_novak Exp $
37  * $Log: scalar_test_poisson.C,v $
38  * Revision 1.5 2016/12/05 16:18:19 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.4 2014/10/13 08:53:47 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.3 2003/10/29 13:14:29 e_gourgoulhon
45  * Change of method name: laplacien --> laplacian.
46  *
47  * Revision 1.2 2003/10/01 13:04:44 e_gourgoulhon
48  * The method Tensor::get_mp() returns now a reference (and not
49  * a pointer) onto a mapping.
50  *
51  * Revision 1.1 2003/09/25 09:13:11 e_gourgoulhon
52  * First version.
53  *
54  *
55  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_test_poisson.C,v 1.5 2016/12/05 16:18:19 j_novak Exp $
56  *
57  */
58 
59 // Headers Lorene
60 #include "tensor.h"
61 
62 namespace Lorene {
63 Tbl Scalar::test_poisson(const Scalar& uu, ostream& ostr, bool detail) const {
64 
65  assert( &(uu.get_mp()) == mp ) ;
66 
67  // Computation of the absolute and relative errors
68  // -----------------------------------------------
69 
70  int dzi ;
71  if ( check_dzpuis(4) ) {
72  dzi = 4 ;
73  }
74  else{
75  if ( check_dzpuis(2) ) {
76  dzi = 2 ;
77  }
78  else{
79  assert( check_dzpuis(0) ) ;
80  dzi = 0 ;
81  }
82  }
83 
84  Tbl tdiff = max( abs( uu.laplacian(dzi) - *this ) ) ;
85 
86  Tbl tmax = max( abs(*this) ) ;
87 
88  int nz = mp->get_mg()->get_nzone() ;
89  int nzm1 = nz - 1 ;
90 
91  Tbl trel(nz) ;
92  trel.set_etat_qcq() ;
93 
94  if ( (dzpuis == 0) || (tmax(nzm1) == double(0)) ) {
95 
96  double s_max = max( tmax ) ;
97 
98  for (int l=0; l<nz; l++) {
99  trel.set(l) = tdiff(l) / s_max ;
100  }
101 
102  }
103  else{
104 
105  double s_max = 0 ;
106  for (int l=0; l<nzm1; l++) {
107  s_max = (tmax(l) > s_max) ? tmax(l) : s_max ;
108  }
109 
110  for (int l=0; l<nzm1; l++) {
111  trel.set(l) = tdiff(l) / s_max ;
112  }
113 
114  trel.set(nzm1) = tdiff(nzm1) / tmax(nzm1) ;
115 
116  }
117 
118  // All the errors are set in the same output 2-D Tbl
119  // -------------------------------------------------
120 
121  Tbl err(3, nz) ;
122 
123  err.set_etat_qcq() ;
124 
125  for(int l=0; l<nz; l++) {
126  err.set(0, l) = trel(l) ;
127  err.set(1, l) = tdiff(l) ;
128  err.set(2, l) = tmax(l) ;
129  }
130 
131  // Display
132  // -------
133 
134  if (detail) {
135  ostr << "Max. source :" ;
136  for (int l=0; l<nz; l++) {
137  ostr << " " << err(2, l) ;
138  }
139 
140  ostr << endl << "Abs. error : " ;
141  for (int l=0; l<nz; l++) {
142  ostr << " " << err(1, l) ;
143  }
144  }
145 
146  ostr << endl << "Rel. error : " ;
147  for (int l=0; l<nz; l++) {
148  ostr << " " << err(0, l) ;
149  }
150 
151  ostr << endl ;
152 
153  return err ;
154 
155 }
156 
157 }
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
friend Tbl max(const Scalar &)
Maximum values of a Scalar in each domain.
Definition: scalar_math.C:614
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external do...
Definition: scalar.h:409
Tbl test_poisson(const Scalar &uu, ostream &ostr, bool detail=false) const
Checks if a Poisson equation with *this as a source has been correctly solved.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Basic array class.
Definition: tbl.h:164
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
friend Scalar abs(const Scalar &)
Absolute value.
Definition: scalar_math.C:526
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874