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 
29 
30 /*
31  * $Id: vector_divfree_A_1z.C,v 1.5 2016/12/05 16:18:18 j_novak Exp $
32  * $Log: vector_divfree_A_1z.C,v $
33  * Revision 1.5 2016/12/05 16:18:18 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.4 2014/10/13 08:53:45 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.3 2014/10/06 15:13:20 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.2 2009/10/23 13:18:46 j_novak
43  * Minor modifications
44  *
45  * Revision 1.1 2008/08/27 09:01:27 jl_cornou
46  * Methods for solving Dirac systems for divergence free vectors
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_1z.C,v 1.5 2016/12/05 16:18:18 j_novak Exp $
50  *
51  */
52 
53 
54 // C headers
55 #include <cstdlib>
56 #include <cassert>
57 #include <cmath>
58 
59 // Lorene headers
60 #include "metric.h"
61 #include "diff.h"
62 #include "proto.h"
63 #include "param.h"
64 
65 //----------------------------------------------------------------------------------
66 //
67 // sol_Dirac_A
68 // 1 seule zone !
69 //----------------------------------------------------------------------------------
70 
71 namespace Lorene {
72 void Vector_divfree::sol_Dirac_A_1z(const Scalar& aaa, Scalar& tilde_vr, Scalar& tilde_eta,
73  const Param* par_bc) const {
74 
75  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
76  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
77 
78  const Mg3d& mgrid = *mp_aff->get_mg() ;
79  assert(mgrid.get_type_r(0) == RARE) ;
80  if (aaa.get_etat() == ETATZERO) {
81  tilde_vr = 0 ;
82  tilde_eta = 0 ;
83  return ;
84  }
85  assert(aaa.get_etat() != ETATNONDEF) ;
86  int nz = mgrid.get_nzone() ;
87  int nzm1 = nz - 1 ;
88  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
89  int n_shell = ced ? nz-2 : nzm1 ;
90  int nz_bc = nzm1 ;
91  if (par_bc != 0x0)
92  if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
93  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
94 //#ifndef NDEBUG
95 // if (!cedbc) {
96 // assert(par_bc != 0x0) ;
97 // assert(par_bc->get_n_tbl_mod() >= 3) ;
98 // }
99 //#endif
100  int nt = mgrid.get_nt(0) ;
101  int np = mgrid.get_np(0) ;
102 
103  Scalar source = aaa ;
104  Scalar source_coq = aaa ;
105  source_coq.annule_domain(0) ;
106  if (ced) source_coq.annule_domain(nzm1) ;
107  source_coq.mult_r() ;
108  source.set_spectral_va().ylm() ;
109  source_coq.set_spectral_va().ylm() ;
110  Base_val base = source.get_spectral_base() ;
111  base.mult_x() ;
112 
113  tilde_vr.annule_hard() ;
114  tilde_vr.set_spectral_base(base) ;
115  tilde_vr.set_spectral_va().set_etat_cf_qcq() ;
116  tilde_vr.set_spectral_va().c_cf->annule_hard() ;
117  tilde_eta.annule_hard() ;
118  tilde_eta.set_spectral_base(base) ;
119  tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
120  tilde_eta.set_spectral_va().c_cf->annule_hard() ;
121 
122  Mtbl_cf sol_vr(mgrid, base) ; sol_vr.annule_hard() ;
123  Mtbl_cf sol_eta(mgrid, base) ; sol_eta.annule_hard() ;
124 
125  int l_q, m_q, base_r ;
126 
127  //---------------
128  //-- NUCLEUS ---
129  //---------------
130  {int lz = 0 ;
131  int nr = mgrid.get_nr(lz) ;
132  double alpha = mp_aff->get_alpha()[lz] ;
133  Matrice ope(2*nr, 2*nr) ;
134  ope.set_etat_qcq() ;
135 
136  for (int k=0 ; k<np+1 ; k++) {
137  for (int j=0 ; j<nt ; j++) {
138  // quantic numbers and spectral bases
139  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
140  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
141  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
142  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
143 
144  for (int lin=0; lin<nr; lin++)
145  for (int col=0; col<nr; col++)
146  ope.set(lin,col) = md(lin,col) + 2*ms(lin,col) ;
147  for (int lin=0; lin<nr; lin++)
148  for (int col=0; col<nr; col++)
149  ope.set(lin,col+nr) = -l_q*(l_q+1)*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) = -ms(lin,col) ;
153  for (int lin=0; lin<nr; lin++)
154  for (int col=0; col<nr; col++)
155  ope.set(lin+nr,col+nr) = md(lin,col) + ms(lin,col) ;
156 
157  ope *= 1./alpha ;
158  int ind1 = nr ;
159 
160  if (l_q==1) {
161  ind1 += 1 ;
162  int pari = 1 ;
163  for (int col=0 ; col<nr; col++) {
164  ope.set(nr-1,col) = pari ;
165  ope.set(nr-1,col+nr) = -pari ;
166  pari = - pari ;
167  }
168  for (int col=0 ; col<nr ; col++) {
169  ope.set(2*nr-1,col+nr)=1 ;
170  }
171  }
172 
173  else{
174  for (int col=0; col<2*nr; col++) {
175  ope.set(ind1+nr-2, col) = 0 ;
176  }
177  for (int col=nr; col<2*nr; col++)
178  ope.set(ind1+nr-2, col) = 1 ;
179  for (int col=0; col<2*nr; col++) {
180  ope.set(nr-1, col) = 0 ;
181  ope.set(2*nr-1, col) = 0 ;
182  }
183  int pari = 1 ;
184  if (base_r == R_CHEBP) {
185  for (int col=0; col<nr; col++) {
186  ope.set(nr-1, col) = pari ;
187  ope.set(2*nr-1, col+nr) = pari ;
188  pari = - pari ;
189  }
190  }
191  else { //In the odd case, the last coefficient must be zero!
192  ope.set(nr-1, nr-1) = 1 ;
193  ope.set(2*nr-1, 2*nr-1) = 1 ;
194  }
195 
196  }
197 
198  ope.set_lu() ;
199 
200  Tbl sec(2*nr) ;
201  sec.set_etat_qcq() ;
202  for (int lin=0; lin<nr; lin++)
203  sec.set(lin) = 0 ;
204  for (int lin=0; lin<nr; lin++)
205  sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
206  (lz, k, j, lin) ;
207  sec.set(2*nr-1) = 0 ;
208 
209 
210 
211 /* // BC is here
212  if ((l_q==2)&&(k==3)) {
213  sec.set(ind1+nr-2) = -5./2. ; }
214  else { sec.set(ind1+nr-2) = 0 ; }*/
215 
216 
217 
218  Tbl sol = ope.inverse(sec) ;
219 
220  for (int i=0; i<nr; i++) {
221  sol_vr.set(lz, k, j, i) = sol(i) ;
222  sol_eta.set(lz, k, j, i) = sol(i+nr) ;
223  }
224  if ((l_q==2)&&(k==3)) {
225  cout << " ========================== " << endl ;
226  cout << " Operateur " << endl ;
227  cout << " ========================== " << endl ;
228  cout << ope << endl ;
229  cout << " ========================== " << endl ;
230  cout << " Second membre " << endl ;
231  cout << " ========================== " << endl ;
232  cout << sec << endl ;
233  cout << " ========================== " << endl ;
234  cout << " Solution " << endl ;
235  cout << " ========================== " << endl ;
236  cout << sol << endl ;
237 
238  }
239  }
240  }
241  }
242  }
243 
244  Mtbl_cf& mvr = *tilde_vr.set_spectral_va().c_cf ;
245  Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
246 
247  Mtbl_cf d_vr = sol_vr ;
248  Mtbl_cf d_eta = sol_eta ;
249 
250 
251  // Loop on l and m
252  //----------------
253  for (int k=0 ; k<np+1 ; k++)
254  for (int j=0 ; j<nt ; j++) {
255  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
256  if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
257  // everything is put to the right place...
258  //----------------------------------------
259  int nr = mgrid.get_nr(0) ; //nucleus
260  for (int i=0 ; i<nr ; i++) {
261  mvr.set(0, k, j, i) = sol_vr(0, k, j, i) ;
262  meta.set(0, k, j, i) = sol_eta(0, k, j, i) ;
263  }
264  } // End of nullite_plm
265  } //End of loop on theta
266 
267 
268  if (tilde_vr.set_spectral_va().c != 0x0)
269  delete tilde_vr.set_spectral_va().c ;
270  tilde_vr.set_spectral_va().c = 0x0 ;
271  tilde_vr.set_spectral_va().ylm_i() ;
272 
273  if (tilde_eta.set_spectral_va().c != 0x0)
274  delete tilde_eta.set_spectral_va().c ;
275  tilde_eta.set_spectral_va().c = 0x0 ;
276  tilde_eta.set_spectral_va().ylm_i() ;
277 
278 }
279 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:103
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
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
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:783
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
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:560
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int get_n_int() const
Returns the number of stored int &#39;s addresses.
Definition: param.C:242
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
#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:309
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
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: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:2048
void mult_x()
The basis is transformed as with a multiplication by .
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
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: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
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:491
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607