Download Computer Modeling Fluid Dynamics - Problem Set 2 | MPO 662 and more Assignments Geology in PDF only on Docsity!
MPO 662 – Problem Set 2
1 Leap-Frog Trapezoidal time stepping
One remedy to the computational mode of the Leap-Frog scheme is to combine it with the trapezoidal step. The combined scheme can be written as follows:
u∗^ = un−^1 + 2∆tf (un, tn) (1)
un+1^ = un^ + ∆tf
( un^ + u∗ 2
, tn+ (^12)
) (2)
The scheme can be viewed as a predictor corrector method. The first step predicts u at time level n + 1 using a leap-frog scheme, while the second one is a trapezoidal step centered at n + 12. Unlike the pure trapezoidal scheme, this one is explicit and uses the average of un^ and
u∗^ as an estimate of un+^
(^12)
. Investigate the order of the scheme for the case f = κu. Derive the expression for the amplification factor. Is there a parasitic mode? and if yes which one? Plot the amplification factor for κ∆t = z = x + iy where |x| < 1 and 0 ≤ y ≤ 2. Compare its stability region to that of the Leap-Frog, RK2 and AB2 schemes.
Hints Plotting stability diagrams is rather cumbersome especially if they involve the square roots of complex numbers. Here is a sample matlab script to achieve that.
M=201; xmin =-1; xmax=1.; dx=(xmax-xmin)/(M-1); x =xmin:dx:xmax; N=201; ymin = 0; ymax=2.; dy=(ymax-ymin)/(N-1); y =ymin:dy:ymax; z= x’ones(size(y)) + i * ones(size(x))’y; % roots of equation a A^2 + bA +c, where a,b and c are functions of z % a=; % b=; % c=; det = sqrt(b.^2-4a.c); A1 = 0.5(-b+det)./(2a); A2 = 0.5(-b-det)./(2*a); A = max(abs(A1),abs(A2)); contour(x,y,A’,[1 1]); % will draw the contour of neutral stability
2 Inertical oscillations
Inertial oscillations on the sphere are described by the following pair of equations.
d dt
( u v
)
( 0 f −f 0
) ( u v
) (3)
- Writing the equations as a system of the form dq/dt = P q, with q =
( u v
) , find the eigenvalues and eigenvectors of the matrix P , where
P =
( 0 f −f 0
)
- Use the eigenvectors to transform the system to an uncoupled set of equations by intro- ducing the new variables u′^ and v′: ( u v
) = R
( u′ v′
) (4)
where R is the matrix whose columns are the eigenvectors of P.
- Based on these eigenvalues which ODE integration schemes are suitable for the numerical solution.
3 RK
Write a program to integrate equations 3 using a TVDRK3 (Gottlieb and Shu, 1998), RK4, and AB3 schemes (see notes for details on the RK4 scheme). The initial conditions are u = 1 and v = 0. Use f = 1. Integrate until time T = 40 using ∆t=4,1.6,0.8,0.4,0.2, and 0.1. Monitor the evolution of u^2 + v^2. Compare the solution at T=40 with the exact solution ue(t) = cos t
and ve(t) = − sin t. Plot the error ǫ =
√ (u − ue)^2 + (v − ve)^2 at T = 40 as a function of ∆t using a log-log scale and estimate the slope of the curve. The TVD-RK3 scheme is a third order Runge-Kutta time-stepping scheme with 3 stages. Its intermediate stages can be seen as successive corrections to un+1^ using a convex combination of weights. In brief the scheme can be described algorithmically as follows:
u(1)^ = un^ + ∆tf (un) (5)
u(2)^ =
un^ +
u(1)^ +
∆tf (u(1)) (6)
un+1^ =
un^ +
u(2)^ +
∆tf (u(2)) (7) (8)
The program should be divided into section. The first one initializes the variables and sets the numerical and physical parameters (like f , ∆t and the number of time steps). The second part consists of the time-loop in which you call a subroutine to advance the solution from time level n to time level n + 1. Within the loop you should design an RK4 subroutine that takes care of the four stages of the RK4 time step. Each of these stages requires its own evaluation of the right hand side P q. A good program would then isolate the computations
program ode use timeint! use time-integration module use rhsmod, only : Diagnostics! use inertial module for diagostics implicit none integer, parameter :: nn=2! 2 dependent variables integer, parameter :: ndiag=2! print diagnostics every ndiag steps real :: q(nn)! q(1) is u and q(2) is v. real :: dt! timestep real :: ntimestep! number of timesteps to take
write(6,)’Enter dt, ntimestep’ read(5,) dt, ntimestep! read parameters call Initialize(q,nn)! initialize q do n = 1,ntimestep! time-loop time = (n-1)dt call RK4(q,time,nn)! do a single RK4 step ! call RK3(q,time,nn)! do a single RK3 step ! call AB3(q,time,nn)! do a single AB3 step if (mod(n,isnap)==0) then time = ndt call Diagnostics(q,err,KE,time,nn) endif enddo stop end program ode module timeint use rhsmod contains subroutine RK4(q,time,nn) implicit none integer, intent(in) :: nn real, intent(in) :: time real, intent(inout) :: q(nn)! on input q is at time n, and on output at time n+ real :: r(nn) call RHS(r,q,time,nn)! calculate the right hand side . return end subroutine RK subroutine RK3(q,time,nn) implicit none integer, intent(in) :: nn real, intent(in) :: time real, intent(inout) :: q(nn)! on input q is at time n, and on output at time n+ . end subroutine RK
end module timeint module rhsmod implicit none contains subroutine RHS(r,q,time,nn)! compute right hand side integer, intent(in) :: nn real, intent(in) :: time real, intent(in) :: q(nn)! function values needed to compute rhs real, intent(in) :: r(nn)! right hand sides . return end subroutine RHS subroutine Diagnostics(q,err,KE,time,nn)! diagnostics . end subroutine Diagnostics subroutine Initialize(q,nn)! Set the initial conditions . end subroutine Initialize end module rhsmod
