############################################################ # History: emkglf_rnpl, ymcn_rnpl # # Spherically symmetric massless scalar collapse # (Einstein/massless-Klein-Gordon) # # See Choptuik PRL 70, 9-12, but note that the # slicing equation used here has had the # a_r / a term eliminated via the Hamiltonian constraint. # # Also note that this version uses Crank-Nicholson time-differencing. # # Equations of Motion: # # pp_t = (l pi / a)_r # pi_t = 3 (r^2 l pi / a)_(r^3) # a_r / a + (a^2 - 1) / (2 r) - 2 PI r (pp^2 + pi^2) = 0 # l_r / l - (a^2 - 1) / (2 r) - 2 PI r (pp^2 + pi^2) = 0 # # where # # phi = scalar field # a = radial metric function # l = lapse # pp = phi_r # pi = a phi_t / l # tmr = 2 m / r = (1 - a^{-2}) # # ph0 = storage for initial scalar field profile # f0, f2, g0 = storage for coefficient functions in # solution of constraint & slicing eqns. # # Difference scheme: # # Scalar field: Crank-Nicholson scheme, Kreiss/Oliger dissipation # (applied symmetrically), regularity at the origin # # Geometric quantities: O(h^2) "j+1/2 centred" differencing # with point-wise Newton iteration for Hamiltonian # constraint. Solvers are hand coded (See calc_al.inc). # # Update scheme: # # Scalar field quantites at level n+1 are updated first # from level n, n-1 values. Geometric variables at level # n+1 are then determined from level n+1 scalar field vars. # # Initial data: phi(r,0) is a Gaussian pulse, phi_t(r,0) # controllable with idsignum: # # idsignum = -1, 0, 1 -> ingoing, time-symm., outgoing ############################################################ # Set the memory size (only needed for Fortran) system parameter int memsiz := 1000000 # Domain extrema parameter float rmin := 0.0 parameter float rmax := 30.0 # Kreiss-Oilger dissipation coefficient parameter float epsdis := 0.5 # Gaussian initial data parameters parameter float amp := 1.0e-3 parameter float r0 := 15.0 parameter float del := 2.0 parameter float p := 2.0 parameter float idsignum := 1.0 # Hamiltonian constraint control and return code parameter int hcmxiter := 50 parameter float hcdlnatol := 1.0e-8 parameter int hcrc := 0 # Flat spacetime switch parameter int flatspace := 0 # Threshold for black hole formation parameter float tmrbh := 0.65 spherical coordinates t,r uniform spherical grid g1 [1:Nr] {rmin:rmax} # See header for grid function glossary float pp on g1 at 0,1 { out_gf := 0 } float pi on g1 at 0,1 { out_gf := 0 } float a on g1 at 0,1 { out_gf := 0 } float l on g1 at 0,1 { out_gf := 0 } float tmr on g1 { out_gf := 1 } float ph0 on g1 float f0 on g1 float f2 on g1 float g0 on g1 attribute int out_gf encodeone := [0 0 0 0 1 0 0 0 0 0 0] operator D_CN(f,t) := (<1>f[0] - <0>f[0]) / dt operator MU(f,t) := (<1>f[0] + <0>f[0]) / 2 operator DISS(f,t) := -epsdis / (16 * dt) * (6 *<0>f[0] + <0>f[2] + <0>f[-2] - 4 * (<0>f[1] + <0>f[-1])) operator D0(f,r) := (<0>f[1] - <0>f[-1]) / (2 * dr) operator D0R3(f,r) := (<0>f[1] - <0>f[-1]) / (r[1]^3 - r[-1]^3) operator D_BW(f,r) := (3*<0>f[0] - 4*<0>f[-1] + <0>f[-2]) / (2 * dr) operator D_FW(f,r) := (-3*<0>f[0] + 4*<0>f[1] - <0>f[2]) / (2 * dr) # O(h^2) quadratic fit operator (for update of Pi(0,t)) operator QFIT(f,r) := +<0>f[0] - 4 * <0>f[1] / 3 + <0>f[2] / 3 # Discretized Klein-Gordon equation evaluate residual pp { [1:1] := <1>pp[0] = 0; [2:2] := D_CN(pp,t) = MU(D0(l * pi / a ,r),t); [3:Nr-2] := D_CN(pp,t) = MU(DISS(pp,t),t) + MU(D0(l * pi / a ,r),t); [Nr-1:Nr-1] := D_CN(pp,t) = MU(D0(l * pi / a ,r),t); [Nr:Nr] := D_CN(pp,t) + MU(D_BW(pp,r) + pp / r,t) = 0; } evaluate residual pi { [1:1] := MU(QFIT(pi,r),t); [2:2] := D_CN(pi,t) = MU(3 * D0R3(r^2 * l * pp / a,r),t); [3:Nr-2] := D_CN(pi,t) = MU(DISS(pi,t),t) + MU(3 * D0R3(r^2 * l * pp / a,r),t); [Nr-1:Nr-1] := D_CN(pi,t) = MU(3 * D0R3(r^2 * l * pp / a,r),t); [Nr:Nr] := D_CN(pi,t) + MU(D_BW(pi,r) + pi / r,t) = 0; } # Gaussian initial data initialize ph0 {[1:Nr] := amp * exp(-abs((r - r0) / del)^p)} initialize pp {[1:1] := 0; [2:Nr-1] := D0(ph0,r) - ph0 / r; [Nr:Nr] := 0 } initialize pi {[1:1] := 0; [2:Nr-1] := idsignum * D0(ph0,r); [Nr:Nr] := 0 } looper iterative ############################################################ # pp and pi are updated using RNPL-generated code, ############################################################ auto update pp, pi ############################################################ # a, l and tmr are updated using code in include file # 'calc_al.inc' ############################################################ # RNPL produces (in 'updates.f') a routine called 'calc_al' # which includes the grid functions and parameters listed # after the keyword HEADER in the function header. The # construct # # a[a,an,anm1] # # causes RNPL to use the names 'a', 'an', and 'anm1' # as formal parameters in the subroutine, rather than # the RNPL defaults 'a_npl', 'a_n', 'a_nm1'. ############################################################ # Both here and in the manual INITIALIZE statement, RNPL # is a little unpredictable in the way that formal # arguments are ordered in the generated file: I can # virtually guarantee that they will *not* be in the # same order as listed in the RNPL source (i.e. here) # so BE CAREFUL since mismatched formal and actual # parameters is a very common error in Fortran, and # often difficult to detect. ############################################################ calc_al.inc calc_al UPDATES a, l, tmr HEADER a[a,anm1], l[l,lnm1], pp[pp,ppnm1], pi[pi,pinm1], tmr, f0, f2, g0, r, t, hcmxiter, hcdlnatol, hcrc, tmrbh, flatspace ############################################################ # ph0, pp, pi are initialized using RNPL code # ****** IMPORTANT!! At this time, if an RNPL code has # *any* manual initializations (as this one does), then # the functions which are to be auto-initialized *must* # be specified in an 'auto initialize' statement. ############################################################ auto initialize ph0, pp, pi ############################################################ # a, l and tmr are initialized using code in the include # file 'init_al.inc' ############################################################ # The syntax and semantics here is completely analagous to # the UPDATES statement, but note that RNPL is currently # a little quirky in that it wants you to initialize the # 'nm1'-level data (rather than the 'n'-level which would # be a little more natural. ############################################################ init_al.inc init_al INITIALIZE a, l, tmr HEADER tmr, a, l, pp, pi, f0, f2, g0, r, t, hcmxiter, hcdlnatol, hcrc, tmrbh, flatspace