############################################################ # RNPL program to solve the 1d wave equation # # phi_tt = phi_xx # # recast in first order form # # pp_t = pi_x # pi_t = pp_x # # where # # pp := phi_x # pi := phi_t # # The program uses Crank-Nicholson differencing with # (implicit )ingoing/outgoing radiation conditions, using # O(h^2) forward and backwards differences, at the boundaries. # # Initial data is a Gaussian for pp with pi = idsignum * pp, # where # # idsignum = -1, 0, 1 -> ingoing, time-symmetric, outgoing # initial data, respectively ############################################################ #----------------------------------------------------------- # Definition of memory size (only needed for Fortran) #----------------------------------------------------------- system parameter int memsiz := 1000000 #----------------------------------------------------------- # Definition of parameters and associated default values. #----------------------------------------------------------- #----------------------------------------------------------- # Note that "xmin" and "xmax" are special in that they are # also implicitly declared to be parameters via the # definition of the coordinate system below. #----------------------------------------------------------- parameter float xmin := 0.0 parameter float xmax := 1.0 #----------------------------------------------------------- # The following four parameters are used in the # specification of the initial data. #----------------------------------------------------------- parameter float amp := 1.0 parameter float xc := 0.5 parameter float xwid := 0.05 parameter float idsignum := 0.0 #----------------------------------------------------------- # Definition of coordinate system, note that the first # coordinate is always assumed to be the time coordinate. #----------------------------------------------------------- rect coordinates t, x #----------------------------------------------------------- # Definition of finite-difference grid: [1:Nx] specifies # the index range, {xmin:xmax} the coordinate range. #----------------------------------------------------------- uniform rect grid g1 [1:Nx] {xmin:xmax} #----------------------------------------------------------- # Definition of grid functions: since a Crank-Nicholson # scheme is being used, the grid functions are defined at # "temporal offsets" 0 and 1, corresponding to "current" and # "advanced" time levels. The directive {out_gf = 1} enables # default output of the grid function (output can be disabled # via {out_gf = 0}, by omitting the directive completely, or # via modification of the file .rnpl.attributes prior to # program invocation). #----------------------------------------------------------- float pp on g1 at 0,1 {out_gf = 1} float pi on g1 at 0,1 {out_gf = 1} #----------------------------------------------------------- # FINITE DIFFERENCE OPERATOR DEFINITIONS #----------------------------------------------------------- #----------------------------------------------------------- # Crank Nicholson time derivative operator (first forward # difference) #----------------------------------------------------------- operator DCN(f,t) := (<1>f[0] - <0>f[0]) / dt #----------------------------------------------------------- # Forward time averaging operator #----------------------------------------------------------- operator MU(f,t) := (<1>f[0] + <0>f[0]) / 2 #----------------------------------------------------------- # O(h^2) centred spatial derivative operator #----------------------------------------------------------- operator D0(f,x) := (<0>f[1] - <0>f[-1]) / (2*dx) #----------------------------------------------------------- # O(h^2) backwards spatial derivative operator #----------------------------------------------------------- operator DB(f,x) := (3*<0>f[0] - 4*<0>f[-1] + <0>f[-2]) / (2*dx) #----------------------------------------------------------- # O(h^2) forwards spatial derivative operator #----------------------------------------------------------- operator DF(f,x) := (-3*<0>f[0] + 4*<0>f[1] - <0>f[2]) / (2*dx) #----------------------------------------------------------- # DIFFERENCE EQUATION DEFINITIONS #----------------------------------------------------------- evaluate residual pp { [1:1] := DCN(pp,t) = MU(DF(pp,x),t); [2:Nx-1] := DCN(pp,t) = MU(D0(pi,x),t); [Nx:Nx] := DCN(pp,t) = MU(-DB(pp,x),t); } evaluate residual pi { [1:1] := DCN(pi,t) = MU(DF(pi,x),t); [2:Nx-1] := DCN(pi,t) = MU(D0(pp,x),t); [Nx:Nx] := DCN(pi,t) = MU(-DB(pi,x),t); } #----------------------------------------------------------- # INITIALIZATION STATEMENTS #----------------------------------------------------------- #----------------------------------------------------------- # Intialize to an ingoing,outgoing or time-symmetric # gaussian pulse (in the spatial derivative of the scalar # field), dependent on the value of idsignum. #----------------------------------------------------------- initialize pp { [1:Nx] := amp*exp(-(x-xc)^2/xwid^2) } initialize pi { [1:Nx] := idsignum * amp*exp(-(x-xc)^2/xwid^2) } auto initialize pp, pi #----------------------------------------------------------- # Definition of type of time stepping algorithm. The use # of "iterative" here, combined with the "auto update ..." # statement below, results in a scheme whereby the # residuals defined above are iteratively relaxed using # a point-wise Newton-Gauss-Seidel technique, until the # residual norms are below a certain threshold. #----------------------------------------------------------- looper iterative #----------------------------------------------------------- # The following statement directs RNPL to automatically # generate code to update the grid functions using the # residual definitions. #----------------------------------------------------------- auto update pp, pi