M578 - Alexiades
                  Lab 7
      Super-Time-Stepping for Diffusion Equation

Implement the Super-Time-Stepping scheme (STS) on (a copy of...)
your explicit code from Lab 3. 
See suggestion below.

Debug using Nsts=1 , ν=0.0 (the explicit scheme itself), 
and  Nsts=2, ν=0.01.

In the following, use M=64, D=0.1, and compare with the exact solution 
to find the max error at time tmax=100.
Along with the STS parameters Nsts and ν, your code should
count and output the number of timesteps and number of supersteps. 
(code runs instantaneously, so use number of timessteps to measure efficiency)

1. Run the explicit scheme itself (Nsts=1 , v=0.0) to record the
   number of timesteps and max error.

2. Choose Nsts = 10 and ν = 0.0,  0.01,  0.04, 0.05,  0.1,  0.2
   and compare the errors (against the exact solution at tmax).
   Which one is most efficient (least total number of timesteps) 
   with error up to 0.001 ?

3. Try Nsts = 20 and various ν, to see which combination performs
   most efficiently within an error of up to 0.001.

4. Also try Nsts = 50 and various ν.  

For each of 2., 3., 4. make an appropriate table of what ν you tried,
and for the most efficient case  compute Speedup as the ratio 
of nsteps for Euler over nsteps for STS:  
       Speedup = nstepsEuler / nstepsSTS 
Turn in:
  A. Table and answer to 2.
  B. Table and answer to 3.
  C. Table and answer to 4.

          STS implementation
Computation of the sub-steps   τi , i=1,...,Nsts   can be done in various ways.
The following routine constructs large-to-small τi (or small-to-large by uncommenting a line).
 subroutine SUPER_TIME_STEP( Nsts, nu, dtEXPL, TAU, durSTS )
            Pi2 = 2*ATAN(1.d0)
         durSTS = 0.d0		!!...initialize duration of superstep
    DO i = 1, Nsts
      !!................. large-to-small ...............: 
           supi = ( 2*i − 1.d0 ) / Nsts
      !!................. small-to-large (uncomment the next line):
        !  supi = 2.d0 − supi
          super = (nu − 1.d0)*COS(supi*Pi2) + nu + 1.d0
         TAU(i) = dtEXPL / super
         durSTS = durSTS + TAU(i)	!!...duration of a superstep 
        print*,' Nsts, nu:',Nsts, nu, ' ==> STS_duration =',durSTS
This needs to be called once to construct the array TAU(i).
Then, in main, instead of using fixed timesteps dt=dtEXPL, use an STS loop:
      nSupSteps = nSupSteps + 1		!!...count supersteps...
    DO ii = 1, Nsts       !!...STS loop...
             dt = TAU(ii)
           time = time + dt
         nsteps = nsteps + 1		!!...count timesteps...
          CALL FLUX( ... )
          CALL PDE ( ... )
    ENDDO                 !!...end of STS loop
    IF( time >= tout ) THEN		!!...now can output
         CALL COMPARE( ... )
         CALL OUTPUT ( ... )
         tout = tout + dtout
Note that only dtEXPL is used in STS, should set factor=1.
Also, comparison is called for at tmax only, can set dtout=tmax.