1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
// Header Lorene:
#include "nbr_spx.h"
#include "utilitaires.h"
#include "graphique.h"
#include "math.h"
#include "metric.h"
#include "param.h"
#include "param_elliptic.h"
#include "tensor.h"
#include "unites.h"
#include "excision_surf.h"
namespace Lorene {
//gives a new value for expansion rescaled with lapse (and its time derivative) obtained by parabolic evolution.
//All manipulated quantities are 2-dimensional.
void Excision_surf::set_expa_parab(double c_theta_lap, double c_theta_fin, Scalar& expa_fin){
// Definition of ff
// ================
// Start Mapping interpolation
if (expa.get_spectral_va().get_etat() == 0)
{
// Scalar theta = sph.theta_plus();
Scalar theta (lapse.get_mp()); theta = 3.; theta.std_spectral_base();
set_expa() = theta;}
Scalar thetaplus = expa;
// Interpolation for the lapse (to get a 2d quantity);
Scalar lapse2(thetaplus.get_mp());
lapse2.annule_hard();
lapse2.std_spectral_base();
int nt = (*lapse.get_mp().get_mg()).get_nt(1);
int np = (*lapse.get_mp().get_mg()).get_np(1);
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
lapse2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
}
lapse2.std_spectral_base();
// End Mapping interpolation
// cout << "convergence?" << endl;
// cout << expa_fin3.val_grid_point(1,0,0,0) << endl;
// cout << theta_plus3.val_grid_point(1,0,0,0) << endl;
Scalar ff = lapse2*(c_theta_lap*thetaplus.lapang() + c_theta_fin*(thetaplus - expa_fin));
// Definition of k_1
Scalar k_1 =delta_t*ff;
// Intermediate value of the expansion, for Runge-Kutta 2nd order scheme
Scalar theta_int = thetaplus + k_1; theta_int.std_spectral_base();
// Recalculation of ff with intermediate values.
Scalar ff_int = lapse2*(c_theta_lap*theta_int.lapang() + c_theta_fin*(theta_int - expa_fin));
ff_int.set_spectral_va().ylm();
// Definition of k_2
Scalar k_2 = delta_t*ff_int;
// Result of RK2 evolution
Scalar bound_theta = thetaplus + k_2;
// Assigns new values to the members expa() and dt_expa().
set_expa()=bound_theta;
set_dt_expa()= ff_int;
return;
}
//gives a new value for the expansion (rescaled with lapse) and its time derivative, obtained by hyperbolic evolution.
// Parameters for the hyperbolic driver are determined by the function Excision_surf::get_evol_params_from_ID()
// so that the expansion stays of regularity $C^{1}$ throughout.
// The scheme used is a RK4
// Final value is here necessarily zero for the expansion
//All manipulated quantities are 2-dimensional.
// Warning: no rescaling of time dimensionality by the lapse yet.
void Excision_surf::set_expa_hyperb(double alpha0, double beta0, double gamma0) {
// Definition of ff
// ================
// Start Mapping interpolation
Scalar thetaplus = expa;
thetaplus.set_spectral_va().ylm();
assert (dt_expa.get_spectral_va().get_etat() != 0);
Scalar d_thetaplus = dt_expa;
d_thetaplus.set_spectral_va().ylm();
//////////////////////////////////////////////////////////////////////:
// Interpolating for the lapse into the 2-dimensional surface (if necessary...)
Scalar lapse2(thetaplus.get_mp());
lapse2.annule_hard();
lapse2.std_spectral_base();
int nt = (*lapse.get_mp().get_mg()).get_nt(1);
int np = (*lapse.get_mp().get_mg()).get_np(1);
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
lapse2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
}
lapse2.std_spectral_base();
/////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////:
// Starting the RK4 1step evolution for the second-order 2d PDE in spherical harmonics.
//////////////////////////////////////////////////////////////////////////////////
//Successive derivative estimates for the scheme
//Step1
Scalar k_1 = d_thetaplus;
Scalar K_1 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus ;
thetaplus = expa + 0.5*delta_t*k_1;
d_thetaplus = dt_expa + 0.5*delta_t*K_1;
//Step2
Scalar k_2 = d_thetaplus;
Scalar K_2 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
thetaplus = expa + 0.5*delta_t*k_2;
d_thetaplus = dt_expa + 0.5*delta_t*K_2;
//Step3
Scalar k_3 = d_thetaplus;
Scalar K_3 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
thetaplus = expa + delta_t*k_3;
d_thetaplus = dt_expa + delta_t*K_3;
//Step4
Scalar k_4 = d_thetaplus;
Scalar K_4 = beta0*d_thetaplus + alpha0*d_thetaplus.lapang() + gamma0*thetaplus;
// Result of RK2 evolution: assignment at evolved time.
thetaplus = expa + (1./6.)*delta_t*(k_1 + 2.*k_2 + 2.*k_3 + k_4);
d_thetaplus = dt_expa + (1./6.)*delta_t*(K_1 + 2.*K_2 + 2.*K_3 + K_4);
set_expa() = thetaplus;
set_dt_expa() = d_thetaplus;
return;
}
}
|
