LORENE: tensorelliptic.C Source File

00001 #include<math.h>
00002 // Lorene headers
00003 #include "metric.h"
00004 #include "nbr_spx.h"
00005 #include "utilitaires.h"
00006 #include "graphique.h"
00007 #include "proto.h"
00008 #include "diff.h"
00009 
00010 
00011  void tensorelliptic( Scalar source, Scalar& resu, double fit, double fit2, double fit0d2, double fit1d2, double fit0d3, double fit1d3) {
00012 
00013   const int nz = (*source.get_mp().get_mg()).get_nzone();   // Number of domains
00014   int nr = (*source.get_mp().get_mg()).get_nr(1);   // Number of collocation points in r in each domain
00015   int nt = (*source.get_mp().get_mg()).get_nt(1);   // Number of collocation points in theta in each domain
00016   int np = (*source.get_mp().get_mg()).get_np(1);   // Number of collocation points in phi in each domain
00017   const Map_af* map = dynamic_cast<const Map_af*>(&source.get_mp()) ;
00018   const Mg3d* mgrid = (*map).get_mg();
00019 
00020     // Some helpful stuff...
00021    const Coord& rr = (*map).r;
00022    Scalar rrr (*map) ; 
00023    rrr = rr ; 
00024    rrr.set_spectral_va().set_base(source.get_spectral_va().base); 
00025    
00026     Scalar source_coq = source ;
00027     source_coq.mult_r() ;
00028     source_coq.mult_r() ;
00029     source.set_spectral_va().ylm() ;
00030     source_coq.set_spectral_va().ylm() ;
00031     Scalar phi(source.get_mp()) ;
00032     phi.annule_hard() ;
00033     //  phi.std_spectral_base();
00034         phi.set_spectral_va().set_base(source.get_spectral_va().base) ;
00035     phi.set_spectral_va().ylm() ;   
00036     Mtbl_cf& sol_coef = (*phi.set_spectral_va().c_cf) ;
00037 
00038     const Base_val& base = source.get_spectral_base() ;
00039     Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
00040     Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
00041     Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
00042 
00043     int l_q, m_q, base_r ;
00044 
00045     { int lz = 0 ;
00046         for (int k=0; k < np; k++)
00047         for (int j=0; j<nt; j++) {
00048             base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00049             if (nullite_plm(j, nt, k, np, base) == 1) {
00050               for (int ii=0 ; ii<nr ; ii++){
00051                 sol_hom1.set(lz, k, j, ii) = 0 ;
00052                 sol_part.set(lz, k, j, ii) = 0 ;
00053               }
00054 
00055             }
00056         }
00057     }
00058 
00059 
00060     { 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.
00061         double alpha = (*map).get_alpha()[lz] ;
00062         double beta = (*map).get_beta()[lz] ;
00063         double ech = beta / alpha ; 
00064         for (int k=0; k < np; k++)
00065           for (int j=0; j<nt; j++) {
00066         base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00067         if (nullite_plm(j, nt, k, np, base) == 1) {
00068           
00069           
00070           
00071           Matrice ope(nr,nr) ;
00072           ope.annule_hard() ;
00073           
00074           Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
00075           Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00076           Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
00077           Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
00078           Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
00079           Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
00080           Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
00081           Diff_x4dsdx2 x4dx2 (base_r, nr); const Matrice& mx4dx2 = x4dx2.get_matrice();
00082           
00083 
00084           ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
00085             - 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) 
00086                       + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2) 
00087                       + 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));
00088           
00089 
00090 
00091 
00092           for (int col=0; col<nr; col++)
00093             ope.set(nr-1, col) = 0 ;
00094           ope.set(nr-1, 0) = 1 ;
00095 
00096           
00097           Tbl rhs(nr); 
00098           rhs.annule_hard() ;
00099           for (int i=0; i<nr; i++)
00100             rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(1, k, j, i) ;
00101           rhs.set(nr-1) = 0 ;
00102           ope.set_lu() ;
00103           Tbl sol = ope.inverse(rhs) ;
00104 
00105           
00106           for (int i=0; i<nr; i++)
00107             sol_part.set(1, k, j, i) = sol(i) ;
00108           
00109           rhs.annule_hard();
00110           rhs.set(nr-1) = 1 ;
00111           sol = ope.inverse(rhs) ;
00112 
00113 
00114           for (int i=0; i<nr; i++)
00115             sol_hom1.set(1, k, j, i) = sol(i) ;              
00116           
00117           
00118         }
00119           }
00120     }
00121     
00122     
00123     // Current implementations only allow grids with more than 3 shells. Zones 2 and 3 are also treated separately, allowing
00124     // A smoother interpolation of the real operator;
00125     
00126           { int lz = 2 ;
00127 
00128         double alpha = (*map).get_alpha()[lz] ;
00129         double beta = (*map).get_beta()[lz] ;
00130         double ech = beta / alpha ; 
00131     
00132       for (int k=0; k < np; k++)
00133         for (int j=0; j<nt; j++) {
00134           base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00135           if (nullite_plm(j, nt, k, np, base) == 1) {
00136         
00137         Matrice ope(nr,nr) ;
00138         ope.annule_hard() ;
00139         
00140         Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
00141         Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00142         Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
00143         Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
00144         Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
00145         Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
00146         Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
00147     
00148             ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
00149         - 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)  ;
00150 
00151 
00152 
00153         for (int col=0; col<nr; col++)
00154           ope.set(nr-1, col) = 0 ;
00155         ope.set(nr-1, 0) = 1 ;
00156         for (int col=0; col<nr; col++) {
00157           ope.set(nr-2, col) = 0 ;
00158         }
00159         ope.set(nr-2, 1) = 1 ;
00160     
00161         Tbl rhs(nr) ;
00162         rhs.annule_hard() ;
00163         for (int i=0; i<nr; i++)
00164           rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
00165         rhs.set(nr-2) = 0 ;
00166         rhs.set(nr-1) = 0 ;
00167         ope.set_lu() ;
00168         Tbl sol = ope.inverse(rhs) ;
00169 
00170         for (int i=0; i<nr; i++)
00171                 sol_part.set(lz, k, j, i) = sol(i) ;
00172         
00173         rhs.annule_hard() ;
00174         rhs.set(nr-2) = 1 ;
00175         sol = ope.inverse(rhs) ;
00176         for (int i=0; i<nr; i++)
00177           sol_hom1.set(lz, k, j, i) = sol(i) ;
00178         
00179         rhs.set(nr-2) = 0 ;
00180         rhs.set(nr-1) = 1 ;
00181         sol = ope.inverse(rhs) ;
00182         for (int i=0; i<nr; i++)
00183           sol_hom2.set(lz, k, j, i) = sol(i) ;
00184         
00185           }
00186         }
00187       
00188       } 
00189 
00190 
00191     
00192           { int lz = 3 ;
00193 
00194         double alpha = (*map).get_alpha()[lz] ;
00195         double beta = (*map).get_beta()[lz] ;
00196         double ech = beta / alpha ; 
00197     
00198       for (int k=0; k < np; k++)
00199         for (int j=0; j<nt; j++) {
00200           base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00201           if (nullite_plm(j, nt, k, np, base) == 1) {
00202         
00203         Matrice ope(nr,nr) ;
00204         ope.annule_hard() ;
00205         
00206         Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
00207         Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00208         Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
00209         Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
00210         Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
00211         Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
00212         Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
00213     
00214             ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
00215         - 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)  ;
00216 
00217         for (int col=0; col<nr; col++)
00218           ope.set(nr-1, col) = 0 ;
00219         ope.set(nr-1, 0) = 1 ;
00220         for (int col=0; col<nr; col++) {
00221           ope.set(nr-2, col) = 0 ;
00222         }
00223         ope.set(nr-2, 1) = 1 ;
00224     
00225         Tbl rhs(nr) ;
00226         rhs.annule_hard() ;
00227         for (int i=0; i<nr; i++)
00228           rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
00229         rhs.set(nr-2) = 0 ;
00230         rhs.set(nr-1) = 0 ;
00231         ope.set_lu() ;
00232         Tbl sol = ope.inverse(rhs) ;
00233 
00234         for (int i=0; i<nr; i++)
00235                 sol_part.set(lz, k, j, i) = sol(i) ;
00236         
00237         rhs.annule_hard() ;
00238         rhs.set(nr-2) = 1 ;
00239         sol = ope.inverse(rhs) ;
00240         for (int i=0; i<nr; i++)
00241           sol_hom1.set(lz, k, j, i) = sol(i) ;
00242         
00243         rhs.set(nr-2) = 0 ;
00244         rhs.set(nr-1) = 1 ;
00245         sol = ope.inverse(rhs) ;
00246         for (int i=0; i<nr; i++)
00247           sol_hom2.set(lz, k, j, i) = sol(i) ;
00248         
00249           }
00250         }
00251       
00252       } 
00253 
00254     
00255        
00256     
00257     for (int lz=4; lz<nz-1; lz++) {
00258       double ech = (*map).get_beta()[lz] / (*map).get_alpha()[lz] ;
00259       for (int k=0; k < np; k++)
00260         for (int j=0; j<nt; j++) {
00261           base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00262           if (nullite_plm(j, nt, k, np, base) == 1) {
00263         
00264         Matrice ope(nr,nr) ;
00265         ope.annule_hard() ;
00266         
00267         Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
00268         Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00269         Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
00270         Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
00271         Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
00272         Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
00273         ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
00274           - l_q*(l_q+1)*mid ;
00275         
00276         for (int col=0; col<nr; col++)
00277           ope.set(nr-1, col) = 0 ;
00278         ope.set(nr-1, 0) = 1 ;
00279         for (int col=0; col<nr; col++) {
00280           ope.set(nr-2, col) = 0 ;
00281         }
00282         ope.set(nr-2, 1) = 1 ;
00283     
00284         Tbl rhs(nr) ;
00285         rhs.annule_hard() ;
00286         for (int i=0; i<nr; i++)
00287           rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
00288         rhs.set(nr-2) = 0 ;
00289         rhs.set(nr-1) = 0 ;
00290         ope.set_lu() ;
00291         Tbl sol = ope.inverse(rhs) ;
00292 
00293         for (int i=0; i<nr; i++)
00294                 sol_part.set(lz, k, j, i) = sol(i) ;
00295         
00296         rhs.annule_hard() ;
00297         rhs.set(nr-2) = 1 ;
00298         sol = ope.inverse(rhs) ;
00299         for (int i=0; i<nr; i++)
00300           sol_hom1.set(lz, k, j, i) = sol(i) ;
00301         
00302         rhs.set(nr-2) = 0 ;
00303         rhs.set(nr-1) = 1 ;
00304         sol = ope.inverse(rhs) ;
00305         for (int i=0; i<nr; i++)
00306           sol_hom2.set(lz, k, j, i) = sol(i) ;
00307         
00308           }
00309         }
00310       
00311     }
00312     { int lz = nz-1 ;
00313     double alpha = (*map).get_alpha()[lz] ;
00314     for (int k=0; k < np; k++)
00315       for (int j=0; j<nt; j++) {
00316         base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
00317         if (nullite_plm(j, nt, k, np, base) == 1) {
00318           
00319           Matrice ope(nr,nr) ;
00320           ope.annule_hard() ;
00321           Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00322           Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
00323           
00324             ope = (mdx2 - l_q*(l_q+1)*ms2)/(alpha*alpha) ;
00325             
00326             for (int i=0; i<nr; i++)
00327               ope.set(nr-1, i) = 0 ;
00328             ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
00329             
00330             for (int i=0; i<nr; i++) {
00331               ope.set(nr-2, i) = 1 ; //for the limit at infinity
00332             }
00333 
00334             if (l_q > 0) {
00335                 for (int i=0; i<nr; i++) {
00336                 ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
00337                 }
00338             }
00339         
00340             Tbl rhs(nr) ;
00341             rhs.annule_hard() ;
00342             for (int i=0; i<nr; i++)
00343                 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, k, j, i) ;
00344             if (l_q>0) rhs.set(nr-3) = 0 ;
00345             rhs.set(nr-2) = 0 ;
00346             rhs.set(nr-1) = 0 ;
00347             ope.set_lu() ;
00348             Tbl sol = ope.inverse(rhs) ;
00349                
00350             for (int i=0; i<nr; i++)
00351                 sol_part.set(lz, k, j, i) = sol(i) ;
00352 
00353             rhs.annule_hard() ;
00354             rhs.set(nr-1) = 1 ;
00355             sol = ope.inverse(rhs) ;
00356             for (int i=0; i<nr; i++)
00357                 sol_hom2.set(lz, k, j, i) = sol(i) ;
00358 
00359             }
00360         }
00361     }
00362 
00363     Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
00364     Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
00365     Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
00366 
00367 
00368     // Now matching the homogeneous solutions between the different domains...
00369     
00370     for (int k=0; k < np; k++)
00371         for (int j=0; j<nt; j++) {
00372         base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
00373         if (nullite_plm(j, nt, k, np, base) == 1) {
00374             Matrice systeme(2*nz-4, 2*nz-4) ;
00375             systeme.annule_hard() ;
00376             Tbl rhs(2*nz-4) ;
00377             rhs.annule_hard() ;
00378 
00379             //First shell 
00380             int lin = 0 ;
00381             int col = 0 ;
00382             
00383             double alpha = (*map).get_alpha()[1] ;
00384             
00385             systeme.set(lin, col) += sol_hom1.val_out_bound_jk(1, j, k) ;
00386             rhs.set(lin) -= sol_part.val_out_bound_jk(1, j, k) ;
00387             
00388             lin++ ;
00389             systeme.set(lin, col) += dhom1.val_out_bound_jk(1, j, k) / alpha ;
00390             rhs.set(lin) -=  dpart.val_out_bound_jk(1, j, k) / alpha ;
00391             col += 1 ;
00392             
00393             //Shells
00394             for (int lz=2; lz<nz-1; lz++) {
00395               alpha = (*map).get_alpha()[lz] ;
00396               lin-- ;
00397               systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, j, k) ;
00398               systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, j, k) ;
00399               rhs.set(lin) += sol_part.val_in_bound_jk(lz, j, k) ;
00400               
00401               lin++ ;
00402               systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, j, k) / alpha ;
00403               systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, j, k) / alpha ;
00404               rhs.set(lin) += dpart.val_in_bound_jk(lz, j, k) / alpha;
00405               
00406               lin++ ;
00407               systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, j, k) ;
00408               systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, j, k) ;
00409               rhs.set(lin) -= sol_part.val_out_bound_jk(lz, j, k) ;
00410               
00411               lin++ ;
00412               systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, j, k) / alpha ;
00413               systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, j, k) / alpha ;
00414               rhs.set(lin) -= dpart.val_out_bound_jk(lz, j, k) / alpha ;
00415               col += 2 ;
00416             }
00417             
00418             //CED
00419             alpha = (*map).get_alpha()[nz-1] ;
00420             lin-- ;
00421             systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, j, k) ;
00422             rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, j, k) ;
00423 
00424             lin++ ;
00425             systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, j, k) ;
00426             rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, j, k) ;
00427 
00428             systeme.set_lu() ;
00429 
00430             Tbl coef = systeme.inverse(rhs);
00431             int indice = 0 ;
00432                         
00433 
00434             for (int i=0; i<(*mgrid).get_nr(1); i++)
00435               sol_coef.set(1, k, j, i) = sol_part(1, k, j, i)
00436             +coef(indice)*sol_hom1(1, k, j, i) ;
00437             indice +=1;
00438             
00439             
00440             for (int lz=2; lz<nz-1; lz++) {
00441               for (int i=0; i<(*mgrid).get_nr(lz); i++)
00442             sol_coef.set(lz, k, j, i) = sol_part(lz, k, j, i)
00443               +coef(indice)*sol_hom1(lz, k, j, i) 
00444               +coef(indice+1)*sol_hom2(lz, k, j, i) ;
00445               indice += 2 ;
00446             }
00447             for (int i=0; i<(*mgrid).get_nr(nz-1); i++)
00448               sol_coef.set(nz-1, k, j, i) = sol_part(nz-1, k, j, i)
00449             +coef(indice)*sol_hom2(nz-1, k, j, i) ;
00450 
00451         }
00452         }
00453     
00454     
00455     delete phi.set_spectral_va().c ;
00456     phi.set_spectral_va().c = 0x0 ;
00457     phi.set_spectral_va().ylm_i() ;
00458     
00459     phi.annule_domain(nz-1);
00460     
00461     resu = phi;
00462 
00463 }