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 
26 char 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 
74 namespace Lorene {
75 void 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  */
136 void 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 
155 void 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  */
195 void 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 
251 void 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 
296 void 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 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64