LORENE
vector_divfree_A_1z.C
1 /*
2  * Methods to impose the Dirac gauge: divergence-free condition.
3  *
4  * (see file sym_tensor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2006 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 char sol_Dirac_A_1z_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_1z.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $" ;
29 
30 /*
31  * $Id: vector_divfree_A_1z.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
32  * $Log: vector_divfree_A_1z.C,v $
33  * Revision 1.4 2014/10/13 08:53:45 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.3 2014/10/06 15:13:20 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.2 2009/10/23 13:18:46 j_novak
40  * Minor modifications
41  *
42  * Revision 1.1 2008/08/27 09:01:27 jl_cornou
43  * Methods for solving Dirac systems for divergence free vectors
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_1z.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
47  *
48  */
49 
50 
51 // C headers
52 #include <cstdlib>
53 #include <cassert>
54 #include <cmath>
55 
56 // Lorene headers
57 #include "metric.h"
58 #include "diff.h"
59 #include "proto.h"
60 #include "param.h"
61 
62 //----------------------------------------------------------------------------------
63 //
64 // sol_Dirac_A
65 // 1 seule zone !
66 //----------------------------------------------------------------------------------
67 
68 namespace Lorene {
69 void Vector_divfree::sol_Dirac_A_1z(const Scalar& aaa, Scalar& tilde_vr, Scalar& tilde_eta,
70  const Param* par_bc) const {
71 
72  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
73  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
74 
75  const Mg3d& mgrid = *mp_aff->get_mg() ;
76  assert(mgrid.get_type_r(0) == RARE) ;
77  if (aaa.get_etat() == ETATZERO) {
78  tilde_vr = 0 ;
79  tilde_eta = 0 ;
80  return ;
81  }
82  assert(aaa.get_etat() != ETATNONDEF) ;
83  int nz = mgrid.get_nzone() ;
84  int nzm1 = nz - 1 ;
85  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
86  int n_shell = ced ? nz-2 : nzm1 ;
87  int nz_bc = nzm1 ;
88  if (par_bc != 0x0)
89  if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
90  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
91 //#ifndef NDEBUG
92 // if (!cedbc) {
93 // assert(par_bc != 0x0) ;
94 // assert(par_bc->get_n_tbl_mod() >= 3) ;
95 // }
96 //#endif
97  int nt = mgrid.get_nt(0) ;
98  int np = mgrid.get_np(0) ;
99 
100  Scalar source = aaa ;
101  Scalar source_coq = aaa ;
102  source_coq.annule_domain(0) ;
103  if (ced) source_coq.annule_domain(nzm1) ;
104  source_coq.mult_r() ;
105  source.set_spectral_va().ylm() ;
106  source_coq.set_spectral_va().ylm() ;
107  Base_val base = source.get_spectral_base() ;
108  base.mult_x() ;
109 
110  tilde_vr.annule_hard() ;
111  tilde_vr.set_spectral_base(base) ;
112  tilde_vr.set_spectral_va().set_etat_cf_qcq() ;
113  tilde_vr.set_spectral_va().c_cf->annule_hard() ;
114  tilde_eta.annule_hard() ;
115  tilde_eta.set_spectral_base(base) ;
116  tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
117  tilde_eta.set_spectral_va().c_cf->annule_hard() ;
118 
119  Mtbl_cf sol_vr(mgrid, base) ; sol_vr.annule_hard() ;
120  Mtbl_cf sol_eta(mgrid, base) ; sol_eta.annule_hard() ;
121 
122  int l_q, m_q, base_r ;
123 
124  //---------------
125  //-- NUCLEUS ---
126  //---------------
127  {int lz = 0 ;
128  int nr = mgrid.get_nr(lz) ;
129  double alpha = mp_aff->get_alpha()[lz] ;
130  Matrice ope(2*nr, 2*nr) ;
131  ope.set_etat_qcq() ;
132 
133  for (int k=0 ; k<np+1 ; k++) {
134  for (int j=0 ; j<nt ; j++) {
135  // quantic numbers and spectral bases
136  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
137  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
138  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
139  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
140 
141  for (int lin=0; lin<nr; lin++)
142  for (int col=0; col<nr; col++)
143  ope.set(lin,col) = md(lin,col) + 2*ms(lin,col) ;
144  for (int lin=0; lin<nr; lin++)
145  for (int col=0; col<nr; col++)
146  ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
147  for (int lin=0; lin<nr; lin++)
148  for (int col=0; col<nr; col++)
149  ope.set(lin+nr,col) = -ms(lin,col) ;
150  for (int lin=0; lin<nr; lin++)
151  for (int col=0; col<nr; col++)
152  ope.set(lin+nr,col+nr) = md(lin,col) + ms(lin,col) ;
153 
154  ope *= 1./alpha ;
155  int ind1 = nr ;
156 
157  if (l_q==1) {
158  ind1 += 1 ;
159  int pari = 1 ;
160  for (int col=0 ; col<nr; col++) {
161  ope.set(nr-1,col) = pari ;
162  ope.set(nr-1,col+nr) = -pari ;
163  pari = - pari ;
164  }
165  for (int col=0 ; col<nr ; col++) {
166  ope.set(2*nr-1,col+nr)=1 ;
167  }
168  }
169 
170  else{
171  for (int col=0; col<2*nr; col++) {
172  ope.set(ind1+nr-2, col) = 0 ;
173  }
174  for (int col=nr; col<2*nr; col++)
175  ope.set(ind1+nr-2, col) = 1 ;
176  for (int col=0; col<2*nr; col++) {
177  ope.set(nr-1, col) = 0 ;
178  ope.set(2*nr-1, col) = 0 ;
179  }
180  int pari = 1 ;
181  if (base_r == R_CHEBP) {
182  for (int col=0; col<nr; col++) {
183  ope.set(nr-1, col) = pari ;
184  ope.set(2*nr-1, col+nr) = pari ;
185  pari = - pari ;
186  }
187  }
188  else { //In the odd case, the last coefficient must be zero!
189  ope.set(nr-1, nr-1) = 1 ;
190  ope.set(2*nr-1, 2*nr-1) = 1 ;
191  }
192 
193  }
194 
195  ope.set_lu() ;
196 
197  Tbl sec(2*nr) ;
198  sec.set_etat_qcq() ;
199  for (int lin=0; lin<nr; lin++)
200  sec.set(lin) = 0 ;
201  for (int lin=0; lin<nr; lin++)
202  sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
203  (lz, k, j, lin) ;
204  sec.set(2*nr-1) = 0 ;
205 
206 
207 
208 /* // BC is here
209  if ((l_q==2)&&(k==3)) {
210  sec.set(ind1+nr-2) = -5./2. ; }
211  else { sec.set(ind1+nr-2) = 0 ; }*/
212 
213 
214 
215  Tbl sol = ope.inverse(sec) ;
216 
217  for (int i=0; i<nr; i++) {
218  sol_vr.set(lz, k, j, i) = sol(i) ;
219  sol_eta.set(lz, k, j, i) = sol(i+nr) ;
220  }
221  if ((l_q==2)&&(k==3)) {
222  cout << " ========================== " << endl ;
223  cout << " Operateur " << endl ;
224  cout << " ========================== " << endl ;
225  cout << ope << endl ;
226  cout << " ========================== " << endl ;
227  cout << " Second membre " << endl ;
228  cout << " ========================== " << endl ;
229  cout << sec << endl ;
230  cout << " ========================== " << endl ;
231  cout << " Solution " << endl ;
232  cout << " ========================== " << endl ;
233  cout << sol << endl ;
234 
235  }
236  }
237  }
238  }
239  }
240 
241  Mtbl_cf& mvr = *tilde_vr.set_spectral_va().c_cf ;
242  Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
243 
244  Mtbl_cf d_vr = sol_vr ;
245  Mtbl_cf d_eta = sol_eta ;
246 
247 
248  // Loop on l and m
249  //----------------
250  for (int k=0 ; k<np+1 ; k++)
251  for (int j=0 ; j<nt ; j++) {
252  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
253  if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
254  // everything is put to the right place...
255  //----------------------------------------
256  int nr = mgrid.get_nr(0) ; //nucleus
257  for (int i=0 ; i<nr ; i++) {
258  mvr.set(0, k, j, i) = sol_vr(0, k, j, i) ;
259  meta.set(0, k, j, i) = sol_eta(0, k, j, i) ;
260  }
261  } // End of nullite_plm
262  } //End of loop on theta
263 
264 
265  if (tilde_vr.set_spectral_va().c != 0x0)
266  delete tilde_vr.set_spectral_va().c ;
267  tilde_vr.set_spectral_va().c = 0x0 ;
268  tilde_vr.set_spectral_va().ylm_i() ;
269 
270  if (tilde_eta.set_spectral_va().c != 0x0)
271  delete tilde_eta.set_spectral_va().c ;
272  tilde_eta.set_spectral_va().c = 0x0 ;
273  tilde_eta.set_spectral_va().ylm_i() ;
274 
275 }
276 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1294
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:100
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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:712
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:392
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_n_int() const
Returns the number of stored int &#39;s addresses.
Definition: param.C:239
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Matrix handling.
Definition: matrice.h:152
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition: scalar.C:797
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
Bases of the spectral expansions.
Definition: base_val.h:322
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
Affine radial mapping.
Definition: map.h:2027
void mult_x()
The basis is transformed as with a multiplication by .
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
void sol_Dirac_A_1z(const Scalar &aaa, Scalar &eta, Scalar &vr, const Param *par_bc=0x0) const
Solves a one-domain system of two-coupled first-order PDEs obtained from the divergence-free conditio...
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Class for the elementary differential operator division by (see the base class Diff )...
Definition: diff.h:329
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601