LORENE
et_rot_mag_equil_plus.C
1 /*
2  * Function et_rot_mag::equilibrium_mag_plus
3  *
4  * Computes rotating equilibirum with a magnetic field with extended features
5  * (see file et_rot_mag.h for documentation)
6  *
7  */
8 
9 /*
10  * Copyright (c) 2012 Pablo Cerda, Michael Gabler
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 /*
33  * $Id: et_rot_mag_equil_plus.C,v 1.5 2016/12/05 16:17:54 j_novak Exp $
34  * $Log: et_rot_mag_equil_plus.C,v $
35  * Revision 1.5 2016/12/05 16:17:54 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.4 2014/10/13 08:52:58 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.3 2014/10/06 15:13:09 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.2 2013/11/25 14:03:55 j_novak
45  * Commented some variables to avoid warnings
46  *
47  * Revision 1.1 2012/08/12 17:48:35 p_cerda
48  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_equil_plus.C,v 1.5 2016/12/05 16:17:54 j_novak Exp $
52  *
53  */
54 
55 // Headers C
56 #include <cmath>
57 
58 // Headers Lorene
59 #include "et_rot_mag.h"
60 #include "param.h"
61 #include "unites.h"
62 #include "graphique.h"
63 
64 namespace Lorene {
66  const Itbl& icontrol, const Tbl& control, Tbl& diff,
67  const int initial_j,
68  const Tbl an_j,
69  Cmp (*f_j)(const Cmp&, const Tbl),
70  Cmp (*)(const Cmp& x, const Tbl),
71  const Tbl bn_j,
72  Cmp (*g_j)(const Cmp&, const Tbl),
73  Cmp (*N_j)(const Cmp& x, const Tbl),
74  const double relax_mag) {
75 
76  // Fundamental constants and units
77  // -------------------------------
78  using namespace Unites_mag ;
79 
80  // Grid parameters
81  // ---------------
82 
83  // const Mg3d* mg = mp.get_mg() ;
84  // int nz = mg->get_nzone() ; // total number of domains
85  // int nzm1 = nz - 1 ;
86 
87  // The following is required to initialize mp_prev as a Map_et:
88  //Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
89 
90 
91  // Parameters to control the iteration
92  // -----------------------------------
93 
94  int mer_max = icontrol(0) ;
95  // int mer_rot = icontrol(1) ;
96  // int mer_change_omega = icontrol(2) ;
97  // int mer_fix_omega = icontrol(3) ;
98  // int mer_mass = icontrol(4) ;
99  int mermax_poisson = icontrol(5) ;
100  // int delta_mer_kep = icontrol(6) ;
101  // int mer_mag = icontrol(7) ;
102  // int mer_change_mag = icontrol(8) ;
103  // int mer_fix_mag = icontrol(9) ;
104 
105 
106  double precis = control(0) ;
107  // double omega_ini = control(1) ;
108  // double relax = control(2) ;
109  // double relax_prev = double(1) - relax ;
110  double relax_poisson = control(3) ;
111  // double thres_adapt = control(4) ;
112  // double precis_adapt = control(5) ;
113  // double Q_ini = control(6) ;
114  // double a_j_ini = control (7) ;
115 
116  // Error indicators
117  // ----------------
118 
119  diff.set_etat_qcq() ;
120  double& diff_A_phi = diff.set(0) ;
121 
122  // Parameters for the function Map_et::adapt
123  // -----------------------------------------
124 
125  int niter ;
126  int adapt_flag = 1 ; // 1 = performs the full computation,
127  // 0 = performs only the rescaling by
128  // the factor alpha_r
129 
130 
131 
132 
133  // Parameters for the Maxwell equations
134  // -------------------------------------
135 
136  double precis_poisson = 1.e-16 ;
137 
138  Param par_poisson_At ; // For scalar At Poisson equation
139  Cmp ssjm1_At(mp) ;
140  ssjm1_At.set_etat_zero() ;
141  par_poisson_At.add_int(mermax_poisson, 0) ; // maximum number of iterations
142  par_poisson_At.add_double(relax_poisson, 0) ; // relaxation parameter
143  par_poisson_At.add_double(precis_poisson, 1) ; // required precision
144  par_poisson_At.add_int_mod(niter, 0) ; // number of iterations actually used
145  par_poisson_At.add_cmp_mod( ssjm1_At ) ;
146 
147  Param par_poisson_Avect ; // For vector Aphi Poisson equation
148 
149  Cmp ssjm1_khi_mag(ssjm1_khi) ;
150  Tenseur ssjm1_w_mag(ssjm1_wshift) ;
151 
152  par_poisson_Avect.add_int(mermax_poisson, 0) ; // maximum number of iterations
153  par_poisson_Avect.add_double(relax_poisson, 0) ; // relaxation parameter
154  par_poisson_Avect.add_double(precis_poisson, 1) ; // required precision
155  par_poisson_Avect.add_cmp_mod( ssjm1_khi_mag ) ;
156  par_poisson_Avect.add_tenseur_mod( ssjm1_w_mag ) ;
157  par_poisson_Avect.add_int_mod(niter, 0) ;
158 
159 
160  // Initializations
161  // ---------------
162 
163  // Initial magnetic quantities
164  a_j = 0;
165 
166  update_metric() ; // update of the metric coefficients
167 
168  equation_of_state() ; // update of the density, pressure, etc...
169 
170  hydro_euler() ; // update of the hydro quantities relative to the
171  // Eulerian observer
172 
173  MHD_comput() ; // update of EM contributions to stress-energy tensor
174 
175 
176  // Output files
177 
178  ofstream fichmulti("multipoles.d") ;
179 
180 
181  ofstream fichconv("convergence.d") ; // Output file for diff_A_phi
182  fichconv << "# diff_A_phi GRV2 " << endl ;
183 
184  ofstream fichfreq("frequency.d") ; // Output file for omega
185  fichfreq << "# f [Hz]" << endl ;
186 
187  ofstream fichevol("evolution.d") ; // Output file for various quantities
188  fichevol <<
189  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
190  << endl ;
191 
192  diff_A_phi = 1 ;
193  // double err_grv2 = 1 ;
194 
195 
196  A_phi = 0. ;
197  A_phi.std_base_scal() ;
198  A_t = 0.;
199  A_t.std_base_scal() ;
200  j_phi = 0.;
201  j_phi.std_base_scal() ;
202 
203  Cmp A_phi_old = A_phi;
204  Cmp A_phi_new = A_phi;
205  Cmp A_t_old = A_t;
206  Cmp A_t_new = A_t;
207  Cmp j_phi_old = j_phi;
208  Cmp j_phi_new = j_phi;
209 
210  //=========================================================================
211  // Start of iteration
212  //=========================================================================
213 
214  for(int mer=0 ; (diff_A_phi > precis) && (mer<mer_max) ; mer++ ) {
215 
216  cout << "-----------------------------------------------" << endl ;
217  cout << "step: " << mer << endl ;
218 
219  fichconv << mer ;
220  fichfreq << mer ;
221  fichevol << mer ;
222 
223  A_t_old = A_t;
224  A_phi_old = A_phi;
225  j_phi_old = j_phi;
226 
227  //-----------------------------------------------
228  // Computation of electromagnetic potentials :
229  // -------------------------------------------
230 
231  magnet_comput_plus(adapt_flag, initial_j,
232  an_j, f_j, bn_j, g_j, N_j, par_poisson_At, par_poisson_Avect) ;
233  A_t_new = A_t;
234  A_phi_new = A_phi;
235  j_phi_new = j_phi;
236 
237  A_t = relax_mag*A_t_new + (1.-relax_mag)*A_t_old ;
238  A_phi = relax_mag*A_phi_new + (1. - relax_mag)*A_phi_old ;
239 
240 
241  double diff_A_phi_n = max(abs((A_phi_new - A_phi_old))).set(0);
242  double max_Aphi = max(abs(A_phi)).set(0);
243  double diff_j_phi_n = max(abs((j_phi_new - j_phi_old))).set(0);
244  double max_jphi = max(abs(j_phi)).set(0);
245 
246  Tbl maphi = A_phi_new.multipole_spectrum();
247  // int nzmax = maphi.get_dim (1) -1;
248 
249  if (max_Aphi == 0) {
250  diff_A_phi = 100.;
251  }else{
252  diff_A_phi = diff_A_phi_n / max_Aphi ;
253  }
254  cout << mer << " "<< diff_A_phi << " " << max(A_phi).set(0) << " " << min(A_phi).set(0) << endl;
255  cout << mer << " "<< diff_j_phi_n << " " << max_jphi << endl;
256 
257 
258  fichmulti << diff_A_phi<< " "
259  <<maphi.set(0,0) << " " <<maphi.set(0,1) << " "
260  <<maphi.set(0,2) << " " <<maphi.set(0,3) << " "
261  <<maphi.set(0,4) << endl;
262 
263 
264  // des_coupe_y(A_phi, 0., nzet, "Magnetic field") ;
265  // des_coupe_y(j_phi, 0., nzet, "Current") ;
266 
267 
268 
269  fichconv << endl ;
270  fichfreq << endl ;
271  fichevol << endl ;
272  fichconv.flush() ;
273  fichfreq.flush() ;
274  fichevol.flush() ;
275 
276  } // End of main loop
277 
278  //=========================================================================
279  // End of iteration
280  //=========================================================================
281 
282  fichconv.close() ;
283  fichfreq.close() ;
284  fichevol.close() ;
285  fichmulti.close();
286 
287 }
288 }
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1145
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: etoile.h:1619
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
virtual void magnet_comput_plus(const int adapt_flag, const int initial_j, const Tbl an_j, Cmp(*f_j)(const Cmp &x, const Tbl), const Tbl bn_j, Cmp(*g_j)(const Cmp &x, const Tbl), Cmp(*N_j)(const Cmp &x, const Tbl), Param &par_poisson_At, Param &par_poisson_Avect)
Computes the electromagnetic quantities solving the Maxwell equations (6) and (7) of [Bocquet...
Lorene prototypes.
Definition: app_hor.h:67
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Definition: etoile.h:1628
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Basic integer array class.
Definition: itbl.h:122
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
Definition: et_rot_mag.h:155
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
void update_metric()
Computes metric coefficients from known potentials.
Definition: et_rot_upmetr.C:72
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:292
double a_j
Amplitude of the curent/charge function.
Definition: et_rot_mag.h:180
void equilibrium_mag_plus(const Itbl &icontrol, const Tbl &control, Tbl &diff, const int initial_j, const Tbl an_j, Cmp(*f_j)(const Cmp &x, const Tbl), Cmp(*M_j)(const Cmp &x, const Tbl), const Tbl bn_j, Cmp(*g_j)(const Cmp &x, const Tbl), Cmp(*N_j)(const Cmp &x, const Tbl), const double relax_mag)
Computes an equilibrium configuration.
Tbl multipole_spectrum()
Gives the spectrum in terms of multipolar modes l .
Definition: cmp.C:765
Parameter storage.
Definition: param.h:125
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:569
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition: et_rot_hydro.C:86
Cmp j_phi
-component of the current 4-vector
Definition: et_rot_mag.h:159
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Definition: et_rot_mag.h:150
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.
Basic array class.
Definition: tbl.h:164
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1007
Standard electro-magnetic units.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388