26char wave_utilities_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/wave_utilities.C,v 1.11 2014/10/13 08:53:31 j_novak Exp $" ;
75void tilde_laplacian(
const Scalar& B_in, Scalar& tilde_lap,
int dl) {
77 if (B_in.get_etat() == ETATZERO) {
78 tilde_lap.set_etat_zero() ;
81 assert(B_in.check_dzpuis(0)) ; ;
83 tilde_lap = B_in.laplacian(0) ;
86 assert(B_in.get_etat() != ETATNONDEF) ;
87 const Map_af* map =
dynamic_cast<const Map_af*
>(&B_in.get_mp()) ;
90 tilde_lap = 2*B_in.dsdr() ;
91 tilde_lap.div_r_dzpuis(3) ;
92 tilde_lap += B_in.dsdr().dsdr() ;
93 tilde_lap.dec_dzpuis() ;
94 tilde_lap.set_spectral_va().ylm() ;
95 Scalar B_over_r2 = B_in ;
96 B_over_r2.div_r_dzpuis(1) ;
97 B_over_r2.div_r_dzpuis(2) ;
98 B_over_r2.set_spectral_va().ylm() ;
100 const Base_val& base = B_in.get_spectral_base() ;
101 const Mg3d& mg = *map->get_mg() ;
102 int nz = mg.get_nzone() ;
103 int l_q, m_q, base_r ;
104 for (
int lz=0; lz<nz; lz++) {
105 if (B_in.domain(lz).get_etat() == ETATZERO) {
106 tilde_lap.set_spectral_va().c_cf->set(lz).set_etat_zero() ;
109 for (
int k=0; k<mg.get_np(lz)+2; k++)
110 for (
int j=0; j<mg.get_nt(lz); j++) {
111 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
114 for (
int i=0; i<mg.get_nr(lz); i++)
115 tilde_lap.set_spectral_va().c_cf->set(lz, k, j, i)
117 (*B_over_r2.get_spectral_va().c_cf)(lz, k, j, i) ;
122 if (tilde_lap.get_spectral_va().c != 0x0) {
123 delete tilde_lap.set_spectral_va().c ;
124 tilde_lap.set_spectral_va().c = 0x0 ;
126 tilde_lap.dec_dzpuis(2) ;
136void runge_kutta3_wave_sys(
double dt,
const Scalar& fff,
const Scalar& phi,
137 Scalar& fnext, Scalar& phinext,
int dl) {
139 const Map& map = fff.get_mp() ;
141 Scalar dk1(map) ; tilde_laplacian(fff, dk1, dl) ;
142 Scalar y1 = fff + 0.5*dt*k1 ;
143 Scalar dy1 = phi + 0.5*dt*dk1 ;
144 Scalar k2 = dy1 ; Scalar dk2(map) ; tilde_laplacian(y1, dk2, dl) ;
145 Scalar y2 = fff - dt*k1 + 2*dt*k2 ;
146 Scalar dy2 = phi - dt*dk1 + 2*dt*dk2 ;
148 Scalar dk3(map) ; tilde_laplacian(y2, dk3, dl) ;
149 fnext = fff + dt*(k1 + 4*k2 + k3)/6. ;
150 phinext = phi + dt*(dk1 + 4*dk2 + dk3)/6. ;
155void initialize_outgoing_BC(
int nz_bound,
const Scalar& phi,
const Scalar& dphi,
158 Scalar source_xi = phi ;
159 source_xi.div_r_dzpuis(2) ;
160 source_xi += phi.dsdr() ;
161 source_xi.dec_dzpuis(2) ;
163 if (source_xi.get_etat() == ETATZERO)
166 source_xi.set_spectral_va().ylm() ;
168 const Base_val& base_x = source_xi.get_spectral_base() ;
169 int np2 = xij.get_dim(1) ;
170 int nt = xij.get_dim(0) ;
171 assert (source_xi.get_mp().get_mg()->get_np(nz_bound) + 2 == np2 ) ;
172 assert (source_xi.get_mp().get_mg()->get_nt(nz_bound) == nt ) ;
174 int l_q, m_q, base_r ;
175 for (
int k=0; k<np2; k++)
176 for (
int j=0; j<nt; j++) {
177 base_x.give_quant_numbers(nz_bound, k, j, m_q, l_q, base_r) ;
179 = source_xi.get_spectral_va().c_cf->val_out_bound_jk(nz_bound, j, k) ;
195void evolve_outgoing_BC(
double dt,
int nz_bound,
const Scalar& phi, Scalar& sphi,
196 Tbl& xij, Tbl& xijm1, Tbl& ccc,
int dl) {
198 const Map* map = &phi.get_mp() ;
199 const Map_af* mp_aff =
dynamic_cast<const Map_af*
>(map) ;
200 assert(mp_aff != 0x0) ;
202 const Mg3d& grid = *mp_aff->get_mg() ;
204 int nz = grid.get_nzone() ;
205 assert(nz_bound < nz) ;
206 assert(phi.get_etat() != ETATZERO) ;
207 assert(sphi.get_etat() != ETATZERO) ;
209 int np2 = grid.get_np(nz_bound) + 2 ;
210 int nt = grid.get_nt(nz_bound) ;
211 assert(xij.get_ndim() == 2) ;
212 assert(xijm1.get_ndim() == 2) ;
213 assert(ccc.get_ndim() == 2) ;
214 assert(xij.get_dim(0) == nt) ;
215 assert(xij.get_dim(1) == np2) ;
216 assert(xijm1.get_dim(0) == nt) ;
217 assert(xijm1.get_dim(1) == np2) ;
218 assert(ccc.get_dim(0) == nt) ;
219 assert(ccc.get_dim(1) == np2) ;
221 double Rmax = mp_aff->get_alpha()[nz_bound] + mp_aff->get_beta()[nz_bound] ;
223 Scalar source_xi = phi ;
224 int dzp = ( source_xi.get_dzpuis() == 0 ? 2 : source_xi.get_dzpuis()+1 ) ;
225 source_xi.div_r_dzpuis(dzp) ;
226 source_xi -= phi.dsdr() ;
227 source_xi.set_spectral_va().ylm() ;
228 sphi.set_spectral_va().ylm() ;
229 const Base_val& base = sphi.get_spectral_base() ;
230 int l_q, m_q, base_r ;
231 for (
int k=0; k<np2; k++)
232 for (
int j=0; j<nt; j++) {
233 base.give_quant_numbers(nz_bound, k, j, m_q, l_q, base_r) ;
236 double fact = 8*Rmax*Rmax + dt*dt*(6+3*l_q*(l_q+1)) + 12*Rmax*dt ;
237 double souphi = -4*dt*dt*l_q*(l_q+1)*
238 source_xi.get_spectral_va().c_cf->val_out_bound_jk(nz_bound, j, k) ;
239 double xijp1 = ( 16*Rmax*Rmax*xij(k,j) -
240 (fact - 24*Rmax*dt)*xijm1(k,j)
242 ccc.set(k, j) = xijp1
243 + sphi.get_spectral_va().c_cf->val_out_bound_jk(nz_bound, j, k) ;
244 xijm1.set(k,j) = xij(k,j) ;
245 xij.set(k,j) = xijp1 ;
251void Dirichlet_BC_AtB(
const Evolution_std<Sym_tensor>& hb_evol,
252 const Evolution_std<Sym_tensor>& dhb_evol, Tbl& ccA, Tbl& ccB) {
254 int iter = hb_evol.j_max() ;
255 assert(dhb_evol.j_max() == iter) ;
257 Scalar mu_ddot = dhb_evol.time_derive(iter,3).mu() ;
259 Tbl ddmu = mu_ddot.tbl_out_bound(0,
true) ;
260 int nt = ddmu.get_dim(0) ;
261 int np2 = ddmu.get_dim(1) ;
262 const Base_val& base = mu_ddot.get_spectral_base() ;
263 int l_q, m_q, base_r ;
265 for (
int k=0; k<np2; k++) {
266 for (
int j=0; j<nt; j++) {
267 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
269 ccA.set(k,j) = ddmu(k,j) / double(l_q*(l_q+1)-2) ;
273 Scalar hrr_ddot = dhb_evol.time_derive(iter,3)(1,1) ;
274 Tbl ddhrr = hrr_ddot.tbl_out_bound(0,
true) ;
275 Scalar eta_ddot = dhb_evol.time_derive(iter,3).eta() ;
276 Tbl ddeta = eta_ddot.tbl_out_bound(0,
true) ;
277 const Base_val& base2 = hrr_ddot.get_spectral_base() ;
279 const Map& map = hrr_ddot.get_mp() ;
280 const Map_radial* mp_rad =
dynamic_cast<const Map_radial*
>(&map) ;
281 assert(mp_rad != 0x0) ;
282 for (
int k=0; k<np2; k++) {
283 for (
int j=0; j<nt; j++) {
284 base2.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
286 ccB.set(k,j) = (double(l_q+1)*ddeta(k,j)
287 + ddhrr(k,j)*mp_rad->val_r_jk(0, 1., j, k))
288 /
double((l_q+1)*(l_q-1)) ;
296void Dirichlet_BC_Amu(
const Evolution_std<Vector>& vb_evol,
297 const Evolution_std<Vector>& dvb_evol, Tbl& ccA, Tbl& ccmu) {
299 int iter = vb_evol.
j_max() ;
300 assert(dvb_evol.j_max() == iter) ;
302 Scalar vr_ddot = dvb_evol.time_derive(iter,3)(1) ;
304 Tbl ddvr = vr_ddot.tbl_out_bound(0,
true) ;
305 int nt = ddvr.get_dim(0) ;
306 int np2 = ddvr.get_dim(1) ;
307 const Base_val& base = vr_ddot.get_spectral_base() ;
308 int l_q, m_q, base_r ;
311 Scalar mu_b = vb_evol[iter].mu();
312 ccmu = mu_b.tbl_out_bound(0,
true);
313 const Map& map = vr_ddot.get_mp();
314 const Map_radial* mp_rad =
dynamic_cast<const Map_radial*
>(&map);
315 assert(mp_rad != 0x0) ;
316 for (
int k=0; k<np2; k++) {
317 for (
int j=0; j<nt; j++) {
318 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
320 ccA.set(k,j) = ddvr(k,j)*mp_rad->val_r_jk(0, 1., j, k) / double(l_q*(l_q+1)) ;
int j_max() const
Returns the larger time step j stored in *this.