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
11namespace 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 prototypes.
Definition app_hor.h:64