LORENE
wave_utilities.C
1/*
2 * Miscellaneous functions for the wave equation
3 *
4 */
5
6/*
7 * Copyright (c) 2008 Jerome Novak
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
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 $" ;
27
28/*
29 * $Id: wave_utilities.C,v 1.11 2014/10/13 08:53:31 j_novak Exp $
30 * $Log: wave_utilities.C,v $
31 * Revision 1.11 2014/10/13 08:53:31 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.10 2008/12/05 13:09:10 j_novak
35 * Minor change in tilde_laplacian.
36 *
37 * Revision 1.9 2008/12/04 18:20:41 j_novak
38 * Better treatment of the case ETATZERO in BC initialization, and of dzpuis for
39 * evolution.
40 *
41 * Revision 1.8 2008/11/27 12:12:38 j_novak
42 * New function to initialize parameters for wave equation.
43 *
44 * Revision 1.7 2008/10/29 08:22:58 jl_cornou
45 * Compatibility conditions in the vector wave-equation case added
46 *
47 * Revision 1.6 2008/10/14 13:10:58 j_novak
48 * New function Dirichlet_BC_AtB, to compute Dirichlet boundary conditions on A and B potentials knowing them on the tensor h^{ij}.
49 *
50 * Revision 1.5 2008/08/27 08:11:47 j_novak
51 * Correction of a mistake in the index in evolve_outgoing_BC.
52 *
53 * Revision 1.4 2008/08/19 06:42:00 j_novak
54 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
55 * cast-type operations, and constant strings that must be defined as const char*
56 *
57 * Revision 1.3 2008/07/18 12:28:41 j_novak
58 * Corrected some mistakes.
59 *
60 * Revision 1.2 2008/07/18 09:17:35 j_novak
61 * New function tilde_laplacian().
62 *
63 * Revision 1.1 2008/07/11 13:20:54 j_novak
64 * Miscellaneous functions for the wave equation.
65 *
66 *
67 * $Header $
68 *
69 */
70
71#include"tensor.h"
72#include"evolution.h"
73
74namespace Lorene {
75void tilde_laplacian(const Scalar& B_in, Scalar& tilde_lap, int dl) {
76
77 if (B_in.get_etat() == ETATZERO) {
78 tilde_lap.set_etat_zero() ;
79 return ;
80 }
81 assert(B_in.check_dzpuis(0)) ; ;
82 if (dl == 0) {
83 tilde_lap = B_in.laplacian(0) ;
84 return ;
85 }
86 assert(B_in.get_etat() != ETATNONDEF) ;
87 const Map_af* map =dynamic_cast<const Map_af*>(&B_in.get_mp()) ;
88 assert(map != 0x0) ;
89
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() ;
99
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() ;
107 }
108 else {
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) ;
112 if (l_q > 1) {
113 l_q += dl ;
114 for (int i=0; i<mg.get_nr(lz); i++)
115 tilde_lap.set_spectral_va().c_cf->set(lz, k, j, i)
116 -= l_q*(l_q+1)*
117 (*B_over_r2.get_spectral_va().c_cf)(lz, k, j, i) ;
118 }
119 }
120 }
121 }
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 ;
125 }
126 tilde_lap.dec_dzpuis(2) ;
127 return ;
128}
129
130/* Performs one time-step integration of the wave equation, using
131 * a third-order Runge-Kutta scheme.
132 * phi = d fff / dt
133 * \Delta fff = d phi / dt
134 * Inputs are dt, fff, phi; outputs fnext, phinext.
135 */
136void runge_kutta3_wave_sys(double dt, const Scalar& fff, const Scalar& phi,
137 Scalar& fnext, Scalar& phinext, int dl) {
138
139 const Map& map = fff.get_mp() ;
140 Scalar k1 = phi ;
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 ;
147 Scalar k3 = dy2 ;
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. ;
151
152 return ;
153}
154
155void initialize_outgoing_BC(int nz_bound, const Scalar& phi, const Scalar& dphi,
156 Tbl& xij)
157{
158 Scalar source_xi = phi ;
159 source_xi.div_r_dzpuis(2) ;
160 source_xi += phi.dsdr() ;
161 source_xi.dec_dzpuis(2) ;
162 source_xi += dphi ;
163 if (source_xi.get_etat() == ETATZERO)
164 xij.annule_hard() ;
165 else {
166 source_xi.set_spectral_va().ylm() ;
167
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 ) ;
173
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) ;
178 xij.set(k, j)
179 = source_xi.get_spectral_va().c_cf->val_out_bound_jk(nz_bound, j, k) ;
180 if (l_q == 0)
181 xij.set(k,j) = 0 ;
182 }
183 }
184}
185
186
187/* Performs one time-step integration of the quantities needed for the
188 * enhanced outgoing-wave boundary condition. It DOES NOT impose the BC
189 * d phi / dr + d phi / dt + phi / r = xi(theta, varphi).
190 * nz_bound: index of the domain on which to impose the BC
191 * phi: the field that should leave the grid
192 * sphi: source of the Robin BC, without xi : a phi + b d phi / dr = sphi + xi
193 * ccc: (output) total source of the Robin BC
194 */
195void evolve_outgoing_BC(double dt, int nz_bound, const Scalar& phi, Scalar& sphi,
196 Tbl& xij, Tbl& xijm1, Tbl& ccc, int dl) {
197
198 const Map* map = &phi.get_mp() ;
199 const Map_af* mp_aff = dynamic_cast<const Map_af*>(map) ;
200 assert(mp_aff != 0x0) ;
201
202 const Mg3d& grid = *mp_aff->get_mg() ;
203#ifndef NDEBUG
204 int nz = grid.get_nzone() ;
205 assert(nz_bound < nz) ;
206 assert(phi.get_etat() != ETATZERO) ;
207 assert(sphi.get_etat() != ETATZERO) ;
208#endif
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) ;
220
221 double Rmax = mp_aff->get_alpha()[nz_bound] + mp_aff->get_beta()[nz_bound] ;
222
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) ;
234 if (l_q > 1) {
235 l_q += dl ;
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)
241 + souphi) / fact ;
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 ;
246 }
247 }
248
249}
250
251void Dirichlet_BC_AtB(const Evolution_std<Sym_tensor>& hb_evol,
252 const Evolution_std<Sym_tensor>& dhb_evol, Tbl& ccA, Tbl& ccB) {
253
254 int iter = hb_evol.j_max() ;
255 assert(dhb_evol.j_max() == iter) ;
256
257 Scalar mu_ddot = dhb_evol.time_derive(iter,3).mu() ;
258
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 ;
264 ccA.annule_hard() ;
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) ;
268 if (l_q>1)
269 ccA.set(k,j) = ddmu(k,j) / double(l_q*(l_q+1)-2) ;
270 }
271 }
272
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() ;
278
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) ;
285 if (l_q>1) {
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)) ;
289 }
290 }
291 }
292
293}
294
295
296void Dirichlet_BC_Amu(const Evolution_std<Vector>& vb_evol,
297 const Evolution_std<Vector>& dvb_evol, Tbl& ccA, Tbl& ccmu) {
298
299 int iter = vb_evol.j_max() ;
300 assert(dvb_evol.j_max() == iter) ;
301
302 Scalar vr_ddot = dvb_evol.time_derive(iter,3)(1) ;
303
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 ;
309 ccA.annule_hard() ;
310 ccmu.annule_hard() ;
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) ;
319 if (l_q>0) {
320 ccA.set(k,j) = ddvr(k,j)*mp_rad->val_r_jk(0, 1., j, k) / double(l_q*(l_q+1)) ;
321 }
322 }
323 }
324 }
325
326
327
328
329
330}
int j_max() const
Returns the larger time step j stored in *this.
Definition evolution.C:482
Lorene prototypes.
Definition app_hor.h:64