LORENE
map_interp_af_and_star.C
1 /*
2  * Copyright (c) 2024 Jerome Novak
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 /*
24  * Methods of classes Map_af and Map_star to interpolate between one another.
25  * (see files map.h and map_star.h)
26  */
27 
28 #include "tensor.h"
29 #include "param.h"
30 #include "proto.h"
31 #include "utilitaires.h"
32 
33 using namespace Lorene ;
34 
35  //---------------------------------------------------------------
36  // Static variables : arrays of coefficient/summation functions
37  //---------------------------------------------------------------
38 namespace
39 {
40  bool first_call = true ;
41 
42  void (*coef_r[MAX_BASE])(const int*, const int*, double*, const int*, double*) ;
43  double (*som_r_1d[MAX_BASE])(const double*, int, double) ;
44 }
45 
46 #include "proto_interp.h"
47 
48  //---------------------------------------------
49  // Spectral summation from Map_af to Map_star
50  //---------------------------------------------
51 
53 {
54 
55  //## Only valid for TYPE_T = SYM & TYPE_P = SYM for the moment (to be improved)
56  assert((mg->get_type_t() == SYM) && (mg->get_type_p() == SYM) ) ;
57 
58  // First call operations
59  if (first_call) {
60  do_first_call_initializations(coef_r, som_r_1d) ;
61  first_call = false ;
62  }
63 
64  // Check whether the input Scalar is defined on a Map_af
65  const Map* p_mp = &f_a.get_mp() ;
66  const Map_af* p_mpa = dynamic_cast<const Map_af*>(p_mp) ;
67  if (p_mpa == 0x0)
68  throw(invalid_argument("Map_star::interpolate_from_map_af: the input Scalar field is not of type Map_af.")) ;
69 
70  const Map_af& mpa = *p_mpa ;
71  const Mg3d* mg_a = mpa.get_mg() ;
72 
73 
74  // Checks that both grids are compatible... (## add symmetries!!)
75  int ndom_a = 0 ;
76  if (!check_grids(mpa, *this, ndom_a)) {
77  throw(invalid_argument("Map_star::interpolate_from_map_af: both mappings are not compatible")) ;
78  }
79  int np0 = mg_a->get_np(0) ;
80  int nt0 = mg_a->get_nt(0) ;
81 
82  Scalar resu(*this) ;
83 
84  switch (f_a.get_etat()) {
85 
86  case ETATNONDEF: {
87  throw(invalid_argument("Map_star::interpolate_from_map_af: the input Scalar field is not defined.")) ;
88  break;
89  }
90 
91  case ETATZERO:{
92  resu.set_etat_zero() ;
93  break ;
94  }
95 
96  case ETATUN:{
97  resu.set_etat_one() ;
98  break ;
99  }
100 
101  case ETATQCQ:{ // General case
102  resu.allocate_all() ;
103  const Base_val& base = f_a.get_spectral_va().base ;
104 
105  //Intermediate array: coefficients in r, values in \theta & \varphi
106  Mtbl interm(mg_a) ;
107  interm.annule_hard() ;
108 
109  // Compute function values # Improve with a transform in theta, phi only
110  f_a.get_spectral_va().coef_i() ;
111 
112  // Makes the coefficient transform in the r variable only
113  //-------------------------------------------------------
114  for (int l=0; l<ndom_a; l++) {
115  const Tbl* f = ((f_a.get_spectral_va().c)->t)[l] ;
116  Tbl* cf = (interm.t)[l];
117  if (f->get_etat() == ETATZERO) {
118  cf->set_etat_zero() ;
119  continue ; // nothing to do if the Tbl = 0
120  }
121  int nr = f->get_dim(0) ;
122  int deg[3] ;
123  deg[0] = np0 ;
124  deg[1] = nt0 ;
125  deg[2] = nr ;
126 
127  // Takes information on the r-basis
128  int base_r = ( base.b[l] & MSQ_R ) >> TRA_R ;
129  assert(base_r < MAX_BASE) ;
130 
131  // Transformation in r:
132  // --------------------
133  if ( nr > 1 ) {
134  assert( admissible_fft(nr-1) || (mg->get_colloc_r(l) != BASE_CHEB) );
135  coef_r[base_r](deg, deg, (f->t), deg, (cf->t)) ;
136  }
137  } // end of loop on domains
138 
139  // Computing values on grid points associated to Map_star coordinates
140  //-------------------------------------------------------------------
141  int nzs = get_mg()->get_nzone() ;
142  Param par_dummy ;
143 
144  Mtbl rad_af = (+r) ;
145  for (int lz=0; lz<nzs; lz++) {
146  for (int k=0; k<np0; k++) {
147  for (int j=0; j<nt0; j++) {
148  int nrs = get_mg()->get_nr(lz) ;
149  for (int i=0; i<nrs; i++) {
150  double radius = rad_af(lz, k, j, i) ;
151  int l_a = 0 ;
152  double xi_a = 0. ;
153  mpa.val_lx_jk(radius, j, k, par_dummy, l_a, xi_a) ;
154  int base_r = ( base.b[l_a] & MSQ_R ) >> TRA_R ;
155  const double* cfr = &interm.set(l_a, k, j, 0) ;
156  int nra = mg_a->get_nr(l_a) ;
157 
158  // Call to the 'som_r_1d' function that makes the spectral summation
159  resu.set_grid_point(lz, k, j, i) = som_r_1d[base_r](cfr, nra, xi_a) ;
160  }
161  }
162  }
163  }
164  break ;
165  } // end of case ETATQCQ
166  default: {
167  throw(invalid_argument("Map_star::interpolate_from_map_af: the input Scalar field is ill-formed.")) ;
168  break ;
169  }
170  } // end of switch
171 
172  return resu ;
173 
174 }
175 
176  //---------------------------------------------
177  // Spectral summation from Map_star to Map_af
178  //---------------------------------------------
179 
181 {
182 
183  //## Only valid for TYPE_T = SYM & TYPE_P = SYM for the moment (to be improved)
184  assert((mg->get_type_t() == SYM) && (mg->get_type_p() == SYM) ) ;
185 
186  // First call operations
187  if (first_call) {
188  do_first_call_initializations(coef_r, som_r_1d) ;
189  first_call = false ;
190  }
191 
192  // Check whether the input Scalar is defined on a Map_star
193  const Map* p_mp = &f_s.get_mp() ;
194  const Map_star* p_mps = dynamic_cast<const Map_star*>(p_mp) ;
195  if (p_mps == 0x0)
196  throw(invalid_argument("Map_af::interpolate_from_map_star: the input Scalar field is not of type Map_star.")) ;
197 
198  const Map_star& mps = *p_mps ;
199  const Mg3d* mg_s = mps.get_mg() ;
200 
201  // Checks that both grids are compatible...
202  int ndom_a = 0 ;
203  if (!check_grids(mps, *this, ndom_a)) {
204  throw(invalid_argument("Map_star::interpolate_from_map_star: both mappings are not compatible")) ;
205  }
206  int np0 = mg_s->get_np(0) ;
207  int nt0 = mg_s->get_nt(0) ;
208 
209  Scalar resu(*this) ;
210 
211  switch (f_s.get_etat()) {
212 
213  case ETATNONDEF: {
214  throw(invalid_argument("Map_af::interpolate_from_map_star: the input Scalar field is not defined.")) ;
215  break;
216  }
217 
218  case ETATZERO:{
219  resu.set_etat_zero() ;
220  break ;
221  }
222 
223  case ETATUN:{
224  resu.set_etat_one() ;
225  break ;
226  }
227 
228  case ETATQCQ:{ // General case
229  resu.allocate_all() ;
230  const Base_val& base = f_s.get_spectral_va().base ;
231 
232  //Intermediate array: coefficients in r, values in \theta & \varphi
233  Mtbl interm(mg_s) ;
234  interm.annule_hard() ;
235 
236  // Compute function values # Improve with a transform in theta, phi only
237  f_s.get_spectral_va().coef_i() ;
238 
239  // Makes the coefficient transform in the r variable only
240  //-------------------------------------------------------
241  int nzs = mg_s->get_nzone() ;
242  for (int l=0; l<nzs; l++) {
243  const Tbl* f = ((f_s.get_spectral_va().c)->t)[l] ;
244  Tbl* cf = (interm.t)[l];
245  if (f->get_etat() == ETATZERO) {
246  cf->set_etat_zero() ;
247  continue ; // nothing to do if the Tbl = 0
248  }
249  int nr = f->get_dim(0) ;
250  int deg[3] ;
251  deg[0] = np0 ;
252  deg[1] = nt0 ;
253  deg[2] = nr ;
254 
255  // Takes information on the r-basis
256  int base_r = ( base.b[l] & MSQ_R ) >> TRA_R ;
257  assert(base_r < MAX_BASE) ;
258 
259  // Transformation in r:
260  // --------------------
261  if ( nr > 1 ) {
262  assert( admissible_fft(nr-1) || (mg->get_colloc_r(l) != BASE_CHEB) );
263  coef_r[base_r](deg, deg, (f->t), deg, (cf->t)) ;
264  }
265  } // end of loop on domains
266 
267  // Computing values on grid points associated to Map_star coordinates
268  //-------------------------------------------------------------------
269  int nza = get_mg()->get_nzone() ;
270  Param par_dummy ;
271  const Coord& r_s = mps.r ;
272  Mtbl rad_star = (+r_s) ;
273  Mtbl rad_af = (+r) ;
274  int nrs_max = mg_s->get_nr(nzs-1) ;
275  for (int lz=0; lz<ndom_a; lz++) {
276  for (int k=0; k<np0; k++) {
277  for (int j=0; j<nt0; j++) {
278  double rmax_s = rad_star(nzs-1, k, j, nrs_max-1) ;
279  int nra = get_mg()->get_nr(lz) ;
280  for (int i=0; i<nra; i++) {
281  double radius_a = rad_af(lz, k, j, i) ;
282  if (radius_a <= rmax_s) {
283  int l_s = 0 ;
284  double xi_s = 0. ;
285  mps.val_lx_jk(radius_a, j, k, par_dummy, l_s, xi_s) ;
286  int base_r = ( base.b[l_s] & MSQ_R ) >> TRA_R ;
287  const double* cfr = &interm.set(l_s, k, j, 0) ;
288  int nrs = mg_s->get_nr(l_s) ;
289 
290  // Call to the 'som_r_1d' function that makes the spectral summation
291  resu.set_grid_point(lz, k, j, i) = som_r_1d[base_r](cfr, nrs, xi_s) ;
292  }
293  else // the point is outside the Map_star grid
294  resu.set_grid_point(lz, k, j, i) = 0. ;
295  }
296  }
297  }
298  }
299  resu.annule(ndom_a, nza-1) ;// Zero outside the domains containing the Map_star
300  break ;
301  } // end of case ETATQCQ
302  default: {
303  throw(invalid_argument("Map_af::interpolate_from_map_star: the input Scalar field is ill-formed.")) ;
304  break ;
305  }
306  } // end of switch
307 
308  return resu ;
309 
310 }
311 
312 namespace Lorene {
313  //------------------------
314  // Miscellaneous
315  //------------------------
316 
317 void do_first_call_initializations(
318  void (*coef_r0[MAX_BASE])
319  (const int*, const int*, double*, const int*, double*),
320  double (*som_r_1d0[MAX_BASE])(const double*, int, double) ) {
321 
322  for (int i=0; i<MAX_BASE; i++) {
323  coef_r0[i] = pasprevu_r ;
324  som_r_1d0[i] = som_r_1d_pas_prevu ;
325  }
326  coef_r0[NONDEF] = base_non_def_r ;
327  coef_r0[R_CHEB >> TRA_R] = cfrcheb_interp ;
328  coef_r0[R_CHEBU >> TRA_R] = cfrcheb_interp ;
329  coef_r0[R_CHEBP >> TRA_R] = cfrchebp_interp ;
330  coef_r0[R_CHEBI >> TRA_R] = cfrchebi_interp ;
331  // coef_r0[R_LEG >> TRA_R] = cfrleg ;
332  // coef_r0[R_LEGP >> TRA_R] = cfrlegp ;
333  // coef_r0[R_LEGI >> TRA_R] = cfrlegi ;
334  // coef_r0[R_JACO02 >> TRA_R] = cfrjaco02 ;
335 
336  som_r_1d0[R_CHEB >> TRA_R] = som_r_1d_cheb ;
337  som_r_1d0[R_CHEBP >> TRA_R] = som_r_1d_chebp ;
338  som_r_1d0[R_CHEBI >> TRA_R] = som_r_1d_chebi ;
339  som_r_1d0[R_CHEBU >> TRA_R] = som_r_1d_cheb ; // no error here...
340  som_r_1d0[R_LEG >> TRA_R] = som_r_1d_leg ;
341  som_r_1d0[R_LEGP >> TRA_R] = som_r_1d_legp ;
342  som_r_1d0[R_LEGI >> TRA_R] = som_r_1d_legi ;
343  som_r_1d0[R_JACO02 >> TRA_R] = som_r_1d_jaco02 ;
344 }
345 
346 bool check_grids(const Map_af& mpa, const Map_star& mps, int& ndom) {
347 
348  const Mg3d* mga = mpa.get_mg() ;
349  const Mg3d* mgs = mps.get_mg() ;
350  int nza = mga->get_nzone() ;
351  int nzs = mgs->get_nzone() ;
352 
353  // First, check that both grids have same symmetries.
354  bool resu = ( mga->get_type_t() == mgs->get_type_t() )
355  && ( mga->get_type_p() == mgs->get_type_p() ) ;
356 
357  // Then, check that the number of points in theta & phi are the same.
358  int np0 = mga->get_np(0) ;
359  int nt0 = mga->get_nt(0) ;
360  for (int l=1; l<nza; l++)
361  resu = resu && (np0 == mga->get_np(l)) && (nt0 == mga->get_nt(l)) ;
362  for (int l=0; l<nzs; l++)
363  resu = resu && (np0 == mgs->get_np(l)) && (nt0 == mgs->get_nt(l)) ;
364 
365  // Finally, looks at the maximal radius of Map_star ...
366  const Coord& rr = mps.r ;
367  Mtbl rad_star = (+rr) ;
368  double rmax = 0. ;
369  int jmax = -1 ;
370  int kmax = -1 ;
371  int nr = mgs->get_nr(nzs-1) ;
372  for (int k=0; k<np0; k++) {
373  for (int j=0; j<nt0; j++) {
374  if (rad_star(nzs-1, k, j, nr-1) > rmax) {
375  rmax = rad_star(nzs-1, k, j, nr-1) ;
376  jmax = j ;
377  kmax = k ;
378  }
379  }
380  }
381  // ... and the corresponding domain for Map_af.
382  double dummy = 0. ;
383  Param par_dummy ;
384  mpa.val_lx_jk(rmax, jmax, kmax, par_dummy, ndom, dummy) ;
385  ndom++ ;
386 
387  return resu ;
388 }
389 
390 bool check_grids(const Map_star& mps, const Map_af& mpa, int& ndom) {
391 
392  const Mg3d* mga = mpa.get_mg() ;
393  const Mg3d* mgs = mps.get_mg() ;
394  int nza = mga->get_nzone() ;
395  int nzs = mgs->get_nzone() ;
396 
397  // First, check that both grids have same symmetries.
398  bool resu = ( mga->get_type_t() == mgs->get_type_t() )
399  && ( mga->get_type_p() == mgs->get_type_p() ) ;
400 
401  // Then, check that the number of points in theta & phi are the same.
402  int np0 = mga->get_np(0) ;
403  int nt0 = mga->get_nt(0) ;
404  for (int l=1; l<nza; l++)
405  resu = resu && (np0 == mga->get_np(l)) && (nt0 == mga->get_nt(l)) ;
406  for (int l=0; l<nzs; l++)
407  resu = resu && (np0 == mgs->get_np(l)) && (nt0 == mgs->get_nt(l)) ;
408 
409  // Finally, looks at the maximal radius of Map_star ...
410  const Coord& rr = mps.r ;
411  Mtbl rad_star = (+rr) ;
412  double rmax = 0. ;
413  int jmax = -1 ;
414  int kmax = -1 ;
415  int nr = mgs->get_nr(nzs-1) ;
416  for (int k=0; k<np0; k++) {
417  for (int j=0; j<nt0; j++) {
418  if (rad_star(nzs-1, k, j, nr-1) > rmax) {
419  rmax = rad_star(nzs-1, k, j, nr-1) ;
420  jmax = j ;
421  kmax = k ;
422  }
423  }
424  }
425  // ... and the corresponding domain for Map_af.
426  double dummy = 0. ;
427  Param par_dummy ;
428  mpa.val_lx_jk(rmax, jmax, kmax, par_dummy, ndom, dummy) ;
429  ndom++ ;
430 
431  return resu ;
432 }
433 
434 }
435 
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:512
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:789
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
void coef_i() const
Computes the physical value of *this.
Base class for coordinate mappings.
Definition: map.h:694
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
int get_colloc_r(int l) const
Returns the type of collocation points used in domain no.
Definition: grilles.h:528
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
double * t
The array of double.
Definition: tbl.h:176
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:334
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
Parameter storage.
Definition: param.h:125
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition: scalar.C:340
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition: mtbl.h:221
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:350
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
void annule_hard()
Sets the Mtbl to zero in a hard way.
Definition: mtbl.C:315
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:700
Multi-domain grid.
Definition: grilles.h:279
Bases of the spectral expansions.
Definition: base_val.h:325
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
Affine radial mapping.
Definition: map.h:2087
#define NONDEF
base inconnue
Definition: type_parite.h:150
bool admissible_fft(int)
Checks whether or not a given number of degrees of freedom is compatible with the FFT transform...
Scalar interpolate_from_map_star(const Scalar &f_s) const
Interpolates from a Scalar defined on a Map_star and returns the new Scalar defined on *this...
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Scalar interpolate_from_map_af(const Scalar &f_a) const
Interpolates from a Scalar defined on a Map_af and returns the new Scalar defined on *this...
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Coord r
r coordinate centered on the grid
Definition: map.h:742
#define R_LEG
base de Legendre ordinaire (fin)
Definition: type_parite.h:182
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166