LORENE: pois_vect_r0.C Source File

00001 /*
00002  *  Solution of the l=0 part of the vector Poisson equation (only r-component)
00003  *
00004  */
00005 
00006 /*
00007  *   Copyright (c) 2007  Jerome Novak
00008  *
00009  *   This file is part of LORENE.
00010  *
00011  *   LORENE is free software; you can redistribute it and/or modify
00012  *   it under the terms of the GNU General Public License version 2
00013  *   as published by the Free Software Foundation.
00014  *
00015  *   LORENE is distributed in the hope that it will be useful,
00016  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00017  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018  *   GNU General Public License for more details.
00019  *
00020  *   You should have received a copy of the GNU General Public License
00021  *   along with LORENE; if not, write to the Free Software
00022  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023  *
00024  */
00025 
00026 char pois_vect_r0_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pois_vect_r0.C,v 1.1 2007/01/23 17:08:46 j_novak Exp $" ;
00027 
00028 /*
00029  * $Id: pois_vect_r0.C,v 1.1 2007/01/23 17:08:46 j_novak Exp $
00030  * $Log: pois_vect_r0.C,v $
00031  * Revision 1.1  2007/01/23 17:08:46  j_novak
00032  * New function pois_vect_r0.C to solve the l=0 part of the vector Poisson
00033  * equation, which involves only the r-component.
00034  *
00035  *
00036  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pois_vect_r0.C,v 1.1 2007/01/23 17:08:46 j_novak Exp $
00037  *
00038  */
00039 
00040 // Lorene headers
00041 #include "metric.h"
00042 #include "proto.h"
00043 #include "diff.h"
00044 
00045 /*
00046  * This function solves for the l=0 component of
00047  * 
00048  *         d2 f   2 df   2f
00049  *         ---- + - -- - --  = source 
00050  *          dr2   r dr   r2
00051  * 
00052  * and returns the soluton f. 
00053  * The input Scalar must have dzpuis = 4.
00054  */
00055 Scalar pois_vect_r0(const Scalar& source) {
00056 
00057     const Map& map0 = source.get_mp() ;
00058     const Map_af* map1 = dynamic_cast<const Map_af*>(&map0) ;
00059     assert(map1 != 0x0) ;
00060     const Map_af& map = *map1 ;
00061 
00062     const Mg3d& mgrid = *map.get_mg() ;
00063     int nz = mgrid.get_nzone() ;
00064 
00065     Scalar resu(map) ;
00066     if (source.get_etat() == ETATZERO) {
00067     resu = 0 ;
00068     return resu ;
00069     }
00070 
00071     resu.annule_hard() ;
00072     resu.std_spectral_base_odd() ;
00073     resu.set_spectral_va().ylm() ;  
00074     Mtbl_cf& sol_coef = (*resu.set_spectral_va().c_cf) ;
00075     const Base_val& base = source.get_spectral_base() ;
00076     assert(resu.get_spectral_base() == base) ;
00077     assert(source.check_dzpuis(4)) ;
00078 
00079     Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
00080     Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
00081     Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
00082 
00083     { int lz = 0 ;
00084       int nr = mgrid.get_nr(lz) ;
00085       double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
00086       assert(mgrid.get_type_r(lz) == RARE) ;
00087       int base_r = R_CHEBI ;
00088       Matrice ope(nr,nr) ;
00089       ope.annule_hard() ;
00090       Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00091       Diff_sxdsdx sdx(base_r, nr) ; const Matrice& msdx = sdx.get_matrice() ;
00092       Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
00093       
00094       for (int lin=0; lin<nr-1; lin++) 
00095       for (int col=0; col<nr; col++)
00096           ope.set(lin,col) = (mdx2(lin,col) + 2*msdx(lin,col) - 2*ms2(lin,col))/alpha2 ;
00097       
00098       ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
00099       for (int i=1; i<nr; i++)
00100       ope.set(nr-1, i) = 0 ;
00101 
00102       Tbl rhs(nr) ;
00103       rhs.annule_hard() ;
00104       for (int i=0; i<nr-1; i++)
00105       rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
00106       rhs.set(nr-1) = 0 ;
00107       ope.set_lu() ;
00108       Tbl sol = ope.inverse(rhs) ;
00109                
00110       for (int i=0; i<nr; i++)
00111       sol_part.set(lz, 0, 0, i) = sol(i) ;
00112       
00113       rhs.annule_hard() ;
00114       rhs.set(nr-1) = 1 ;
00115       sol = ope.inverse(rhs) ;
00116       for (int i=0; i<nr; i++)
00117       sol_hom1.set(lz, 0, 0, i) = sol(i) ;
00118       
00119     }
00120 
00121     for (int lz=1; lz<nz-1; lz++) {
00122     int nr = mgrid.get_nr(lz) ;
00123     double alpha = map.get_alpha()[lz] ;
00124     double beta = map.get_beta()[lz] ;
00125     double ech = beta / alpha ;
00126     assert(mgrid.get_type_r(lz) == FIN) ;
00127     int base_r = R_CHEB ;
00128             
00129     Matrice ope(nr,nr) ;
00130     ope.annule_hard() ;
00131     
00132     Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
00133     Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00134     Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
00135     Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
00136     Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
00137     Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
00138     
00139     for (int lin=0; lin<nr-2; lin++) 
00140         for (int col=0; col<nr; col++) 
00141         ope.set(lin, col) = mx2dx2(lin, col) + 2*ech*mxdx2(lin, col) + ech*ech*mdx2(lin, col) 
00142             + 2*(mxdx(lin, col) + ech*mdx(lin, col)) - 2*mid(lin, col) ;
00143     
00144     ope.set(nr-2, 0) = 0 ;
00145     ope.set(nr-2, 1) = 1 ;
00146     for (int col=2; col<nr; col++) {
00147         ope.set(nr-2, col) = 0 ;
00148     }
00149     ope.set(nr-1, 0) = 1 ;
00150     for (int col=1; col<nr; col++)
00151         ope.set(nr-1, col) = 0 ;
00152     
00153     Tbl src(nr) ;
00154     src.set_etat_qcq() ;
00155     for (int i=0; i<nr; i++) 
00156         src.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
00157     Tbl xsrc = src ; multx_1d(nr, &xsrc.t, base_r) ;
00158     Tbl x2src = src ; multx2_1d(nr, &x2src.t, base_r) ;
00159     Tbl rhs(nr) ;
00160     rhs.set_etat_qcq() ;
00161     for (int i=0; i<nr-2; i++) 
00162         rhs.set(i) = beta*beta*src(i) + 2*beta*alpha*xsrc(i) + alpha*alpha*x2src(i) ;
00163     rhs.set(nr-2) = 0 ;
00164     rhs.set(nr-1) = 0 ;
00165     ope.set_lu() ;
00166     Tbl sol = ope.inverse(rhs) ;
00167     
00168     for (int i=0; i<nr; i++)
00169         sol_part.set(lz, 0, 0, i) = sol(i) ;
00170     
00171     rhs.annule_hard() ;
00172     rhs.set(nr-2) = 1 ;
00173     sol = ope.inverse(rhs) ;
00174     for (int i=0; i<nr; i++)
00175         sol_hom1.set(lz, 0, 0, i) = sol(i) ;
00176     
00177     rhs.set(nr-2) = 0 ;
00178     rhs.set(nr-1) = 1 ;
00179     sol = ope.inverse(rhs) ;
00180     for (int i=0; i<nr; i++)
00181         sol_hom2.set(lz, 0, 0, i) = sol(i) ;
00182     
00183     }
00184 
00185     { int lz = nz-1 ;
00186       int nr = mgrid.get_nr(lz) ;
00187       double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
00188       assert(mgrid.get_type_r(lz) == UNSURR) ;
00189       int base_r = R_CHEBU ;
00190             
00191       Matrice ope(nr,nr) ;
00192       ope.annule_hard() ;
00193       Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
00194       Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
00195       
00196       for (int lin=0; lin<nr-3; lin++) 
00197       for (int col=0; col<nr; col++) 
00198           ope.set(lin, col) = (mdx2(lin, col) - 2*ms2(lin, col))/alpha2 ;
00199       
00200       for (int i=0; i<nr; i++) {
00201       ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
00202       }
00203       
00204       for (int i=0; i<nr; i++) {
00205       ope.set(nr-2, i) = 1 ; //for the limit at inifinity
00206       }
00207       
00208       ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
00209       for (int i=1; i<nr; i++)
00210       ope.set(nr-1, i) = 0 ;
00211       
00212       Tbl rhs(nr) ;
00213       rhs.annule_hard() ;
00214       for (int i=0; i<nr-3; i++)
00215       rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
00216       rhs.set(nr-3) = 0 ;
00217       rhs.set(nr-2) = 0 ;
00218       rhs.set(nr-1) = 0 ;
00219       ope.set_lu() ;
00220       Tbl sol = ope.inverse(rhs) ;
00221       for (int i=0; i<nr; i++)
00222       sol_part.set(lz, 0, 0, i) = sol(i) ;
00223       
00224       rhs.annule_hard() ;
00225       rhs.set(nr-1) = 1 ;
00226       sol = ope.inverse(rhs) ;
00227       for (int i=0; i<nr; i++)
00228       sol_hom2.set(lz, 0, 0, i) = sol(i) ;
00229       
00230     }
00231 
00232     Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
00233     Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
00234     Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
00235     
00236     Matrice systeme(2*(nz-1), 2*(nz-1)) ;
00237     systeme.annule_hard() ;
00238     Tbl rhs(2*(nz-1)) ;
00239     rhs.annule_hard() ;
00240 
00241     //Nucleus
00242     int lin = 0 ;
00243     int col = 0 ;
00244     double alpha = map.get_alpha()[0] ;
00245     systeme.set(lin,col) = sol_hom1.val_out_bound_jk(0, 0, 0) ;
00246     rhs.set(lin) -= sol_part.val_out_bound_jk(0, 0, 0) ;
00247     lin++ ;
00248     systeme.set(lin,col) = dhom1.val_out_bound_jk(0, 0, 0) / alpha ;
00249     rhs.set(lin) -= dpart.val_out_bound_jk(0, 0, 0) / alpha ;
00250     col++ ;
00251     
00252     //Shells
00253     for (int lz=1; lz<nz-1; lz++) {
00254     alpha = map.get_alpha()[lz] ;
00255     lin-- ;
00256     systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, 0, 0) ;
00257     systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, 0, 0) ;
00258     rhs.set(lin) += sol_part.val_in_bound_jk(lz, 0, 0) ;
00259     
00260     lin++ ;
00261     systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, 0, 0) / alpha ;
00262     systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, 0, 0) / alpha ;
00263     rhs.set(lin) += dpart.val_in_bound_jk(lz, 0, 0) / alpha;
00264     
00265     lin++ ;
00266     systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, 0, 0) ;
00267     systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, 0, 0) ;
00268     rhs.set(lin) -= sol_part.val_out_bound_jk(lz, 0, 0) ;
00269     
00270     lin++ ;
00271     systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, 0, 0) / alpha ;
00272     systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, 0, 0) / alpha ;
00273     rhs.set(lin) -= dpart.val_out_bound_jk(lz, 0, 0) / alpha ;
00274     col += 2 ;
00275     }
00276     
00277     //CED
00278     alpha = map.get_alpha()[nz-1] ;
00279     lin-- ;
00280     systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, 0, 0) ;
00281     rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, 0, 0) ;
00282     
00283     lin++ ;
00284     systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, 0, 0) ;
00285     rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, 0, 0) ;
00286     
00287     systeme.set_lu() ;
00288     Tbl coef = systeme.inverse(rhs) ;
00289     int indice = 0 ;
00290     
00291     for (int i=0; i<mgrid.get_nr(0); i++)
00292     sol_coef.set(0, 0, 0, i) = sol_part(0, 0, 0, i) 
00293         + coef(indice)*sol_hom1(0, 0, 0, i) ;
00294     indice++ ;
00295     for (int lz=1; lz<nz-1; lz++) {
00296     for (int i=0; i<mgrid.get_nr(lz); i++)
00297         sol_coef.set(lz, 0, 0, i) = sol_part(lz, 0, 0, i)
00298         +coef(indice)*sol_hom1(lz, 0, 0, i) 
00299         +coef(indice+1)*sol_hom2(lz, 0, 0, i) ;
00300     indice += 2 ;
00301     }
00302     for (int i=0; i<mgrid.get_nr(nz-1); i++)
00303     sol_coef.set(nz-1, 0, 0, i) = sol_part(nz-1, 0, 0, i)
00304         +coef(indice)*sol_hom2(nz-1, 0, 0, i) ;
00305 
00306     delete resu.set_spectral_va().c ;
00307     resu.set_spectral_va().c = 0x0 ;
00308     resu.set_dzpuis(0) ;
00309     resu.set_spectral_va().ylm_i() ;
00310     
00311     return resu ; 
00312 }