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 char scalar_test_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_test_poisson.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $" ;
34 
35 /*
36  * $Id: scalar_test_poisson.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
37  * $Log: scalar_test_poisson.C,v $
38  * Revision 1.4 2014/10/13 08:53:47 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.3 2003/10/29 13:14:29 e_gourgoulhon
42  * Change of method name: laplacien --> laplacian.
43  *
44  * Revision 1.2 2003/10/01 13:04:44 e_gourgoulhon
45  * The method Tensor::get_mp() returns now a reference (and not
46  * a pointer) onto a mapping.
47  *
48  * Revision 1.1 2003/09/25 09:13:11 e_gourgoulhon
49  * First version.
50  *
51  *
52  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_test_poisson.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
53  *
54  */
55 
56 // Headers Lorene
57 #include "tensor.h"
58 
59 namespace Lorene {
60 Tbl Scalar::test_poisson(const Scalar& uu, ostream& ostr, bool detail) const {
61 
62  assert( &(uu.get_mp()) == mp ) ;
63 
64  // Computation of the absolute and relative errors
65  // -----------------------------------------------
66 
67  int dzi ;
68  if ( check_dzpuis(4) ) {
69  dzi = 4 ;
70  }
71  else{
72  if ( check_dzpuis(2) ) {
73  dzi = 2 ;
74  }
75  else{
76  assert( check_dzpuis(0) ) ;
77  dzi = 0 ;
78  }
79  }
80 
81  Tbl tdiff = max( abs( uu.laplacian(dzi) - *this ) ) ;
82 
83  Tbl tmax = max( abs(*this) ) ;
84 
85  int nz = mp->get_mg()->get_nzone() ;
86  int nzm1 = nz - 1 ;
87 
88  Tbl trel(nz) ;
89  trel.set_etat_qcq() ;
90 
91  if ( (dzpuis == 0) || (tmax(nzm1) == double(0)) ) {
92 
93  double s_max = max( tmax ) ;
94 
95  for (int l=0; l<nz; l++) {
96  trel.set(l) = tdiff(l) / s_max ;
97  }
98 
99  }
100  else{
101 
102  double s_max = 0 ;
103  for (int l=0; l<nzm1; l++) {
104  s_max = (tmax(l) > s_max) ? tmax(l) : s_max ;
105  }
106 
107  for (int l=0; l<nzm1; l++) {
108  trel.set(l) = tdiff(l) / s_max ;
109  }
110 
111  trel.set(nzm1) = tdiff(nzm1) / tmax(nzm1) ;
112 
113  }
114 
115  // All the errors are set in the same output 2-D Tbl
116  // -------------------------------------------------
117 
118  Tbl err(3, nz) ;
119 
120  err.set_etat_qcq() ;
121 
122  for(int l=0; l<nz; l++) {
123  err.set(0, l) = trel(l) ;
124  err.set(1, l) = tdiff(l) ;
125  err.set(2, l) = tmax(l) ;
126  }
127 
128  // Display
129  // -------
130 
131  if (detail) {
132  ostr << "Max. source :" ;
133  for (int l=0; l<nz; l++) {
134  ostr << " " << err(2, l) ;
135  }
136 
137  ostr << endl << "Abs. error : " ;
138  for (int l=0; l<nz; l++) {
139  ostr << " " << err(1, l) ;
140  }
141  }
142 
143  ostr << endl << "Rel. error : " ;
144  for (int l=0; l<nz; l++) {
145  ostr << " " << err(0, l) ;
146  }
147 
148  ostr << endl ;
149 
150  return err ;
151 
152 }
153 
154 }
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Scalar::check_dzpuis
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:873
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Scalar::abs
friend Scalar abs(const Scalar &)
Absolute value.
Definition: scalar_math.C:523
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Scalar::max
friend Tbl max(const Scalar &)
Maximum values of a Scalar in each domain.
Definition: scalar_math.C:611
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Scalar::test_poisson
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.
Definition: scalar_test_poisson.C:60
Lorene::Scalar::dzpuis
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external d...
Definition: scalar.h:403