LORENE
vector_df_poisson.C
1 /*
2  * Resolution of the divergence-free vector Poisson equation
3  *
4  * (see file vector.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: vector_df_poisson.C,v 1.16 2016/12/05 16:18:18 j_novak Exp $
32  * $Log: vector_df_poisson.C,v $
33  * Revision 1.16 2016/12/05 16:18:18 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.15 2014/10/13 08:53:44 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.14 2014/10/06 15:13:20 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.13 2005/02/15 09:45:22 j_novak
43  * Correction of an error in the matching.
44  *
45  * Revision 1.12 2005/02/09 16:53:11 j_novak
46  * Now V^r and eta are matched across domains, but not any of their derivatives.
47  *
48  * Revision 1.11 2005/02/09 14:52:01 j_novak
49  * Better solution in the shells.
50  *
51  * Revision 1.10 2005/02/09 13:20:27 j_novak
52  * Completely new way of solving the vector Poisson equation in the Div_free
53  * case: the system is inverted "as a whole" for V^r and eta. This works only
54  * with Map_af...
55  *
56  *
57  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.16 2016/12/05 16:18:18 j_novak Exp $
58  *
59  */
60 
61 // C headers
62 #include <cassert>
63 #include <cmath>
64 
65 // Lorene headers
66 #include "tensor.h"
67 #include "diff.h"
68 #include "proto.h"
69 #include "param.h"
70 
71 namespace Lorene {
73 
74  // All this has a meaning only for spherical components:
75 #ifndef NDEBUG
76  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
77  assert(bvs != 0x0) ;
78 #endif
79 
80  int nitermax = par.get_int() ;
81  int& niter = par.get_int_mod() ;
82  double relax = par.get_double() ;
83  double precis = par.get_double(1) ;
84  Cmp& ss_khi = par.get_cmp_mod(0) ;
85  Cmp& ss_mu = par.get_cmp_mod(1) ;
86 
87  // Solution for the r-component
88  // ----------------------------
89 
90  Scalar source_r = *(cmp[0]) ;
91  source_r.mult_r() ;
92 
93  Param par_khi ;
94  par_khi.add_int(nitermax, 0) ;
95  par_khi.add_int_mod(niter, 0) ;
96  par_khi.add_double(relax, 0) ;
97  par_khi.add_double(precis, 1) ;
98  par_khi.add_cmp_mod(ss_khi, 0) ;
99 
100  Scalar khi (*mp) ;
101  khi.set_etat_zero() ;
102 
103  source_r.poisson(par_khi, khi) ;
104  khi.div_r() ; // khi now contains V^r
105 
106  // Solution for mu
107  // ---------------
108 
109  Param par_mu ;
110  par_mu.add_int(nitermax, 0) ;
111  par_mu.add_int_mod(niter, 0) ;
112  par_mu.add_double(relax, 0) ;
113  par_mu.add_double(precis, 1) ;
114  par_mu.add_cmp_mod(ss_mu, 0) ;
115 
116  Scalar mu_resu (*mp) ;
117  mu_resu.set_etat_zero() ;
118 
119  mu().poisson(par_mu, mu_resu) ;
120 
121  // Final result
122  // ------------
123 
124  Vector_divfree resu(*mp, *triad, *met_div) ;
125 
126  resu.set_vr_mu(khi, mu_resu) ;
127 
128  return resu ;
129 
130 }
131 
132 /*
133  * In the case without parameters, first is solved the equation for mu and then
134  * the system of equations for (eta, V^r) is inverted as a whole:
135  * d2 eta / dr2 + 2/r d eta / dr - 1/r d V^r / dr = S(eta)
136  * d V^r / dr + 2/r V^r - l(l+1)/r eta = 0 (div free condition)
137  *
138  * There is no l=0 contribution (divergence free in all space!).
139  * In the nucleus and the CED the system is inverted for h(r) and v(r) ,
140  * such that eta = r^2 h and V^r = r^2 v in the nucleus,
141  * in the compactified domain one has eta = u^2 h and V^r = u^2 v (where u=1/r);
142  * In the shells, both equations are multiplied by r.
143  * These methods are used only to get particular solutions.
144  *
145  * Homogeneous solutions are known analitically: r^(l-1) and/or 1/r^(l+2)
146  * It is then only possible to match eta and V^r, but not their derivatives,
147  * due to the div-free condition.
148  */
150 
151  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
152 #ifndef NDEBUG
153  for (int i=0; i<3; i++)
154  assert(cmp[i]->check_dzpuis(4)) ;
155  // All this has a meaning only for spherical components:
156  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
157  assert(bvs != 0x0) ;
158  //## ... and affine mapping, for the moment!
159  assert( mpaff != 0x0) ;
160 #endif
161 
162  // Final result
163  // ------------
164  Vector_divfree resu(*mpaff, *triad, *met_div) ;
165 
166  // Solution for mu
167  // ---------------
168  Scalar mu_resu = mu().poisson() ;
169 
170  Scalar f_r(*mpaff) ;
171  if (cmp[0]->get_etat() == ETATZERO) { // no work needed ...
172  f_r.set_etat_zero() ;
173  resu.set_vr_mu(f_r, mu_resu) ;
174  return resu ;
175  }
176 
177  // Some working objects
178  //---------------------
179  const Mg3d& mg = *mpaff->get_mg() ;
180  int nz = mg.get_nzone() ; int nzm1 = nz - 1;
181  assert(mg.get_type_r(0) == RARE) ;
182  assert(mg.get_type_r(nzm1) == UNSURR) ;
183  Scalar S_r = *cmp[0] ;
184  Scalar S_eta = eta() ;
185  S_r.set_spectral_va().ylm() ;
186  S_eta.set_spectral_va().ylm() ;
187  const Base_val& base = S_eta.get_spectral_va().base ;
188  Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
189  Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
190  Mtbl_cf solution_hom_un(mg, base) ; solution_hom_un.annule_hard() ;
191  Mtbl_cf solution_hom_deux(mg, base) ; solution_hom_deux.annule_hard() ;
192 
193  // Build-up & inversion of the system for (eta, V^r) in each domain
194  //-----------------------------------------------------------------
195 
196  // Nucleus
197  //--------
198  int nr = mg.get_nr(0) ;
199  int nt = mg.get_nt(0) ;
200  int np = mg.get_np(0) ;
201  double alpha = mpaff->get_alpha()[0] ;
202  double beta = mpaff->get_beta()[0] ;
203  int l_q = 0 ; int m_q = 0; int base_r = 0 ;
204  int nr0 = nr - 1 ; //one degree of freedom less because of division by r^2
205 
206  // Loop on l and m
207  //----------------
208  for (int k=0 ; k<np+1 ; k++) {
209  for (int j=0 ; j<nt ; j++) {
210  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
211  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
212  int dege1 = (l_q <3 ? 0 : 1) ; //degeneracy of eq.1
213  int dege2 = 0 ; //degeneracy of eq.2
214  int nr_eq1 = nr0 - dege1 ; //Eq.1 is for h (eta/r^2)
215  int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
216  int nrtot = nr_eq1 + nr_eq2 ;
217  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
218  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
219  Diff_x2dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
220  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
221  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
222 
223  // Building the operator
224  //----------------------
225  for (int lin=0; lin<nr_eq1; lin++) { //eq.1
226  for (int col=dege1; col<nr0; col++)
227  oper.set(lin,col-dege1)
228  = md2(lin,col) + 6*mxd(lin,col) + 6*mid(lin,col) ;
229  for (int col=dege2; col<nr0; col++)
230  oper.set(lin,col-dege2+nr_eq1) = -mxd(lin,col) -2*mid(lin,col) ;
231  }
232  for (int lin=0; lin<nr0; lin++) { //eq.2
233  for (int col=dege1; col<nr0; col++)
234  oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
235  for (int col=dege2; col<nr0; col++)
236  oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) + 4*mid(lin,col) ;
237  }
238  oper.set_lu() ;
239 
240  // Filling the r.h.s
241  //------------------
242  for (int i=0; i<nr_eq1; i++) //eq.1
243  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
244  for (int i=0; i<nr0; i++) //eq.2
245  sec_membre.set(i+nr_eq1) = 0 ;
246 
247  // Inversion of the "big" operator
248  //--------------------------------
249  Tbl big_res = oper.inverse(sec_membre) ;
250  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
251  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
252 
253  // Putting coefficients of h and v to individual arrays
254  //-----------------------------------------------------
255  for (int i=0; i<dege1; i++)
256  res_eta.set(i) = 0 ;
257  for (int i=dege1; i<nr0; i++)
258  res_eta.set(i) = big_res(i-dege1) ;
259  res_eta.set(nr0) = 0 ;
260  for (int i=0; i<dege2; i++)
261  res_vr.set(i) = 0 ;
262  for (int i=dege2; i<nr0; i++)
263  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
264  res_vr.set(nr0) = 0 ;
265 
266  // Multiplication by xi^2 (the alpha^2 factor comes soon)
267  //-------------------------------------------------------
268  multx2_1d(nr, &res_eta.t, base_r) ;
269  multx2_1d(nr, &res_vr.t, base_r) ;
270 
271  // Homogeneous solution (only r^(l-1) in the nucleus)
272  Tbl sol_hom = solh(nr, l_q-1, 0., base_r) ;
273  for (int i=0 ; i<nr ; i++) {
274  sol_part_eta.set(0, k, j, i) = alpha*alpha*res_eta(i) ;
275  sol_part_vr.set(0, k, j, i) = alpha*alpha*res_vr(i) ;
276  solution_hom_un.set(0, k, j, i) = sol_hom(i) ;
277  solution_hom_deux.set(0, k, j, i) = 0. ;
278  }
279  }
280  }
281  }
282 
283 
284  // Shells
285  //-------
286  for (int zone=1 ; zone<nzm1 ; zone++) {
287  nr = mg.get_nr(zone) ;
288  assert (nr > 5) ;
289  assert(nt == mg.get_nt(zone)) ;
290  assert(np == mg.get_np(zone)) ;
291  alpha = mpaff->get_alpha()[zone] ;
292  beta = mpaff->get_beta()[zone] ;
293  double ech = beta / alpha ;
294 
295  // Loop on l and m
296  //----------------
297  for (int k=0 ; k<np+1 ; k++) {
298  for (int j=0 ; j<nt ; j++) {
299  base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
300  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
301  int dege1 = 2 ; //degeneracy of eq.1
302  int dege2 = 0 ; //degeneracy of eq.2
303  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
304  int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
305  int nrtot = nr_eq1 + nr_eq2 + 1;
306  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
307  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
308  Diff_x2dsdx2 x2d2(base_r, nr+1); const Matrice& m2d2 = x2d2.get_matrice() ;
309  Diff_xdsdx2 xd2(base_r, nr+1) ; const Matrice& mxd2 = xd2.get_matrice() ;
310  Diff_dsdx2 d2(base_r, nr+1) ; const Matrice& md2 = d2.get_matrice() ;
311  Diff_xdsdx xd(base_r, nr+1) ; const Matrice& mxd = xd.get_matrice() ;
312  Diff_dsdx d1(base_r, nr+1) ; const Matrice& md = d1.get_matrice() ;
313  Diff_id id(base_r, nr+1) ; const Matrice& mid = id.get_matrice() ;
314 
315  // Building the operator
316  //----------------------
317  for (int lin=0; lin<nr_eq1; lin++) {
318  for (int col=dege1; col<nr; col++)
319  oper.set(lin,col-dege1)
320  = mxd2(lin,col) + ech*md2(lin,col) + 2*md(lin,col) ;
321  for (int col=dege2; col<nr+1; col++)
322  oper.set(lin,col-dege2+nr_eq1) = -md(lin,col) ;
323  }
324  for (int lin=0; lin<nr_eq2; lin++) {
325  for (int col=dege1; col<nr; col++)
326  oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
327  for (int col=dege2; col<nr+1; col++)
328  oper.set(lin+nr_eq1, col-dege2+nr_eq1) =
329  mxd(lin,col) + ech*md(lin,col) + 2*mid(lin,col) ;
330  }
331  //Additional line to avoid spurious homogeneous solutions
332  //this line is the first one of the V^r eq.
333  for (int col=dege1; col<nr; col++)
334  oper.set(nrtot-1, col-dege1) = 0 ;
335  for (int col=dege2; col<nr+1; col++)
336  oper.set(nrtot-1, col-dege2+nr_eq1) =
337  m2d2(0,col) + ech*(2*mxd2(0,col) + ech*md2(0,col))
338  +4*(mxd(0,col) +ech*md(0,col))
339  +(2 - l_q*(l_q+1))*mid(0,col) ;
340  oper.set_lu() ;
341 
342  // Filling the r.h.s
343  //------------------
344  Tbl sr(5) ; sr.set_etat_qcq() ;
345  Tbl seta(nr) ; seta.set_etat_qcq() ;
346  for (int i=0; i<5; i++) {
347  sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
348  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
349  }
350  for (int i=5; i<nr; i++)
351  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
352  Tbl xsr= sr ; Tbl x2sr= sr ;
353  Tbl xseta= seta ;
354  multx2_1d(5, &x2sr.t, base_r) ; multx_1d(5, &xsr.t, base_r) ;
355  multx_1d(nr, &xseta.t, base_r) ;
356 
357  for (int i=0; i<nr_eq1; i++)
358  sec_membre.set(i) = alpha*(alpha*xseta(i) + beta*seta(i)) ;
359  for (int i=0; i<nr_eq2; i++)
360  sec_membre.set(i+nr_eq1) = 0 ;
361  sec_membre.set(nr_eq1+nr_eq2) = alpha*alpha*x2sr(0) + 2*alpha*beta*xsr(0)
362  + beta*beta*sr(0) ;
363 
364  // Inversion of the "big" operator
365  //--------------------------------
366  Tbl big_res = oper.inverse(sec_membre) ;
367  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
368  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
369 
370  // Putting coefficients of h and v to individual arrays
371  //-----------------------------------------------------
372  for (int i=0; i<dege1; i++)
373  res_eta.set(i) = 0 ;
374  for (int i=dege1; i<nr; i++)
375  res_eta.set(i) = big_res(i-dege1) ;
376  for (int i=0; i<dege2; i++)
377  res_vr.set(i) = 0 ;
378  for (int i=dege2; i<nr; i++)
379  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
380 
381  //homogeneous solutions
382  Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
383  Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
384  for (int i=0 ; i<nr ; i++) {
385  sol_part_eta.set(zone, k, j, i) = res_eta(i) ;
386  sol_part_vr.set(zone, k, j, i) = res_vr(i) ;
387  solution_hom_un.set(zone, k, j, i) = sol_hom1(0,i) ;
388  solution_hom_deux.set(zone, k, j, i) = sol_hom2(1,i) ;
389  }
390  }
391  }
392  }
393  }
394 
395  // Compactified external domain
396  //-----------------------------
397  nr = mg.get_nr(nzm1) ;
398  assert(nt == mg.get_nt(nzm1)) ;
399  assert(np == mg.get_np(nzm1)) ;
400  alpha = mpaff->get_alpha()[nzm1] ;
401  assert (nr > 4) ;
402  nr0 = nr - 2 ; //two degrees of freedom less because of division by r^2
403 
404  // Loop on l and m
405  //----------------
406  for (int k=0 ; k<np+1 ; k++) {
407  for (int j=0 ; j<nt ; j++) {
408  base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
409  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
410  int dege1 = 0; //degeneracy of eq.1
411  int dege2 = 1; //degeneracy of eq.2
412  int nr_eq1 = nr0 - dege1 ; //Eq.1 is for eta
413  int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
414  int nrtot = nr_eq1 + nr_eq2 ;
415  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
416  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
417  Diff_x2dsdx2 x2d2(base_r, nr) ; const Matrice& m2d2 = x2d2.get_matrice() ;
418  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
419  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
420 
421  // Building the operator
422  //----------------------
423  for (int lin=0; lin<nr_eq1; lin++) {
424  for (int col=dege1; col<nr0; col++)
425  oper.set(lin,col-dege1)
426  = m2d2(lin,col) + 4*mxd(lin,col) + 2*mid(lin,col) ;
427  for (int col=dege2; col<nr0; col++)
428  oper.set(lin,col-dege2+nr_eq1) =
429  mxd(lin,col) + 2*mid(lin,col) ;
430  }
431  for (int lin=0; lin<nr_eq2; lin++) {
432  for (int col=dege1; col<nr0; col++)
433  oper.set(lin+nr_eq1,col-dege1) = l_q*(l_q+1)*mid(lin,col) ;
434  for (int col=dege2; col<nr0; col++)
435  oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) ;
436  }
437  oper.set_lu() ;
438 
439  // Filling the r.h.s
440  //------------------
441  for (int i=0; i<nr_eq1; i++)
442  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
443  for (int i=0; i<nr_eq2; i++)
444  sec_membre.set(i+nr_eq1) = 0 ;
445  Tbl big_res = oper.inverse(sec_membre) ;
446  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
447  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
448 
449  // Putting coefficients of h and v to individual arrays
450  //-----------------------------------------------------
451  for (int i=0; i<dege1; i++)
452  res_eta.set(i) = 0 ;
453  for (int i=dege1; i<nr0; i++)
454  res_eta.set(i) = big_res(i-dege1) ;
455  res_eta.set(nr0) = 0 ;
456  res_eta.set(nr0+1) = 0 ;
457  for (int i=0; i<dege2; i++)
458  res_vr.set(i) = 0 ;
459  for (int i=dege2; i<nr0; i++)
460  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
461  res_vr.set(nr0) = 0 ;
462  res_vr.set(nr0+1) = 0 ;
463 
464  // Multiplication by r^2
465  //-----------------------
466  Tbl res_v2(nr) ; res_v2.set_etat_qcq() ;
467  Tbl res_e2(nr) ; res_e2.set_etat_qcq() ;
468  mult2_xm1_1d_cheb(nr, res_eta.t, res_e2.t) ;
469  mult2_xm1_1d_cheb(nr, res_vr.t, res_v2.t) ;
470 
471  // Homogeneous solution (only 1/r^(l+2) in the CED)
472  Tbl sol_hom = solh(nr, l_q+1, 0., base_r) ;
473  for (int i=0 ; i<nr ; i++) {
474  sol_part_eta.set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
475  sol_part_vr.set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
476  solution_hom_un.set(nzm1, k, j, i) = sol_hom(i) ;
477  solution_hom_deux.set(nzm1, k, j, i) = 0. ;
478  }
479  }
480  }
481  }
482 
483  // Now let's match everything ...
484  //-------------------------------
485 
486  // Resulting V^r & eta
487  Scalar vr(*mpaff) ; vr.set_etat_qcq() ;
488  vr.set_spectral_base(base) ;
490  Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
491  cf_vr.annule_hard() ;
492  Scalar het(*mpaff) ; het.set_etat_qcq() ;
493  het.set_spectral_base(base) ;
495  Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
496  cf_eta.annule_hard() ;
497  int taille = 2*nzm1 ;
498  Tbl sec_membre(taille) ;
499  Matrice systeme(taille, taille) ;
500  systeme.set_etat_qcq() ;
501  int ligne ; int colonne ;
502 
503  // Loop on l and m
504  //----------------
505  for (int k=0 ; k<np+1 ; k++)
506  for (int j=0 ; j<nt ; j++) {
507  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
508  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
509 
510  ligne = 0 ;
511  colonne = 0 ;
512  sec_membre.annule_hard() ;
513  for (int l=0; l<taille; l++)
514  for (int c=0; c<taille; c++)
515  systeme.set(l,c) = 0 ;
516  //Nucleus
517  nr = mg.get_nr(0) ;
518  alpha = mpaff->get_alpha()[0] ;
519  // value of x^(l-1) at 1 ...
520  systeme.set(ligne, colonne) = 1. ;
521  for (int i=0 ; i<nr ; i++)
522  sec_membre.set(ligne) -= sol_part_eta(0, k, j, i) ;
523  ligne++ ;
524  // ... and of its couterpart for V^r
525  systeme.set(ligne, colonne) = l_q;
526  for (int i=0; i<nr; i++)
527  sec_membre.set(ligne) -= sol_part_vr(0,k,j,i) ;
528  colonne++ ;
529  //shells
530  for (int zone=1 ; zone<nzm1 ; zone++) {
531  nr = mg.get_nr(zone) ;
532  alpha = mpaff->get_alpha()[zone] ;
533  double echelle = mpaff->get_beta()[zone]/alpha ;
534  ligne -- ;
535  //value of (x+echelle)^(l-1) at -1
536  systeme.set(ligne, colonne) = -pow(echelle-1., double(l_q-1)) ;
537  // value of 1/(x+echelle) ^(l+2) at -1
538  systeme.set(ligne, colonne+1) = -1/pow(echelle-1., double(l_q+2)) ;
539  for (int i=0 ; i<nr ; i++)
540  if (i%2 == 0)
541  sec_membre.set(ligne) += sol_part_eta(zone, k, j, i) ;
542  else sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
543  ligne++ ;
544  // ... and their couterparts for V^r
545  systeme.set(ligne, colonne) = -l_q*pow(echelle-1., double(l_q-1)) ;
546  systeme.set(ligne, colonne+1) = (l_q+1)/pow(echelle-1., double(l_q+2));
547  for (int i=0 ; i<nr ; i++)
548  if (i%2 == 0)
549  sec_membre.set(ligne) += sol_part_vr(zone, k, j, i) ;
550  else sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
551  ligne++ ;
552 
553  //value of (x+echelle)^(l-1) at 1 :
554  systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
555  // value of 1/(x+echelle)^(l+2) at 1 :
556  systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
557  for (int i=0 ; i<nr ; i++)
558  sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
559  ligne ++ ;
560  //... and their couterparts for V^r
561  systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
562  systeme.set(ligne, colonne+1) = -double(l_q+1)
563  / pow(echelle+1., double(l_q+2)) ;
564  for (int i=0 ; i<nr ; i++)
565  sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i);
566  colonne += 2 ;
567  }
568  //Compactified external domain
569  nr = mg.get_nr(nzm1) ;
570 
571  alpha = mpaff->get_alpha()[nzm1] ;
572  ligne -- ;
573  //value of (x-1)^(l+2) at -1 :
574  systeme.set(ligne, colonne) = -pow(-2, double(l_q+2)) ;
575  for (int i=0 ; i<nr ; i++)
576  if (i%2 == 0) sec_membre.set(ligne) += sol_part_eta(nzm1, k, j, i) ;
577  else sec_membre.set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
578  //... and of its couterpart for V^r
579  systeme.set(ligne+1, colonne) = double(l_q+1)*pow(-2, double(l_q+2)) ;
580  for (int i=0 ; i<nr ; i++)
581  if (i%2 == 0) sec_membre.set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
582  else sec_membre.set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
583 
584  // Solution of the system giving the coefficients for the homogeneous
585  // solutions
586  //-------------------------------------------------------------------
587  if (taille > 2) systeme.set_band(2,2) ;
588  else systeme.set_band(1,1) ;
589  systeme.set_lu() ;
590  Tbl facteurs(systeme.inverse(sec_membre)) ;
591  int conte = 0 ;
592 
593  // everything is put to the right place, the same combination of hom.
594  // solutions (with some l or -(l+1) factors) must be used for V^r
595  //-------------------------------------------------------------------
596  nr = mg.get_nr(0) ; //nucleus
597  for (int i=0 ; i<nr ; i++) {
598  cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
599  +facteurs(conte)*solution_hom_un(0, k, j, i) ;
600  cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
601  +double(l_q)*facteurs(conte)*solution_hom_un(0, k, j, i) ;
602  }
603  conte++ ;
604  for (int zone=1 ; zone<nzm1 ; zone++) { //shells
605  nr = mg.get_nr(zone) ;
606  for (int i=0 ; i<nr ; i++) {
607  cf_eta.set(zone, k, j, i) =
608  sol_part_eta(zone, k, j, i)
609  +facteurs(conte)*solution_hom_un(zone, k, j, i)
610  +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
611  cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
612  +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
613  -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
614  }
615  conte+=2 ;
616  }
617  nr = mg.get_nr(nz-1) ; //compactified external domain
618  for (int i=0 ; i<nr ; i++) {
619  cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
620  +facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
621  cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
622  -double(l_q+1)*facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
623  }
624  } // End of nullite_plm
625  } //End of loop on theta
626 
627  vr.set_spectral_va().ylm_i() ;
628  het.set_spectral_va().ylm_i() ;
629 
630  resu.set_vr_eta_mu(vr, het, mu_resu) ;
631 
632  return resu ;
633 
634 }
635 
636 }
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_x2dsdx2.C:98
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:172
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:604
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
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 set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition: valeur.C:715
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:395
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
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
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:777
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:97
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
const Metric *const met_div
Metric with respect to which the divergence is defined.
Definition: vector.h:730
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1052
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
void div_r()
Division by r everywhere; dzpuis is not changed.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx2.C:100
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:490
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
Matrix handling.
Definition: matrice.h:152
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx2.C:94
double * t
The array of double.
Definition: tbl.h:176
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:608
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:433
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source: .
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:367
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:803
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 annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:315
Affine radial mapping.
Definition: map.h:2042
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
void set_vr_mu(const Scalar &vr_i, const Scalar &mu_i)
Sets the angular potentials (see member p_mu ), and the component of the vector.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:107
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:531
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
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
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1007
Divergence-free vectors.
Definition: vector.h:724
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:101
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through , and .
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607