LORENE
set_expa_evol.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 "unites.h"
12#include "excision_surf.h"
13
14
15
16
17namespace Lorene {
18//gives a new value for expansion rescaled with lapse (and its time derivative) obtained by parabolic evolution.
19//All manipulated quantities are 2-dimensional.
20void Excision_surf::set_expa_parab(double c_theta_lap, double c_theta_fin, Scalar& expa_fin){
21
22 // Definition of ff
23 // ================
24
25 // Start Mapping interpolation
26
27 if (expa.get_spectral_va().get_etat() == 0)
28 {
29 // Scalar theta = sph.theta_plus();
30 Scalar theta (lapse.get_mp()); theta = 3.; theta.std_spectral_base();
31
32 set_expa() = theta;}
33
34
35
36 Scalar thetaplus = expa;
37
38 // Interpolation for the lapse (to get a 2d quantity);
39 Scalar lapse2(thetaplus.get_mp());
40 lapse2.annule_hard();
41 lapse2.std_spectral_base();
42
43 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
44 int np = (*lapse.get_mp().get_mg()).get_np(1);
45 for (int k=0; k<np; k++)
46 for (int j=0; j<nt; j++) {
47
48
49 lapse2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
50
51 }
52
53 lapse2.std_spectral_base();
54
55
56 // End Mapping interpolation
57// cout << "convergence?" << endl;
58// cout << expa_fin3.val_grid_point(1,0,0,0) << endl;
59// cout << theta_plus3.val_grid_point(1,0,0,0) << endl;
60
61
62 Scalar ff = lapse2*(c_theta_lap*thetaplus.lapang() + c_theta_fin*(thetaplus - expa_fin));
63
64
65 // Definition of k_1
66 Scalar k_1 =delta_t*ff;
67
68 // Intermediate value of the expansion, for Runge-Kutta 2nd order scheme
69 Scalar theta_int = thetaplus + k_1; theta_int.std_spectral_base();
70
71 // Recalculation of ff with intermediate values.
72 Scalar ff_int = lapse2*(c_theta_lap*theta_int.lapang() + c_theta_fin*(theta_int - expa_fin));
73 ff_int.set_spectral_va().ylm();
74
75 // Definition of k_2
76 Scalar k_2 = delta_t*ff_int;
77
78 // Result of RK2 evolution
79 Scalar bound_theta = thetaplus + k_2;
80
81
82 // Assigns new values to the members expa() and dt_expa().
83 set_expa()=bound_theta;
84 set_dt_expa()= ff_int;
85
86 return;
87}
88
89
90
91//gives a new value for the expansion (rescaled with lapse) and its time derivative, obtained by hyperbolic evolution.
92// Parameters for the hyperbolic driver are determined by the function Excision_surf::get_evol_params_from_ID()
93// so that the expansion stays of regularity $C^{1}$ throughout.
94// The scheme used is a RK4
95// Final value is here necessarily zero for the expansion
96//All manipulated quantities are 2-dimensional.
97// Warning: no rescaling of time dimensionality by the lapse yet.
98
99void Excision_surf::set_expa_hyperb(double alpha0, double beta0, double gamma0) {
100
101 // Definition of ff
102 // ================
103
104 // Start Mapping interpolation
105 Scalar thetaplus = expa;
106 thetaplus.set_spectral_va().ylm();
107 assert (dt_expa.get_spectral_va().get_etat() != 0);
108 Scalar d_thetaplus = dt_expa;
109 d_thetaplus.set_spectral_va().ylm();
110
112 // Interpolating for the lapse into the 2-dimensional surface (if necessary...)
113 Scalar lapse2(thetaplus.get_mp());
114 lapse2.annule_hard();
115 lapse2.std_spectral_base();
116
117 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
118 int np = (*lapse.get_mp().get_mg()).get_np(1);
119 for (int k=0; k<np; k++)
120 for (int j=0; j<nt; j++) {
121 lapse2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
122 }
123
124 lapse2.std_spectral_base();
125
127
129 // Starting the RK4 1step evolution for the second-order 2d PDE in spherical harmonics.
131
132
133 //Successive derivative estimates for the scheme
134
135 //Step1
136 Scalar k_1 = d_thetaplus;
137 Scalar K_1 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus ;
138
139 thetaplus = expa + 0.5*delta_t*k_1;
140 d_thetaplus = dt_expa + 0.5*delta_t*K_1;
141
142 //Step2
143 Scalar k_2 = d_thetaplus;
144 Scalar K_2 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
145
146 thetaplus = expa + 0.5*delta_t*k_2;
147 d_thetaplus = dt_expa + 0.5*delta_t*K_2;
148
149 //Step3
150 Scalar k_3 = d_thetaplus;
151 Scalar K_3 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
152
153 thetaplus = expa + delta_t*k_3;
154 d_thetaplus = dt_expa + delta_t*K_3;
155
156 //Step4
157 Scalar k_4 = d_thetaplus;
158 Scalar K_4 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
159
160
161 // Result of RK2 evolution: assignment at evolved time.
162 thetaplus = expa + (1./6.)*delta_t*(k_1 + 2.*k_2 + 2.*k_3 + k_4);
163 d_thetaplus = dt_expa + (1./6.)*delta_t*(K_1 + 2.*K_2 + 2.*K_3 + K_4);
164
165
166 set_expa() = thetaplus;
167 set_dt_expa() = d_thetaplus;
168
169
170
171 return;
172
173}
174
175}
Scalar & set_expa()
Sets the expansion function on the surface at time t (considering to protect this function)
Scalar lapse
The lapse defined on the 3 slice.
Scalar & set_dt_expa()
Sets the time derivative of the expansion function on the surface at time t (considering to protect t...
void set_expa_parab(double c_theta_lap, double c_theta_fin, Scalar &expa_fin)
Sets a new value for expansion rescaled over lapse (and its derivative), obtained by parabolic evolut...
double delta_t
The time step for evolution in parabolic drivers.
Scalar expa
The 2d expansion, directly evolved from the initial excision with Einstein Equations.
Scalar dt_expa
The time derivative of the expansion, derived from Einstein equations and arbitrary evolution.
void set_expa_hyperb(double alph0, double beta0, double gamma0)
Sets a new value for expansion rescaled over lapse (and its derivative), obtained by hyperbolic evolu...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
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
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
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
int get_etat() const
Returns the logical state.
Definition valeur.h:726
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64