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