LORENE
cmp_test_poisson.C
1 /*
2  * Method of class Cmp 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) 2000-2001 Eric Gourgoulhon
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 char cmp_test_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_test_poisson.C,v 1.3 2014/10/13 08:52:48 j_novak Exp $" ;
32 
33 /*
34  * $Id: cmp_test_poisson.C,v 1.3 2014/10/13 08:52:48 j_novak Exp $
35  * $Log: cmp_test_poisson.C,v $
36  * Revision 1.3 2014/10/13 08:52:48 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.2 2002/10/16 14:36:34 j_novak
40  * Reorganization of #include instructions of standard C++, in order to
41  * use experimental version 3 of gcc.
42  *
43  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
44  * LORENE
45  *
46  * Revision 2.0 2000/10/05 14:22:24 eric
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_test_poisson.C,v 1.3 2014/10/13 08:52:48 j_novak Exp $
51  *
52  */
53 
54 // Headers Lorene
55 #include "cmp.h"
56 
57 namespace Lorene {
58 Tbl Cmp::test_poisson(const Cmp& uu, ostream& ostr, bool detail) const {
59 
60  assert( uu.get_mp() == mp ) ;
61 
62  // Computation of the absolute and relative errors
63  // -----------------------------------------------
64 
65  int dzi ;
66  if ( check_dzpuis(4) ) {
67  dzi = 4 ;
68  }
69  else{
70  if ( check_dzpuis(2) ) {
71  dzi = 2 ;
72  }
73  else{
74  assert( check_dzpuis(0) ) ;
75  dzi = 0 ;
76  }
77  }
78 
79  Tbl tdiff = max( abs( uu.laplacien(dzi) - *this ) ) ;
80 
81  Tbl tmax = max( abs(*this) ) ;
82 
83  int nz = mp->get_mg()->get_nzone() ;
84  int nzm1 = nz - 1 ;
85 
86  Tbl trel(nz) ;
87  trel.set_etat_qcq() ;
88 
89  if ( (dzpuis == 0) || (tmax(nzm1) == double(0)) ) {
90 
91  double s_max = max( tmax ) ;
92 
93  for (int l=0; l<nz; l++) {
94  trel.set(l) = tdiff(l) / s_max ;
95  }
96 
97  }
98  else{
99 
100  double s_max = 0 ;
101  for (int l=0; l<nzm1; l++) {
102  s_max = (tmax(l) > s_max) ? tmax(l) : s_max ;
103  }
104 
105  for (int l=0; l<nzm1; l++) {
106  trel.set(l) = tdiff(l) / s_max ;
107  }
108 
109  trel.set(nzm1) = tdiff(nzm1) / tmax(nzm1) ;
110 
111  }
112 
113  // All the errors are set in the same output 2-D Tbl
114  // -------------------------------------------------
115 
116  Tbl err(3, nz) ;
117 
118  err.set_etat_qcq() ;
119 
120  for(int l=0; l<nz; l++) {
121  err.set(0, l) = trel(l) ;
122  err.set(1, l) = tdiff(l) ;
123  err.set(2, l) = tmax(l) ;
124  }
125 
126  // Display
127  // -------
128 
129  if (detail) {
130  ostr << "Max. source :" ;
131  for (int l=0; l<nz; l++) {
132  ostr << " " << err(2, l) ;
133  }
134 
135  ostr << endl << "Abs. error : " ;
136  for (int l=0; l<nz; l++) {
137  ostr << " " << err(1, l) ;
138  }
139  }
140 
141  ostr << endl << "Rel. error : " ;
142  for (int l=0; l<nz; l++) {
143  ostr << " " << err(0, l) ;
144  }
145 
146  ostr << endl ;
147 
148  return err ;
149 
150 }
151 
152 }
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Cmp::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: cmp.C:715
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::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Lorene::Cmp::test_poisson
Tbl test_poisson(const Cmp &uu, ostream &ostr, bool detail=false) const
Checks if a Poisson equation with *this as a source has been correctly solved.
Definition: cmp_test_poisson.C:58
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Cmp::dzpuis
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified z...
Definition: cmp.h:461
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Cmp::mp
const Map * mp
Reference mapping.
Definition: cmp.h:451
Lorene::Cmp::laplacien
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
Definition: cmp_deriv.C:242
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901