set_expa_evol.C
// 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;
}
}
Lorene::Excision_surf::set_expa
Scalar & set_expa()
Sets the expansion function on the surface at time t (considering to protect this function).
Definition excision_surf.h:220
Lorene::Excision_surf::lapse
Scalar lapse
The lapse defined on the 3 slice.
Definition excision_surf.h:56
Lorene::Excision_surf::set_dt_expa
Scalar & set_dt_expa()
Sets the time derivative of the expansion function on the surface at time t (considering to protect t...
Definition excision_surf.h:224
Lorene::Excision_surf::set_expa_parab
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...
Definition set_expa_evol.C:20
Lorene::Excision_surf::delta_t
double delta_t
The time step for evolution in parabolic drivers.
Definition excision_surf.h:68
Lorene::Excision_surf::expa
Scalar expa
The 2d expansion, directly evolved from the initial excision with Einstein Equations.
Definition excision_surf.h:74
Lorene::Excision_surf::dt_expa
Scalar dt_expa
The time derivative of the expansion, derived from Einstein equations and arbitrary evolution.
Definition excision_surf.h:77
Lorene::Excision_surf::set_expa_hyperb
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...
Definition set_expa_evol.C:99
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::lapang
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition scalar_deriv.C:461
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::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::set_grid_point
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
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765