############################################################ # 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. # # 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: O(h^2) leap-frog with Kreiss/Oliger # dissipation, regularity at the origin and outgoing # radiation conditions. # 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 -1,0,1 float pi on g1 at -1,0,1 float a on g1 at -1,0,1 float l on g1 at -1,0,1 float tmr on g1 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] # O(h^2) centred time derivatives, without and with dissipation operator D0(f,t) := (<1>f[0] - <-1>f[0]) / (2 * dt) operator D0DISS(f,t) := (<1>f[0] - <-1>f[0] + epsdis / 16 * (6 *<-1>f[0] + <-1>f[2] + <-1>f[-2] - 4 * (<-1>f[1] + <-1>f[-1]))) / (2 * dt) # O(h^2) backward time derivative operator DB(f,t) := (3 * <1>f[0] - 4 * <0>f[0] + <-1>f[0]) / (2 * dt) # O(h^2) centred space derivative and O(h^2) r^3 derivative 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) # O(h^2) backward and forward space derivatives operator DB(f,r) := (3 * <1>f[0] - 4 * <1>f[-1] + <1>f[-2]) / (2 * dr) operator DF(f,r) := (-3 * <1>f[0] + 4 * <1>f[1] - <1>f[2]) / (2 * dr) # O(h^2) quadratic fit operator (for update of Pi(0,t)) operator QFIT(f,r) := +<1>f[0] - 4 * <1>f[1] / 3 + <1>f[2] / 3 # Discretized Klein-Gordon equation evaluate residual pp { [1:1] := <1>pp[0] = 0; [2:2] := D0(pp,t) = D0(l * pi / a ,r); [3:Nr-2] := D0DISS(pp,t) = D0(l * pi / a ,r); [Nr-1:Nr-1] := D0(pp,t) = D0(l * pi / a ,r); [Nr:Nr] := DB(pp,t) + (DB(pp,r) + pp / r) = 0; } ############################################################ # Note the ordering of residual evaluation in the following, # which ensures that the "quadratic fit" determining # the advanced central value uses properly updated (from # the interior eqns.) neighboring values. If the # [1:1] := QFIT(pi,r); statement were placed first, the # update scheme requires an extra iteration after scalar field # hits the origin (try it!) ############################################################ evaluate residual pi { [2:2] := D0(pi,t) = 3 * D0R3(r^2 * l * pp / a,r); [3:Nr-2] := D0DISS(pi,t) = 3 * D0R3(r^2 * l * pp / a,r); [Nr-1:Nr-1] := D0(pi,t) = 3 * D0R3(r^2 * l * pp / a,r); [1:1] := QFIT(pi,r); [Nr:Nr] := DB(pi,t) + (DB(pi,r) + pi / r) = 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,an,anm1], l[l,ln,lnm1], pp[pp,ppn,ppnm1], pi[pi,pin,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