######################################################################## # # History: axi.8_rnpl (Hirschmann/Liebling) # # This code uses leap-frog evolution with some implicitness in # the outgoing boundary conditions. Dissipation in both the rho # and z directions a la Kreiss and Oliger is added to all three # basic dynamical variables. Long-time behaviour is considerably # improved over axi.8_rnpl, although, at least on coarse grids, # (level 0) scheme is still unstable with epsdis=0.5. # ######################################################################## # Pi_t = Omega_z + ( 2 d/d rho**2 )( rho Phi) # Phi_t = Pi_rho # Omega_t = Pi_z # # First order form: # Pi := phi_t # Phi := phi_rho # Omega := phi_z # # Copyright 1996,1997: Steve Liebling, Eric Hirschmann, Matthew Choptuik # UT Austin ######################################################################## ######################################################################## # Parameters, grid and grid functions ######################################################################## # This is how to set the memory size system parameter int memsiz := 4000000 parameter float rhomin parameter float rhomax parameter float zmin parameter float zmax parameter float amp parameter float delta parameter float r0 parameter float epsdis := 0 parameter float signum := 0 parameter float rho_center := 0 parameter float z_center := 0 parameter float a := 10 parameter float b := 10 rec coordinates t,rho,z uniform rec grid g1 [1:Nrho][1:Nz] {rhomin:rhomax} {zmin:zmax} float Phi on g1 at -1,0,1 float pi on g1 at -1,0,1 float Phirho on g1 at -1,0,1 float Phiz on g1 at -1,0,1 #These don't get evolved: float r on g1 at 0 attribute int out_gf encodeone := [ 1 0 0 0 0 0 0 0 0] ######################################################################## # Difference operators ######################################################################## #Leap Frog Difference Operators operator D_LF(f,t) := (<1>f[0][0] - <-1>f[0][0])/(2 * dt) operator D_LFDISSRHO(f,t) := (<1>f[0][0] - <-1>f[0][0] + epsdis/16 * (6 * <-1>f[0][0] + <-1>f[2][0] + <-1>f[-2][0] - 4 * (<-1>f[1][0] + <-1>f[-1][0])))/(2 * dt) operator D_LFDISSZ(f,t) := (<1>f[0][0] - <-1>f[0][0] + epsdis/16 * (6 * <-1>f[0][0] + <-1>f[0][2] + <-1>f[0][-2] - 4 * (<-1>f[0][1] + <-1>f[0][-1])))/(2 * dt) operator D_LFDISSRHOZ(f,t) := ( <1>f[0][0] - <-1>f[0][0] + epsdis/16 * (6 * <-1>f[0][0] + <-1>f[2][0] + <-1>f[-2][0] - 4 * (<-1>f[1][0] + <-1>f[-1][0])) + epsdis/16 * (6 * <-1>f[0][0] + <-1>f[0][2] + <-1>f[0][-2] - 4 * (<-1>f[0][1] + <-1>f[0][-1])) ) / (2 * dt) operator D_BW2(f,t) := ( 3 * <1>f[0][0] - 4 * <0>f[0][0] + <-1>f[0][0])/(2 * dt) operator D_LF(f,rho) := (<0>f[1][0] - <0>f[-1][0])/(2 * drho) operator D_LF(f,z) := (<0>f[0][1] - <0>f[0][-1])/(2 * dz) operator D_LFRET(f,rho) := (<-1>f[1][0] - <-1>f[-1][0])/(2 * drho) operator D_LFRET(f,z) := (<-1>f[0][1] - <-1>f[0][-1])/(2 * dz) operator D_FW2ADV(f,z ) := (-3 * <1>f[0][0] + 4 * <1>f[0][1] - <1>f[0][2])/(2 * dz) operator D_FW2ADV(f,rho) := (-3 * <1>f[0][0] + 4 * <1>f[1][0] - <1>f[2][0])/(2 * drho) operator D_FW2(f,rho) := (-3 * <0>f[0][0] + 4 * <0>f[1][0] - <0>f[2][0])/(2 * drho) operator D_FW2(f,z ) := (-3 * <0>f[0][0] + 4 * <0>f[0][1] - <0>f[0][2])/(2 * dz) operator D_FW2RET(f,rho) := (-3 * <-1>f[0][0] + 4 * <-1>f[1][0] - <-1>f[2][0])/(2 * drho) #Quadratic fit operator QFIT(f,rho) := (<1>f[0][0] - 4 * <1>f[1][0]/3 + <1>f[2][0]/3) operator QFITDISSZ(f,rho) := (<1>f[0][0] - 4 * <1>f[1][0]/3 + <1>f[2][0]/3 + epsdis/16 * (6 * <-1>f[0][0] + <-1>f[0][2] + <-1>f[0][-2] - 4 * (<-1>f[0][1] + <-1>f[0][-1]))) operator D_BW2ADV(f,rho) := ( 3 * <1>f[0][0] - 4 * <1>f[-1][0] + <1>f[-2][0])/(2 * drho) operator D_BW2ADV(f,z ) := ( 3 * <1>f[0][0] - 4 * <1>f[0][-1] + <1>f[0][-2])/(2 * dz) operator D_BW2(f,rho) := ( 3 * <0>f[0][0] - 4 * <0>f[-1][0] + <0>f[-2][0])/(2 * drho) operator D_BW2(f,z ) := ( 3 * <0>f[0][0] - 4 * <0>f[0][-1] + <0>f[0][-2])/(2 * dz) #Center Difference at Advanced time operator D_CADV(f,rho) := ( <1>f[1][0] - <1>f[-1][0])/(2 * drho) operator D_CADV(f,z ) := ( <1>f[0][1] - <1>f[0][-1])/(2 * dz) operator D_ADVAVG(f,z ) := ( -<1>f[0][0] + (<1>f[0][1] + <1>f[0][-1])/2 ) operator D_FEX(f,rho) := ( -<0>f[0][0] + 2 * <0>f[1][0] - <0>f[2][0] ) operator D_FEXRET(f,rho) := ( -<-1>f[0][0] + 2 * <-1>f[1][0] - <-1>f[2][0] ) operator D_rho_LF(f,rho) := ( rho[1]*<0>f[1][0] - rho[-1]*<0>f[-1][0] )/( rho[1]^2 - rho[-1]^2 ) ######################################################################## # Residual definitions (equations of motion) ######################################################################## evaluate residual pi { [2:2][2:2] := D_LF(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [3:Nrho-2][2:2] := D_LFDISSRHO(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [Nrho-1:Nrho-1][2:2] := D_LF(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [2:2][3:Nz-2] := D_LFDISSZ(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [3:Nrho-2][3:Nz-2] := D_LFDISSRHOZ(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [Nrho-1:Nrho-1][3:Nz-2] := D_LFDISSZ(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [2:2][Nz-1:Nz-1] := D_LF(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [3:Nrho-2][Nz-1:Nz-1] := D_LFDISSRHO(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [Nrho-1:Nrho-1][Nz-1:Nz-1] := D_LF(pi,t) - D_LF(Phiz,z) - 2.0*D_rho_LF(Phirho,rho); [2:Nrho-1][Nz:Nz] := r * D_BW2(pi,t) + rho * D_CADV(pi,rho) + z * D_BW2ADV(pi,z) + <1>Pi[0][0]; [2:Nrho-1][1:1] := r * D_BW2(pi,t) + rho * D_CADV(pi,rho) + z * D_FW2ADV(pi,z) + <1>Pi[0][0]; [Nrho:Nrho][2:Nz-1] := r * D_BW2(pi,t) + rho * D_BW2ADV(pi,rho) + z * D_CADV(pi,z) + <1>Pi[0][0]; [Nrho:Nrho][1:1] := r * D_BW2(pi,t) + rho * D_BW2ADV(pi,rho) + z * D_FW2ADV(pi,z) + <1>Pi[0][0]; [Nrho:Nrho][Nz:Nz] := r * D_BW2(pi,t) + rho * D_BW2ADV(pi,rho) + z * D_BW2ADV(pi,z) + <1>Pi[0][0]; [1:1][1:Nz] := QFIT(pi,rho); } evaluate residual Phirho { [2:Nrho-1][Nz-1:Nz] := D_LF(Phirho,t) - D_LF(Pi,rho); [2:Nrho-1][1:2] := D_LF(Phirho,t) - D_LF(Pi,rho); [3:Nrho-2][3:Nz-2] := D_LFDISSRHOZ(Phirho,t) - D_LF(Pi,rho); [2:2][3:Nz-2] := D_LF(Phirho,t) - D_LF(Pi,rho); [Nrho-1:Nrho-1][3:Nz-2] := D_LF(Phirho,t) - D_LF(Pi,rho); [1:1][1:Nz] := <1>Phirho[0][0]; [Nrho:Nrho][1:Nz] := D_BW2(Phirho,t) - D_BW2ADV(Pi,rho); } evaluate residual Phi { [1:Nrho][1:Nz] := D_LF(Phi,t) - <0>Pi[0][0]; } evaluate residual Phiz { [1:2][2:Nz-1] := D_LF(Phiz,t) - D_LF(Pi,z); [Nrho-1:Nrho][2:Nz-1] := D_LF(Phiz,t) - D_LF(Pi,z); [3:Nrho-2][2:Nz-1] := D_LFDISSRHOZ(Phiz,t) - D_LF(Pi,z); [3:Nrho-2][2:2] := D_LF(Phiz,t) - D_LF(Pi,z); [3:Nrho-2][Nz-1:Nz-1] := D_LF(Phiz,t) - D_LF(Pi,z); [1:Nrho][1:1] := D_BW2(Phiz,t) - D_FW2ADV(Pi,z); [1:Nrho][Nz:Nz] := D_BW2(Phiz,t) - D_BW2ADV(Pi,z); } ######################################################################## # Initializations and update structure ######################################################################## initialize r { [1:Nrho][1:Nz] := sqrt(rho^2 + z^2) } initialize Phi { [1:Nrho][1:Nz] := amp * exp(-((sqrt(((rho-rho_center)/a)^2 + ((z-z_center)/b)^2)-r0)/delta)^2) } initialize Pi { [2:Nrho-1][2:Nz-1] := signum * (<-1>Phi[0][0] + rho * D_LFRET(Phi,rho) + z * D_LFRET(Phi,z))/r; [1:1][1:Nz] := D_FEXRET(Pi,rho); [Nrho:Nrho][1:Nz] := 0 } initialize Phirho { [2:Nrho-1][1:Nz] := D_LFRET(Phi,rho); [1:1][1:Nz] := 0; [Nrho:Nrho][1:Nz] := 0 } initialize Phiz { [1:Nrho][2:Nz-1] := D_LFRET(Phi,z); [1:Nrho][1:1] := 0; [1:Nrho][Nz:Nz] := 0 } looper iterative auto update pi, Phi, Phirho, Phiz