LORENE
scalar_exp_filter.C
1 /*
2  * Function applying an exponential filter to a Scalar:
3  * sigma(n/N) = exp(alpha*(n/N)^(2p)). See scalar.h for documentation.
4  */
5 
6 /*
7  * Copyright (c) 2007 Jerome Novak
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 
27 
28 /*
29  * $Id: scalar_exp_filter.C,v 1.9 2023/09/01 11:40:14 g_servignat Exp $
30  * $Log: scalar_exp_filter.C,v $
31  * Revision 1.9 2023/09/01 11:40:14 g_servignat
32  * Absolute value of m_q is taken in phi filter
33  *
34  * Revision 1.8 2023/08/31 08:27:26 g_servignat
35  * Added the possibility to filter in the r direction within the ylm filter. An order filtering of 0 means no filtering (for all 3 directions).
36  *
37  * Revision 1.7 2023/08/28 09:53:33 g_servignat
38  * Added ylm filter for Tensor and Scalar in theta + phi directions
39  *
40  * Revision 1.6 2016/12/05 16:18:18 j_novak
41  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42  *
43  * Revision 1.5 2014/10/13 08:53:46 j_novak
44  * Lorene classes and functions now belong to the namespace Lorene.
45  *
46  * Revision 1.4 2014/10/06 15:16:15 j_novak
47  * Modified #include directives to use c++ syntax.
48  *
49  * Revision 1.3 2012/01/17 10:29:27 j_penner
50  * added two routines to handle generalized exponential filtering
51  *
52  * Revision 1.2 2007/10/31 10:50:16 j_novak
53  * Testing whether the coefficients are zero in a given domain.
54  *
55  * Revision 1.1 2007/10/31 10:33:13 j_novak
56  * Added exponential filters to smooth Gibbs-type phenomena.
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.9 2023/09/01 11:40:14 g_servignat Exp $
60  *
61  */
62 
63 // C headers
64 #include <cassert>
65 #include <cmath>
66 
67 // Lorene headers
68 #include "tensor.h"
69 #include "proto.h"
70 
71 namespace Lorene {
72 void Scalar::exponential_filter_r(int lzmin, int lzmax, int p,
73  double alpha) {
74  assert(lzmin >= 0) ;
75  const Mg3d& mgrid = *mp->get_mg() ;
76 #ifndef NDEBUG
77  int nz = mgrid.get_nzone() ;
78 #endif
79  assert(lzmax < nz) ;
80  assert(etat != ETATNONDEF) ;
81  if (etat == ETATZERO) return ;
82  va.coef() ;
83  assert(va.c_cf != 0x0) ;
84  assert(alpha < 0.) ;
85  double alp = log(pow(10., alpha)) ;
86 
87  for (int lz=lzmin; lz<=lzmax; lz++) {
88  if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
89  for (int k=0; k<mgrid.get_np(lz); k++)
90  for (int j=0; j<mgrid.get_nt(lz); j++) {
91  int nr = mgrid.get_nr(lz) ;
92  for (int i=0; i<nr; i++) {
93  double eta = double(i)/double(nr-1) ;
94  va.c_cf->set(lz, k, j, i) *= exp(alp*pow(eta, 2*p)) ;
95  }
96  }
97  }
98  if (va.c != 0x0) {
99  delete va.c ;
100  va.c = 0x0 ;
101  }
102  va.del_deriv() ;
103  del_deriv() ;
104 
105  return ;
106 }
107 
108 void Scalar::sarra_filter_r(int lzmin, int lzmax, double p,
109  double alpha) {
110  assert(lzmin >= 0) ;
111  const Mg3d& mgrid = *mp->get_mg() ;
112 #ifndef NDEBUG
113  int nz = mgrid.get_nzone() ;
114 #endif
115  assert(lzmax < nz) ;
116  assert(etat != ETATNONDEF) ;
117  if (etat == ETATZERO) return ;
118  va.coef() ;
119  assert(va.c_cf != 0x0) ;
120  assert(alpha < 0.) ;
121 
122  for (int lz=lzmin; lz<=lzmax; lz++) {
123  if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
124  for (int k=0; k<mgrid.get_np(lz); k++)
125  for (int j=0; j<mgrid.get_nt(lz); j++) {
126  int nr = mgrid.get_nr(lz) ;
127  for (int i=0; i<nr; i++) {
128  double eta = double(i)/double(nr) ;
129  va.c_cf->set(lz, k, j, i) *= exp(alpha*pow(eta, p)) ;
130  }
131  }
132  }
133  if (va.c != 0x0) {
134  delete va.c ;
135  va.c = 0x0 ;
136  }
137  va.del_deriv() ;
138  del_deriv() ;
139 
140  return ;
141 }
142 void exp_filter_r_all_domains( Scalar& ss, int p, double alpha ) {
143  int nz = ss.get_mp().get_mg()->get_nzone() ;
144  ss.exponential_filter_r(0, nz-1, p, alpha) ;
145  return ;
146 }
147 
148 void Scalar::sarra_filter_r_all_domains( double p, double alpha ) {
149  int nz = get_mp().get_mg()->get_nzone() ;
150  sarra_filter_r(0, nz-1, p, alpha) ;
151  return ;
152 }
153 
154 void Scalar::exponential_filter_ylm(int lzmin, int lzmax, int p,
155  double alpha) {
156  assert(lzmin >= 0) ;
157  const Mg3d& mgrid = *mp->get_mg() ;
158  #ifndef NDEBUG
159  int nz = mgrid.get_nzone() ;
160  #endif
161  assert(lzmax < nz) ;
162  assert(etat != ETATNONDEF) ;
163  if (etat == ETATZERO) return ;
164  double alp = log(pow(10., alpha)) ;
165  va.ylm() ;
166  assert(va.c_cf != 0x0) ;
167  const Base_val& base = va.base ;
168  int l_q, m_q, base_r ;
169 
170  for (int lz=lzmin; lz<=lzmax; lz++)
171  if ((*va.c_cf)(lz).get_etat() == ETATQCQ) {
172  int np = mgrid.get_np(lz) ;
173  int nt = mgrid.get_nt(lz) ;
174  int nr = mgrid.get_nr(lz) ;
175  int lmax = base.give_lmax(mgrid, lz) ;
176  for (int k=0; k<np; k++)
177  for (int j=0; j<nt; j++) {
178  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
179  if (nullite_plm(j, nt, k, np, base) == 1 ) {
180  double eta = double(l_q) / double(lmax) ;
181  double sigma = exp(alp*pow(eta, 2*p)) ;
182  for (int i=0; i<nr; i++)
183  va.c_cf->set(lz, k, j, i) *= sigma ;
184  }
185  }
186  }
187 
188  va.ylm_i() ;
189  if (va.c != 0x0) {
190  delete va.c ;
191  va.c = 0x0 ;
192  }
193  va.del_deriv() ;
194  del_deriv() ;
195  return ;
196  }
197 
198 void Scalar::exponential_filter_ylm_phi(int lzmin, int lzmax, int p_r, int p_tet, int p_phi,
199  double alpha) {
200  assert(lzmin >= 0) ;
201  const Mg3d& mgrid = *mp->get_mg() ;
202 #ifndef NDEBUG
203  int nz = mgrid.get_nzone() ;
204 #endif
205  assert(lzmax < nz) ;
206  assert(etat != ETATNONDEF) ;
207  if (etat == ETATZERO) return ;
208  double alp = log(pow(10., alpha)) ;
209  va.ylm() ;
210  assert(va.c_cf != 0x0) ;
211  const Base_val& base = va.base ;
212  int l_q, m_q, base_r ;
213 
214  for (int lz=lzmin; lz<=lzmax; lz++)
215  if ((*va.c_cf)(lz).get_etat() == ETATQCQ) {
216  int np = mgrid.get_np(lz) ;
217  int nt = mgrid.get_nt(lz) ;
218  int nr = mgrid.get_nr(lz) ;
219  int lmax = base.give_lmax(mgrid, lz) ;
220  for (int k=0; k<np+1; k++)
221  for (int j=0; j<nt; j++) {
222  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
223  if (nullite_plm(j, nt, k, np, base) == 1 ) {
224  double eta_theta = double(l_q) / double(lmax) ;
225  double sigma_theta = (p_tet != 0) ? exp(alp*pow(eta_theta, 2*p_tet)) : 1. ;
226 
227  double eta_phi = (np>1) ? double(abs(m_q)) / double(np) : 0. ;
228  double sigma_phi = (p_phi != 0) ? exp(alp*pow(eta_phi, 2*p_phi)) : 1. ;
229  for (int i=0; i<nr; i++)
230  {
231  double eta_r = double(i) / double(nr-1) ;
232  double sigma_r = (p_r != 0) ? exp(alp*pow(eta_r, 2*p_r)) : 1. ;
233  va.c_cf->set(lz, k, j, i) *= sigma_r * sigma_theta * sigma_phi ;
234  }
235  }
236  }
237  }
238  va.ylm_i() ;
239  if (va.c != 0x0) {
240  delete va.c ;
241  va.c = 0x0 ;
242  }
243  va.del_deriv() ;
244  del_deriv() ;
245  return ;
246 }
247 
248 void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha ) {
249  int nz = ss.get_mp().get_mg()->get_nzone() ;
250  ss.exponential_filter_ylm(0, nz-1, p, alpha) ;
251  return ;
252 }
253 
254 void exp_filter_ylm_all_domains_phi(Scalar& ss, int p_r, int p_tet, int p_phi, double alpha ) {
255  int nz = ss.get_mp().get_mg()->get_nzone() ;
256  ss.exponential_filter_ylm_phi(0, nz-1, p_r, p_tet, p_phi, alpha) ;
257  return ;
258 }
259 
260 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
friend Scalar exp(const Scalar &)
Exponential.
Definition: scalar_math.C:326
Lorene prototypes.
Definition: app_hor.h:67
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual void exponential_filter_ylm_phi(int lzmin, int lzmax, int p_r, int p_tet, int p_phi, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
void exp_filter_ylm_all_domains(Scalar &ss, int p, double alpha=-16.)
Applies an exponential filter in angular directions in all domains.
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:454
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition: scalar.C:293
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Valeur va
The numerical value of the Scalar.
Definition: scalar.h:411
void sarra_filter_r_all_domains(double p, double alpha=1E-16)
Applies an exponential filter in radial direction in all domains for the case where p is a double (se...
friend Scalar log(const Scalar &)
Neperian logarithm.
Definition: scalar_math.C:388
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
Bases of the spectral expansions.
Definition: base_val.h:325
void sarra_filter_r(int lzmin, int lzmax, double p, double alpha=-1E-16)
Applies an exponential filter to the spectral coefficients in the radial direction.
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Definition: scalar.h:402
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
void del_deriv()
Logical destructor of the derivatives.
Definition: valeur.C:641
friend Scalar abs(const Scalar &)
Absolute value.
Definition: scalar_math.C:526
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
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