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
28char 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
68namespace Lorene {
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() ;
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 ;
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 ;
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}
Bases of the spectral expansions.
Definition base_val.h:322
void mult_x()
The basis is transformed as with a multiplication by .
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
Multi-domain grid.
Definition grilles.h:273
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Basic array class.
Definition tbl.h:161
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...
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64