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
26char scalar_exp_filter_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $" ;
27
28/*
29 * $Id: scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
30 * $Log: scalar_exp_filter.C,v $
31 * Revision 1.5 2014/10/13 08:53:46 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.4 2014/10/06 15:16:15 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.3 2012/01/17 10:29:27 j_penner
38 * added two routines to handle generalized exponential filtering
39 *
40 * Revision 1.2 2007/10/31 10:50:16 j_novak
41 * Testing whether the coefficients are zero in a given domain.
42 *
43 * Revision 1.1 2007/10/31 10:33:13 j_novak
44 * Added exponential filters to smooth Gibbs-type phenomena.
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_exp_filter.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
48 *
49 */
50
51// C headers
52#include <cassert>
53#include <cmath>
54
55// Lorene headers
56#include "tensor.h"
57#include "proto.h"
58
59namespace Lorene {
61 double alpha) {
62 assert(lzmin >= 0) ;
63 const Mg3d& mgrid = *mp->get_mg() ;
64#ifndef NDEBUG
65 int nz = mgrid.get_nzone() ;
66#endif
67 assert(lzmax < nz) ;
68 assert(etat != ETATNONDEF) ;
69 if (etat == ETATZERO) return ;
70 va.coef() ;
71 assert(va.c_cf != 0x0) ;
72 assert(alpha < 0.) ;
73 double alp = log(pow(10., alpha)) ;
74
75 for (int lz=lzmin; lz<=lzmax; lz++) {
76 if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
77 for (int k=0; k<mgrid.get_np(lz); k++)
78 for (int j=0; j<mgrid.get_nt(lz); j++) {
79 int nr = mgrid.get_nr(lz) ;
80 for (int i=0; i<nr; i++) {
81 double eta = double(i)/double(nr) ;
82 va.c_cf->set(lz, k, j, i) *= exp(alp*pow(eta, 2*p)) ;
83 }
84 }
85 }
86 if (va.c != 0x0) {
87 delete va.c ;
88 va.c = 0x0 ;
89 }
90 va.del_deriv() ;
91 del_deriv() ;
92
93 return ;
94}
95
96void Scalar::sarra_filter_r(int lzmin, int lzmax, double p,
97 double alpha) {
98 assert(lzmin >= 0) ;
99 const Mg3d& mgrid = *mp->get_mg() ;
100#ifndef NDEBUG
101 int nz = mgrid.get_nzone() ;
102#endif
103 assert(lzmax < nz) ;
104 assert(etat != ETATNONDEF) ;
105 if (etat == ETATZERO) return ;
106 va.coef() ;
107 assert(va.c_cf != 0x0) ;
108 assert(alpha < 0.) ;
109
110 for (int lz=lzmin; lz<=lzmax; lz++) {
111 if ((*va.c_cf)(lz).get_etat() == ETATQCQ)
112 for (int k=0; k<mgrid.get_np(lz); k++)
113 for (int j=0; j<mgrid.get_nt(lz); j++) {
114 int nr = mgrid.get_nr(lz) ;
115 for (int i=0; i<nr; i++) {
116 double eta = double(i)/double(nr) ;
117 va.c_cf->set(lz, k, j, i) *= exp(alpha*pow(eta, p)) ;
118 }
119 }
120 }
121 if (va.c != 0x0) {
122 delete va.c ;
123 va.c = 0x0 ;
124 }
125 va.del_deriv() ;
126 del_deriv() ;
127
128 return ;
129}
130void exp_filter_r_all_domains( Scalar& ss, int p, double alpha ) {
131 int nz = ss.get_mp().get_mg()->get_nzone() ;
132 ss.exponential_filter_r(0, nz-1, p, alpha) ;
133 return ;
134}
135
136void Scalar::sarra_filter_r_all_domains( double p, double alpha ) {
137 int nz = get_mp().get_mg()->get_nzone() ;
138 sarra_filter_r(0, nz-1, p, alpha) ;
139 return ;
140}
141
143 double alpha) {
144 assert(lzmin >= 0) ;
145 const Mg3d& mgrid = *mp->get_mg() ;
146#ifndef NDEBUG
147 int nz = mgrid.get_nzone() ;
148#endif
149 assert(lzmax < nz) ;
150 assert(etat != ETATNONDEF) ;
151 if (etat == ETATZERO) return ;
152 double alp = log(pow(10., alpha)) ;
153 va.ylm() ;
154 assert(va.c_cf != 0x0) ;
155 const Base_val& base = va.base ;
156 int l_q, m_q, base_r ;
157
158 for (int lz=lzmin; lz<=lzmax; lz++)
159 if ((*va.c_cf)(lz).get_etat() == ETATQCQ) {
160 int np = mgrid.get_np(lz) ;
161 int nt = mgrid.get_nt(lz) ;
162 int nr = mgrid.get_nr(lz) ;
163 int lmax = base.give_lmax(mgrid, lz) ;
164 for (int k=0; k<np; k++)
165 for (int j=0; j<nt; j++) {
166 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
167 if (nullite_plm(j, nt, k, np, base) == 1 ) {
168 double eta = double(l_q) / double(lmax) ;
169 double sigma = exp(alp*pow(eta, 2*p)) ;
170 for (int i=0; i<nr; i++)
171 va.c_cf->set(lz, k, j, i) *= sigma ;
172 }
173 }
174 }
175
176 va.ylm_i() ;
177 if (va.c != 0x0) {
178 delete va.c ;
179 va.c = 0x0 ;
180 }
181 va.del_deriv() ;
182 del_deriv() ;
183 return ;
184}
185
186void exp_filter_ylm_all_domains(Scalar& ss, int p, double alpha ) {
187 int nz = ss.get_mp().get_mg()->get_nzone() ;
188 ss.exponential_filter_ylm(0, nz-1, p, alpha) ;
189 return ;
190}
191
192}
Bases of the spectral expansions.
Definition base_val.h:322
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Multi-domain grid.
Definition grilles.h:273
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
friend Scalar exp(const Scalar &)
Exponential.
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
friend Scalar log(const Scalar &)
Neperian logarithm.
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.
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.
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
friend Scalar pow(const Scalar &, int)
Power .
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...
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
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 ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void del_deriv()
Logical destructor of the derivatives.
Definition valeur.C:638
void exp_filter_ylm_all_domains(Scalar &ss, int p, double alpha=-16.)
Applies an exponential filter in angular directions in all domains.
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64