LORENE
boundfree_tensBC.C
1
2// Header Lorene:
3#include "nbr_spx.h"
4#include "utilitaires.h"
5#include "graphique.h"
6#include "math.h"
7#include "metric.h"
8#include "param.h"
9#include "param_elliptic.h"
10#include "tensor.h"
11#include "diff.h"
12#include "proto.h"
13
14namespace Lorene {
15// Resolution of tensorial equation (N^2/Psi^4)Delta(hij) - LbLbhij = Sij, using degenerate elliptic solver.
16// Here assumption is made that no boundary condition has to be enforced, mainly Beta^i*s_i = N/psi^2
17
18Sym_tensor boundfree_tensBC ( Sym_tensor source, Vector Beta, Scalar Psi, Scalar Nn, Sym_tensor hij_guess, double precision, int loopmax) {
19 cout << "================================================================" << endl;
20 cout << "STARTING THE SUBITERATION FOR THE CONFORMAL METRIC" << endl;
21 cout << " iteration parameters are the following: " << endl;
22 cout << " convergence precision required:" << precision << endl;
23 cout << " max number of global steps :" << loopmax << endl;
24 cout << "================================================================" << endl;
25
26 // Construction of a multi-grid (Mg3d)
27 // -----------------------------------
28 const int nz = (*source.get_mp().get_mg()).get_nzone(); // Number of domains
29 int nr = (*source.get_mp().get_mg()).get_nr(1); // Number of collocation points in r in each domain
30 const Coord& rr = source.get_mp().r;
31 Scalar rrr (source.get_mp()) ;
32 rrr = rr ;
33 rrr.std_spectral_base();
34
35 const Metric_flat& mets = (source.get_mp()).flat_met_spher() ;
36
38
39 Sym_tensor hij = hij_guess;
40 Sym_tensor hij_new = hij;
41
42 Scalar n2sp4 = (Nn/(Psi*Psi))*(Nn/(Psi*Psi)) ; // Scale factor in front of Poisson Equation.
43 n2sp4.std_spectral_base();
44
45 // Resolution variables
46 Scalar khi = hij(1,1);
47 if (khi.get_etat() == ETATZERO) {
48 khi.annule_hard() ;
49 khi.set_dzpuis(4) ;
50 khi.std_spectral_base() ;
51 }
52 khi.set_spectral_va().ylm();
53 khi.mult_r();
54 khi.mult_r_dzpuis(1);
55 Scalar mmu = hij.mu();
56 if (mmu.get_etat() == ETATZERO) {
57 mmu.annule_hard() ;
58 mmu.std_spectral_base() ;
59 }
60 mmu.set_spectral_va().ylm();
61 Scalar etta = hij.eta();
62 if (etta.get_etat() == ETATZERO) {
63 etta.annule_hard() ;
64 etta.std_spectral_base() ;
65 }
66 etta.set_spectral_va().ylm();
67 Scalar Aa = hij.compute_A();
68 if (Aa.get_etat() == ETATZERO) {
69 Aa.annule_hard() ;
70 Aa.std_spectral_base() ;
71 }
72 Aa.set_spectral_va().ylm();
73 Scalar Bt = hij.compute_tilde_B();
74 if (Bt.get_etat() == ETATZERO) {
75 Bt.annule_hard() ;
76 Bt.std_spectral_base() ;
77 }
78 Bt.set_spectral_va().ylm();
79 Scalar hh = hij.trace(mets);
80 if (hh.get_etat() == ETATZERO) {
81 hh.annule_hard() ;
82 hh.std_spectral_base() ;
83 }
84
85 //Fitting scalar
86 Scalar fit1(source.get_mp()); fit1.set_etat_qcq();
87
88 // For storing the result of inversion
89 Scalar Aanew (source.get_mp()); Aanew.annule_hard(); Aanew.std_spectral_base();
90 Scalar Btnew (source.get_mp()); Btnew.annule_hard(); Btnew.std_spectral_base();
91
92
93 // Construction of sources for the next iteration
94
95 Sym_tensor LbLbhij = (hij.derive_lie(Beta)).derive_lie(Beta);
96 hij.annule_domain(0);
97 LbLbhij.annule_domain(0);
98 LbLbhij.inc_dzpuis(1);
99
100 Sym_tensor source_hij = source/n2sp4;
101 Scalar sourcetrace = source_hij.trace(mets);
102 Scalar Bttrace = source_hij.compute_tilde_B();
103
104 Sym_tensor source_hij2 = LbLbhij/n2sp4;
105
106 Scalar Btsource2 = source_hij2.compute_tilde_B();
107
108 Scalar source_Bt2 = Bttrace + Btsource2;
109
110 source_hij = source_hij + source_hij2;
111
112 Scalar r2LbLbrr = LbLbhij(1,1);
113 r2LbLbrr.mult_r();
114 r2LbLbrr.mult_r();
115 Scalar LbLbmu = LbLbhij.mu();
116
117 Scalar source_khi2 = source_hij(1,1);
118 source_khi2.mult_r();
119 source_khi2.mult_r();
120
121 Scalar source_mu2 = source_hij.mu();
122 Scalar source_eta2 = source_hij.eta();
123
124 Scalar source_A2 = source_hij.compute_A();
125
126
127 source_khi2.annule_domain(0);
128 source_mu2.annule_domain(0);
129 source_A2.annule_domain(0);
130 source_Bt2.annule_domain(0);
131 source_eta2.annule_domain(0);
132
133 // source_A2.set_spectral_va().set_base( Aa.set_spectral_va().get_base());
134 // source_Bt2.set_spectral_va().set_base( Bt.set_spectral_va().get_base());
135
136
137
139 Scalar Betacarre = (Beta(1)*Beta(1))/n2sp4 ;
140
141 double fitd1 = (Betacarre.val_grid_point(1,0,0,nr-1) - Betacarre.val_grid_point(1,0,0,0))/(rrr.val_grid_point(1,0,0,nr-1) - rrr.val_grid_point(1,0,0,0)) ;
142
143// double error = 0.; // Voluntary error on fit.
144// fitd1 += error;
145
146
147 int nrint = (nr-1)/2 ;
148 double ampl_r = (rrr.val_grid_point(1,0,0, nr -1) - rrr.val_grid_point(1,0,0 ,0))/2.;
149
150 Scalar approx(source.get_mp());
151 approx.annule_hard();
152 approx.std_spectral_base();
153
154
155
156 fit1 = fitd1*(rrr-1.) +1.;
157 // First order approximation
158 approx.set_domain(1)= fit1.set_domain(1); // MAKE PARTICULAR FOR ETATUN; DECLARE FIT1?2?3 FIRST TO BE CLEAN.
159
160
161
162 //Second order approximation
163 Scalar firststep = Betacarre - approx;
164
165 double ampli = firststep.val_grid_point(1,0,0,nrint);
166 double fit2d1 = - ampli/(ampl_r* ampl_r);
167
168
169
170 approx.set_domain(1) += (fit2d1*(rrr - 1.)*(rrr - rrr.val_grid_point(1,0,0,nr-1))).set_domain(1);
171
172
173 double fit0d2 = approx.val_grid_point(1,0,0,nr -1);
174 double fit1d2 = (Betacarre.val_grid_point(2,0,0,nr-1) - fit0d2)/(rrr.val_grid_point(2,0,0,nr-1)- rrr.val_grid_point(2,0,0,0));
175 double fit0d3 = Betacarre.val_grid_point(3,0,0,0);
176 double fit1d3 = ( - fit0d3)/(rrr.val_grid_point(3,0,0,nr-1)- rrr.val_grid_point(3,0,0,0));
177
178 approx.set_domain(2) = (fit0d2 + fit1d2*(rrr - rrr.val_grid_point(2,0,0,0))).set_domain(2);
179 approx.set_domain(3) = (fit0d3 + fit1d3*(rrr - rrr.val_grid_point(3,0,0,0))).set_domain(3);
180
181
182
183 for(int ii=1; ii<=3; ii++){
184
185 source_khi2.set_domain(ii) += (-approx*(khi.dsdr().dsdr())).set_domain(ii);
186
187 source_mu2.set_domain(ii) += (-approx*(mmu.dsdr().dsdr())).set_domain(ii);
188
189 source_eta2.set_domain(ii) += (-approx*(etta.dsdr().dsdr())).set_domain(ii);
190
191 source_A2.set_domain(ii) += (-approx*(Aa.dsdr().dsdr())).set_domain(ii);
192
193 source_Bt2.set_domain(ii) += (-approx*(Bt.dsdr().dsdr())).set_domain(ii);
194
195 }
196
197
198
199 // Convergence markers
200 Scalar Aa_old(source.get_mp());
201 Scalar Bt_old(source.get_mp());
202
203
204 // Parameters for the iteration
205 cout <<"==================================================================================" << endl;
206 cout << "amplitude for the tensor equation source (used as scaling for convergence marker)" << endl;
207 cout <<"==================================================================================" << endl;
208 double scale = max(maxabs(source));
209 double diff_ent = 0.15 ; // Initialisation of the difference marker between two iterations on some value
210
215
216 for (int mer = 0; (diff_ent>precision) && (mer< loopmax) ; mer++) {
217 double relax = 0.15; // Global relaxation parameter
218
219// Résolution for variables A and tilde(B)
220 tensorelliptic (source_A2, Aanew , fitd1, fit2d1, fit0d2, fit1d2, fit0d3, fit1d3);
221 tensorellipticBt (source_Bt2, Btnew, fitd1, fit2d1, fit0d2, fit1d2, fit0d3, fit1d3);
222
223
224 Aa = Aanew ;
225 Bt = Btnew ; // No relaxation on these variables; it will be done globally on the tensor.
226
229
230 Scalar essaiA = Aanew.laplacian() - approx*(Aanew.dsdr().dsdr());
231 essaiA.annule_domain(0);
232 Scalar bobofA = essaiA -source_A2 ;
233 bobofA.set_spectral_va().ylm() ;
234 // maxabs (bobofA);
235 // bobofA.spectral_display();
236 bobofA.set_spectral_va().ylm_i();
237
239
240 //-------------------------------------------
241 // Retrieving the tensor hij by Dirac inversion
242 //------------------------------------------
243
245 // Magnetic part
246
247 Scalar mmuA = mmu;
248 Scalar mmuAsr = mmuA; mmuAsr.div_r();
249 Scalar xxA(source.get_mp()); xxA.annule_hard(); xxA.std_spectral_base();
250
251 const Scalar AA = Aa;
252 Param* par_bc = 0x0;
253 Sym_tensor_trans hijt(source.get_mp(), source.get_mp().get_bvect_spher(), mets);
254 hijt = hij;
255 hijt.sol_Dirac_A2 ( AA , mmuAsr, xxA, source_mu2, par_bc);
256
257
258 // Monitoring A and Hmu;
259 Scalar musrsr = mmuAsr; musrsr.div_r_dzpuis(2);
260 Aanew = xxA.dsdr() - musrsr;
261 Scalar xxsr = xxA; xxsr.div_r_dzpuis(2);
262 // Scalar Hmut = 3.*musrsr + mmuAsr.dsdr() + 2.*xxsr + xxsr.lapang();
263
265 //Electric part
266
267 Scalar hrrBC (source.get_mp()); hrrBC.annule_hard();
268 Scalar wwBC(source.get_mp()); wwBC.annule_hard();
269 Scalar tilde_etaBC(source.get_mp()); tilde_etaBC.annule_hard();
270
271 Scalar source_khi3 = source_khi2;
272 source_khi3.annule_domain(nz-1);
273 source_khi3 += 2*hh;
274 source_khi3.set_spectral_va().ylm();
275
276 hijt.sol_Dirac_BC3( Bt, hh, hrrBC , tilde_etaBC , wwBC , source_khi3, source_eta2, par_bc, par_bc) ;
277
278 // Tensor reconstruction
279 hij_new.set_auxiliary(hrrBC, tilde_etaBC, mmuAsr, wwBC, xxA, hh -hrrBC);
280
281 hij= relax*hij_new + (1 - relax)*hij ; // Global relaxation (opposite to relaxation on resolution variables).
282
283 // Calculation of updated trace.
284 hh = (1 + hij(1,1))*( hij(2,3)*hij(2,3) - hij(2,2)*hij(3,3) )
285 + hij(1,2)*hij(1,2)*(1 + hij(3,3))
286 + hij(1,3)*hij(1,3)*(1 + hij(2,2))
287 - hij(1,1)*(hij(2,2) + hij(3,3)) - 2*hij(1,2)*hij(1,3)*hij(2,3) ;
288
289 khi = hij(1,1);
290 khi.mult_r();
291 khi.mult_r_dzpuis(0);
292 mmu = hij.mu();
293 etta = hij.eta();
294
295 Aa = hij.compute_A();
296 Bt = hij.compute_tilde_B();
297
298 if (mer >=1){
299 Aa.set_spectral_va().ylm();
300 Bt.set_spectral_va().ylm();
301 Aa_old.set_spectral_va().ylm();
302 Bt_old.set_spectral_va().ylm();
303
304 diff_ent = max(maxabs (Bt - Bt_old))/scale; // Convergence marker
305 }
306 Aa_old = Aa;
307 Bt_old = Bt;
308
309
310 // Update of sources for the next loop
311 LbLbhij = (hij.derive_lie(Beta)).derive_lie(Beta);
312 LbLbhij.inc_dzpuis(1);
313 r2LbLbrr = LbLbhij(1,1);
314 r2LbLbrr.mult_r();
315 r2LbLbrr.mult_r();
316 LbLbmu = LbLbhij.mu();
317
318
319 source_hij = source/n2sp4;
320
321 Bttrace = source_hij.compute_tilde_B();
322
323 source_hij2 = LbLbhij/n2sp4;
324
325 Btsource2 = source_hij2.compute_tilde_B();
326
327 source_Bt2 = Bttrace + Btsource2;
328
329 source_hij = source_hij + source_hij2;
330
331 source_khi2 = source_hij(1,1);
332 source_khi2.mult_r();
333 source_khi2.mult_r();
334 source_mu2 = source_hij.mu();
335 source_eta2 = source_hij.eta();
336 source_A2= source_hij.compute_A();
337
338
339 source_A2.set_spectral_va().set_base( Aa.set_spectral_va().get_base());
340 source_Bt2.set_spectral_va().set_base( Bt.set_spectral_va().get_base());
341
342 for(int ii=1; ii<=3; ii++){
343 source_khi2.set_domain(ii) += (-approx*(khi.dsdr().dsdr())).set_domain(ii);
344 source_mu2.set_domain(ii) += (-approx*(mmu.dsdr().dsdr())).set_domain(ii);
345 source_eta2.set_domain(ii) += (-approx*(etta.dsdr().dsdr())).set_domain(ii);
346 source_A2.set_domain(ii) += (-approx*(Aa.dsdr().dsdr())).set_domain(ii);
347 source_Bt2.set_domain(ii) += (-approx*(Bt.dsdr().dsdr())).set_domain(ii);
348 }
349 source_khi2.annule_domain(0);
350 source_mu2.annule_domain(0);
351 source_eta2.annule_domain(0);
352 source_A2.annule_domain(0);
353 source_Bt2.annule_domain(0);
354
355
356
357// cout << "real A resolution?" << endl;
358// Sym_tensor mucorrect = LbLbhij + source;
359// mucorrect = mucorrect/n2sp4;
360// Scalar AAs = mucorrect.compute_A();
361// AAs.annule_domain(nz-1);
362// Aa.annule_domain(nz-1);
363// Scalar Aa2 = Aa.laplacian();
364// Scalar voirA = AAs - Aa2;
365
366// voirA.set_spectral_va().ylm_i();
367// // des_meridian(voirA, 1.00001*Rhor, Rout,"diffA", 1);
368// voirA.set_spectral_va().ylm();
369
370
371// maxabs (voirA);
372// voirA.spectral_display("voirA", 1.e-9);
373
374
375 cout << "diff_ent" << diff_ent << endl;
376 cout << mer << endl;
377
378
379 }
380
381 return hij;
382 }
383
384
385
386
387
388
389}
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:147
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition scalar_deriv.C:436
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Lorene::Scalar::div_r
void div_r()
Division by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:133
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Sym_tensor_trans
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:608
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Sym_tensor::derive_lie
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
Lorene::maxabs
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Definition tensor_calculus_ext.C:606
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::dsdr
virtual void dsdr(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765