LORENE
tensorelliptic.C
1 #include<math.h>
2 // Lorene headers
3 #include "metric.h"
4 #include "nbr_spx.h"
5 #include "utilitaires.h"
6 #include "graphique.h"
7 #include "proto.h"
8 #include "diff.h"
9 
10 
11 namespace Lorene {
12  void tensorelliptic( Scalar source, Scalar& resu, double fit, double fit2, double fit0d2, double fit1d2, double fit0d3, double fit1d3) {
13 
14  const int nz = (*source.get_mp().get_mg()).get_nzone(); // Number of domains
15  int nr = (*source.get_mp().get_mg()).get_nr(1); // Number of collocation points in r in each domain
16  int nt = (*source.get_mp().get_mg()).get_nt(1); // Number of collocation points in theta in each domain
17  int np = (*source.get_mp().get_mg()).get_np(1); // Number of collocation points in phi in each domain
18  const Map_af* map = dynamic_cast<const Map_af*>(&source.get_mp()) ;
19  const Mg3d* mgrid = (*map).get_mg();
20 
21  // Some helpful stuff...
22  const Coord& rr = (*map).r;
23  Scalar rrr (*map) ;
24  rrr = rr ;
25  rrr.set_spectral_va().set_base(source.get_spectral_va().base);
26 
27  Scalar source_coq = source ;
28  source_coq.mult_r() ;
29  source_coq.mult_r() ;
30  source.set_spectral_va().ylm() ;
31  source_coq.set_spectral_va().ylm() ;
32  Scalar phi(source.get_mp()) ;
33  phi.annule_hard() ;
34  // phi.std_spectral_base();
35  phi.set_spectral_va().set_base(source.get_spectral_va().base) ;
36  phi.set_spectral_va().ylm() ;
37  Mtbl_cf& sol_coef = (*phi.set_spectral_va().c_cf) ;
38 
39  const Base_val& base = source.get_spectral_base() ;
40  Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
41  Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
42  Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
43 
44  int l_q, m_q, base_r ;
45 
46  { int lz = 0 ;
47  for (int k=0; k < np; k++)
48  for (int j=0; j<nt; j++) {
49  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
50  if (nullite_plm(j, nt, k, np, base) == 1) {
51  for (int ii=0 ; ii<nr ; ii++){
52  sol_hom1.set(lz, k, j, ii) = 0 ;
53  sol_part.set(lz, k, j, ii) = 0 ;
54  }
55 
56  }
57  }
58  }
59 
60 
61  { int lz = 1 ; // The first shell is a really particular case, where the operator is different, and homogeneous solutions have to be handled really carefully.
62  double alpha = (*map).get_alpha()[lz] ;
63  double beta = (*map).get_beta()[lz] ;
64  double ech = beta / alpha ;
65  for (int k=0; k < np; k++)
66  for (int j=0; j<nt; j++) {
67  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
68  if (nullite_plm(j, nt, k, np, base) == 1) {
69 
70 
71 
72  Matrice ope(nr,nr) ;
73  ope.annule_hard() ;
74 
75  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
76  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
77  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
78  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
79  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
80  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
81  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
82  Diff_x4dsdx2 x4dx2 (base_r, nr); const Matrice& mx4dx2 = x4dx2.get_matrice();
83 
84 
85  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
86  - l_q*(l_q+1)*mid - ((mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) -(fit)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2)
87  + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
88  + fit2*alpha*(beta-1.)*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*(beta- rrr.val_grid_point(1,0,0, nr-1))*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2)+ fit2*(beta-1.)*(beta- rrr.val_grid_point(1,0,0,nr -1))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2));
89 
90 
91 
92 
93  for (int col=0; col<nr; col++)
94  ope.set(nr-1, col) = 0 ;
95  ope.set(nr-1, 0) = 1 ;
96 
97 
98  Tbl rhs(nr);
99  rhs.annule_hard() ;
100  for (int i=0; i<nr; i++)
101  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(1, k, j, i) ;
102  rhs.set(nr-1) = 0 ;
103  ope.set_lu() ;
104  Tbl sol = ope.inverse(rhs) ;
105 
106 
107  for (int i=0; i<nr; i++)
108  sol_part.set(1, k, j, i) = sol(i) ;
109 
110  rhs.annule_hard();
111  rhs.set(nr-1) = 1 ;
112  sol = ope.inverse(rhs) ;
113 
114 
115  for (int i=0; i<nr; i++)
116  sol_hom1.set(1, k, j, i) = sol(i) ;
117 
118 
119  }
120  }
121  }
122 
123 
124  // Current implementations only allow grids with more than 3 shells. Zones 2 and 3 are also treated separately, allowing
125  // A smoother interpolation of the real operator;
126 
127  { int lz = 2 ;
128 
129  double alpha = (*map).get_alpha()[lz] ;
130  double beta = (*map).get_beta()[lz] ;
131  double ech = beta / alpha ;
132 
133  for (int k=0; k < np; k++)
134  for (int j=0; j<nt; j++) {
135  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
136  if (nullite_plm(j, nt, k, np, base) == 1) {
137 
138  Matrice ope(nr,nr) ;
139  ope.annule_hard() ;
140 
141  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
142  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
143  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
144  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
145  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
146  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
147  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
148 
149  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
150  - l_q*(l_q+1)*mid - fit0d2*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit1d2)*(rrr.val_grid_point(lz, 0, 0, 0))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d2)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d2)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) ;
151 
152 
153 
154  for (int col=0; col<nr; col++)
155  ope.set(nr-1, col) = 0 ;
156  ope.set(nr-1, 0) = 1 ;
157  for (int col=0; col<nr; col++) {
158  ope.set(nr-2, col) = 0 ;
159  }
160  ope.set(nr-2, 1) = 1 ;
161 
162  Tbl rhs(nr) ;
163  rhs.annule_hard() ;
164  for (int i=0; i<nr; i++)
165  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
166  rhs.set(nr-2) = 0 ;
167  rhs.set(nr-1) = 0 ;
168  ope.set_lu() ;
169  Tbl sol = ope.inverse(rhs) ;
170 
171  for (int i=0; i<nr; i++)
172  sol_part.set(lz, k, j, i) = sol(i) ;
173 
174  rhs.annule_hard() ;
175  rhs.set(nr-2) = 1 ;
176  sol = ope.inverse(rhs) ;
177  for (int i=0; i<nr; i++)
178  sol_hom1.set(lz, k, j, i) = sol(i) ;
179 
180  rhs.set(nr-2) = 0 ;
181  rhs.set(nr-1) = 1 ;
182  sol = ope.inverse(rhs) ;
183  for (int i=0; i<nr; i++)
184  sol_hom2.set(lz, k, j, i) = sol(i) ;
185 
186  }
187  }
188 
189  }
190 
191 
192 
193  { int lz = 3 ;
194 
195  double alpha = (*map).get_alpha()[lz] ;
196  double beta = (*map).get_beta()[lz] ;
197  double ech = beta / alpha ;
198 
199  for (int k=0; k < np; k++)
200  for (int j=0; j<nt; j++) {
201  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
202  if (nullite_plm(j, nt, k, np, base) == 1) {
203 
204  Matrice ope(nr,nr) ;
205  ope.annule_hard() ;
206 
207  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
208  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
209  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
210  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
211  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
212  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
213  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
214 
215  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
216  - l_q*(l_q+1)*mid - fit0d3*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit1d3)*(rrr.val_grid_point(lz, 0, 0, 0))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d3)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d3)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) ;
217 
218  for (int col=0; col<nr; col++)
219  ope.set(nr-1, col) = 0 ;
220  ope.set(nr-1, 0) = 1 ;
221  for (int col=0; col<nr; col++) {
222  ope.set(nr-2, col) = 0 ;
223  }
224  ope.set(nr-2, 1) = 1 ;
225 
226  Tbl rhs(nr) ;
227  rhs.annule_hard() ;
228  for (int i=0; i<nr; i++)
229  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
230  rhs.set(nr-2) = 0 ;
231  rhs.set(nr-1) = 0 ;
232  ope.set_lu() ;
233  Tbl sol = ope.inverse(rhs) ;
234 
235  for (int i=0; i<nr; i++)
236  sol_part.set(lz, k, j, i) = sol(i) ;
237 
238  rhs.annule_hard() ;
239  rhs.set(nr-2) = 1 ;
240  sol = ope.inverse(rhs) ;
241  for (int i=0; i<nr; i++)
242  sol_hom1.set(lz, k, j, i) = sol(i) ;
243 
244  rhs.set(nr-2) = 0 ;
245  rhs.set(nr-1) = 1 ;
246  sol = ope.inverse(rhs) ;
247  for (int i=0; i<nr; i++)
248  sol_hom2.set(lz, k, j, i) = sol(i) ;
249 
250  }
251  }
252 
253  }
254 
255 
256 
257 
258  for (int lz=4; lz<nz-1; lz++) {
259  double ech = (*map).get_beta()[lz] / (*map).get_alpha()[lz] ;
260  for (int k=0; k < np; k++)
261  for (int j=0; j<nt; j++) {
262  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
263  if (nullite_plm(j, nt, k, np, base) == 1) {
264 
265  Matrice ope(nr,nr) ;
266  ope.annule_hard() ;
267 
268  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
269  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
270  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
271  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
272  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
273  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
274  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
275  - l_q*(l_q+1)*mid ;
276 
277  for (int col=0; col<nr; col++)
278  ope.set(nr-1, col) = 0 ;
279  ope.set(nr-1, 0) = 1 ;
280  for (int col=0; col<nr; col++) {
281  ope.set(nr-2, col) = 0 ;
282  }
283  ope.set(nr-2, 1) = 1 ;
284 
285  Tbl rhs(nr) ;
286  rhs.annule_hard() ;
287  for (int i=0; i<nr; i++)
288  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
289  rhs.set(nr-2) = 0 ;
290  rhs.set(nr-1) = 0 ;
291  ope.set_lu() ;
292  Tbl sol = ope.inverse(rhs) ;
293 
294  for (int i=0; i<nr; i++)
295  sol_part.set(lz, k, j, i) = sol(i) ;
296 
297  rhs.annule_hard() ;
298  rhs.set(nr-2) = 1 ;
299  sol = ope.inverse(rhs) ;
300  for (int i=0; i<nr; i++)
301  sol_hom1.set(lz, k, j, i) = sol(i) ;
302 
303  rhs.set(nr-2) = 0 ;
304  rhs.set(nr-1) = 1 ;
305  sol = ope.inverse(rhs) ;
306  for (int i=0; i<nr; i++)
307  sol_hom2.set(lz, k, j, i) = sol(i) ;
308 
309  }
310  }
311 
312  }
313  { int lz = nz-1 ;
314  double alpha = (*map).get_alpha()[lz] ;
315  for (int k=0; k < np; k++)
316  for (int j=0; j<nt; j++) {
317  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
318  if (nullite_plm(j, nt, k, np, base) == 1) {
319 
320  Matrice ope(nr,nr) ;
321  ope.annule_hard() ;
322  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
323  Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
324 
325  ope = (mdx2 - l_q*(l_q+1)*ms2)/(alpha*alpha) ;
326 
327  for (int i=0; i<nr; i++)
328  ope.set(nr-1, i) = 0 ;
329  ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
330 
331  for (int i=0; i<nr; i++) {
332  ope.set(nr-2, i) = 1 ; //for the limit at infinity
333  }
334 
335  if (l_q > 0) {
336  for (int i=0; i<nr; i++) {
337  ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
338  }
339  }
340 
341  Tbl rhs(nr) ;
342  rhs.annule_hard() ;
343  for (int i=0; i<nr; i++)
344  rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, k, j, i) ;
345  if (l_q>0) rhs.set(nr-3) = 0 ;
346  rhs.set(nr-2) = 0 ;
347  rhs.set(nr-1) = 0 ;
348  ope.set_lu() ;
349  Tbl sol = ope.inverse(rhs) ;
350 
351  for (int i=0; i<nr; i++)
352  sol_part.set(lz, k, j, i) = sol(i) ;
353 
354  rhs.annule_hard() ;
355  rhs.set(nr-1) = 1 ;
356  sol = ope.inverse(rhs) ;
357  for (int i=0; i<nr; i++)
358  sol_hom2.set(lz, k, j, i) = sol(i) ;
359 
360  }
361  }
362  }
363 
364  Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
365  Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
366  Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
367 
368 
369  // Now matching the homogeneous solutions between the different domains...
370 
371  for (int k=0; k < np; k++)
372  for (int j=0; j<nt; j++) {
373  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
374  if (nullite_plm(j, nt, k, np, base) == 1) {
375  Matrice systeme(2*nz-4, 2*nz-4) ;
376  systeme.annule_hard() ;
377  Tbl rhs(2*nz-4) ;
378  rhs.annule_hard() ;
379 
380  //First shell
381  int lin = 0 ;
382  int col = 0 ;
383 
384  double alpha = (*map).get_alpha()[1] ;
385 
386  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(1, j, k) ;
387  rhs.set(lin) -= sol_part.val_out_bound_jk(1, j, k) ;
388 
389  lin++ ;
390  systeme.set(lin, col) += dhom1.val_out_bound_jk(1, j, k) / alpha ;
391  rhs.set(lin) -= dpart.val_out_bound_jk(1, j, k) / alpha ;
392  col += 1 ;
393 
394  //Shells
395  for (int lz=2; lz<nz-1; lz++) {
396  alpha = (*map).get_alpha()[lz] ;
397  lin-- ;
398  systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, j, k) ;
399  systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, j, k) ;
400  rhs.set(lin) += sol_part.val_in_bound_jk(lz, j, k) ;
401 
402  lin++ ;
403  systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, j, k) / alpha ;
404  systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, j, k) / alpha ;
405  rhs.set(lin) += dpart.val_in_bound_jk(lz, j, k) / alpha;
406 
407  lin++ ;
408  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, j, k) ;
409  systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, j, k) ;
410  rhs.set(lin) -= sol_part.val_out_bound_jk(lz, j, k) ;
411 
412  lin++ ;
413  systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, j, k) / alpha ;
414  systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, j, k) / alpha ;
415  rhs.set(lin) -= dpart.val_out_bound_jk(lz, j, k) / alpha ;
416  col += 2 ;
417  }
418 
419  //CED
420  alpha = (*map).get_alpha()[nz-1] ;
421  lin-- ;
422  systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, j, k) ;
423  rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, j, k) ;
424 
425  lin++ ;
426  systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, j, k) ;
427  rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, j, k) ;
428 
429  systeme.set_lu() ;
430 
431  Tbl coef = systeme.inverse(rhs);
432  int indice = 0 ;
433 
434 
435  for (int i=0; i<(*mgrid).get_nr(1); i++)
436  sol_coef.set(1, k, j, i) = sol_part(1, k, j, i)
437  +coef(indice)*sol_hom1(1, k, j, i) ;
438  indice +=1;
439 
440 
441  for (int lz=2; lz<nz-1; lz++) {
442  for (int i=0; i<(*mgrid).get_nr(lz); i++)
443  sol_coef.set(lz, k, j, i) = sol_part(lz, k, j, i)
444  +coef(indice)*sol_hom1(lz, k, j, i)
445  +coef(indice+1)*sol_hom2(lz, k, j, i) ;
446  indice += 2 ;
447  }
448  for (int i=0; i<(*mgrid).get_nr(nz-1); i++)
449  sol_coef.set(nz-1, k, j, i) = sol_part(nz-1, k, j, i)
450  +coef(indice)*sol_hom2(nz-1, k, j, i) ;
451 
452  }
453  }
454 
455 
456  delete phi.set_spectral_va().c ;
457  phi.set_spectral_va().c = 0x0 ;
458  phi.set_spectral_va().ylm_i() ;
459 
460  phi.annule_domain(nz-1);
461 
462  resu = phi;
463 
464 }
465 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67