added PROGRAM1 - translated from Fortran.

Author: alsam

examples/scalar_law/PROGRAM1/consdiff.jl

@@ -0,0 +1,55 @@
+#***********************************************************************
+#  Copyright 2006 John A. Trangenstein
+#
+#  This software is made available for research and instructional use 
+#  only. 
+#  You may copy and use this software without charge for these 
+#  non-commercial purposes, provided that the copyright notice and 
+#  associated text is reproduced on all copies.  
+#  For all other uses (including distribution of modified versions), 
+#  please contact the author at
+#    John A. Trangenstein
+#    Department of Mathematics
+#    Duke University
+#    Durham, NC 27708-0320
+#    USA
+#  or
+#    [email protected]
+#  
+#  This software is made available "as is" without any assurance that it
+#  is completely correct, or that it will work for your purposes.  
+#  Use the software at your own risk.
+#***********************************************************************
+
+function consdiff(fic::Int,lac::Int,fif::Int,laf::Int,fix::Int,lax::Int,
+                  ifirst::Int,ilast::Int,
+                  x::OffsetArray{Float64,1},flux::OffsetArray{Float64,1},
+                  conservd::OffsetArray{Float64,1})
+#   ******************************************************************
+#   update conservd to new time
+#   ******************************************************************
+    @unsafe for ic=ifirst:ilast
+        conservd[ic] -= (flux[ic+1]-flux[ic]) / (x[ic+1]-x[ic])
+    end
+end
+
+function stabledt(fc::Int,lc::Int,fm::Int,lm::Int,ifirst::Int,ilast::Int,
+                  conserved::OffsetArray{Float64,1},
+                  lambda_cell::OffsetArray{Float64,1},
+                  x::OffsetArray{Float64,1})
+#   ******************************************************************
+#   compute stable timestep
+#   ******************************************************************
+    sdt=huge
+    @unsafe for ic=ifirst:ilast
+        abs_lambda=abs(lambda_cell[ic])
+        dx = x[ic+1]-x[ic]
+        if abs_lambda > roundoff*dx
+            sdt=min(sdt,dx/abs_lambda)
+        end
+    end
+    @assert sdt > 0.0
+    sdt
+end
+
+

examples/scalar_law/PROGRAM1/linaddmain.jl

@@ -0,0 +1,102 @@
+#***********************************************************************
+#  Copyright 2006 John A. Trangenstein
+#
+#  This software is made available for research and instructional use 
+#  only. 
+#  You may copy and use this software without charge for these 
+#  non-commercial purposes, provided that the copyright notice and 
+#  associated text is reproduced on all copies.  
+#  For all other uses (including distribution of modified versions), 
+#  please contact the author at
+#    John A. Trangenstein
+#    Department of Mathematics
+#    Duke University
+#    Durham, NC 27708-0320
+#    USA
+#  or
+#    [email protected]
+#  
+#  This software is made available "as is" without any assurance that it
+#  is completely correct, or that it will work for your purposes.  
+#  Use the software at your own risk.
+#***********************************************************************
+
+using OffsetArrays
+
+const roundoff  =  1.0-14
+const small     =  1.0-20
+const huge      =  1.0e300
+const undefind  =  1.0e300
+const lnrndoff  =  14.0e0
+const lnsmall   =  20.0e0
+
+# problem-specific parameters:
+const jump      =  0.0
+const x_left    = -0.2
+const x_right   =  1.0
+const statelft  =  2.0
+const statergt  =  0.0
+const velocity  =  1.0
+
+include("riemprob.jl")
+include("linearad.jl")
+include("consdiff.jl")
+include("upwind.jl")
+
+function main()
+    const ncells = 10000
+
+    u    = OffsetArray(Float64, -2:ncells+1)
+    x    = OffsetArray(Float64,  0:ncells)
+    flux = OffsetArray(Float64,  0:ncells)
+    dfdu = OffsetArray(Float64, -2:ncells+1)
+
+    const nsteps   =  100
+    #const nsteps   =  10000
+    const tmax     =  0.8
+    const cfl      =  0.9
+
+    # array bounds:
+    const fc=-2
+    const lc=ncells+1
+    const fm=0
+    const lm=ncells-1
+    const fs=0
+    const ls=ncells-1
+    const ifirst=0
+    const ilast=ncells-1
+
+    initsl(ncells,fc,lc,fm,lm,ifirst,ilast, u,x)
+
+    #println("x : $x u : $u")
+    bcmesh(fm,lm,ncells, x)
+    bccells(fc,lc,ncells, u)
+    fluxderv(fc,lc,fc,lc, u, dfdu)
+
+    istep=0
+    t=0.0
+    while istep < nsteps && t < tmax
+        bccells(fc,lc,ncells, u)
+        fluxderv(fc,lc,fc,lc, u, dfdu)
+        dt=cfl*stabledt(fc,lc,fm,lm,ifirst,ilast, u, dfdu,x)
+        method(dt,fc,lc,fm,lm,fs,ls,ifirst,ilast, u, flux)
+        consdiff(fc,lc,fs,ls,fm,lm,ifirst,ilast, x,flux, u)
+
+        t=t+dt
+        istep=istep+1
+    end
+
+    # write final results (plot later)
+    #@unsafe for ic=0:ncells-1
+    for ic=0:ncells-1
+        xc = (x[ic]+x[ic+1])*0.5
+        uc = u[ic]
+        @printf("%e %e\n",xc,uc)
+    end
+
+
+end # main
+
+main()
+
+

examples/scalar_law/PROGRAM1/linearad.jl

@@ -0,0 +1,46 @@
+#***********************************************************************
+#  Copyright 2006 John A. Trangenstein
+#
+#  This software is made available for research and instructional use 
+#  only. 
+#  You may copy and use this software without charge for these 
+#  non-commercial purposes, provided that the copyright notice and 
+#  associated text is reproduced on all copies.  
+#  For all other uses (including distribution of modified versions), 
+#  please contact the author at
+#    John A. Trangenstein
+#    Department of Mathematics
+#    Duke University
+#    Durham, NC 27708-0320
+#    USA
+#  or
+#    [email protected]
+#  
+#  This software is made available "as is" without any assurance that it
+#  is completely correct, or that it will work for your purposes.  
+#  Use the software at your own risk.
+#***********************************************************************
+
+function fluxderv(fc::Int,lc::Int,ifirst::Int,ilast::Int,
+                  conservd::OffsetArray{Float64,1}, dfdu::OffsetArray{Float64,1})
+#   ******************************************************************
+#   compute derivative of flux with respect to conserved variable
+#   ******************************************************************
+    @unsafe for ic=ifirst:ilast
+        dfdu[ic]=velocity
+    end
+end
+
+function riemann(fs::Int,ls::Int,ifirst::Int,ilast::Int,
+                 left::OffsetArray{Float64,1},right::OffsetArray{Float64,1},
+                 flux::OffsetArray{Float64,1})
+#   ******************************************************************
+#   the following riemann solver is hard-wired for linear advection
+#   with velocity > 0
+#   ******************************************************************
+    @unsafe for ie=ifirst:ilast+1
+        flux[ie]=velocity*left[ie]
+    end
+end 
+
+

examples/scalar_law/PROGRAM1/riemprob.jl

@@ -0,0 +1,91 @@
+#***********************************************************************
+#  Copyright 2006 John A. Trangenstein
+#
+#  This software is made available for research and instructional use 
+#  only. 
+#  You may copy and use this software without charge for these 
+#  non-commercial purposes, provided that the copyright notice and 
+#  associated text is reproduced on all copies.  
+#  For all other uses (including distribution of modified versions), 
+#  please contact the author at
+#    John A. Trangenstein
+#    Department of Mathematics
+#    Duke University
+#    Durham, NC 27708-0320
+#    USA
+#  or
+#    [email protected]
+#  
+#  This software is made available "as is" without any assurance that it
+#  is completely correct, or that it will work for your purposes.  
+#  Use the software at your own risk.
+#***********************************************************************
+
+function initsl(ncells::Int, fc::Int, lc::Int, fx::Int, lx::Int, ifirst::Int, ilast::Int,
+                conservd::OffsetArray{Float64,1}, x::OffsetArray{Float64,1})
+#   ******************************************************************
+#   interior values only; others defined by boundary conditions
+#   ******************************************************************
+    dx=(x_right-x_left)/ncells
+    @unsafe for ie=ifirst:ilast+1
+        x[ie]=x_left+ie*dx
+    end
+    ijump=max(ifirst-1,min(convert(Int,round(ncells*(jump-x_left)/(x_right-x_left))),ilast+1))
+    @unsafe for ic=ifirst:ijump-1
+        conservd[ic]=statelft
+    end
+    frac=(jump-x_left-ijump*dx)/(x_right-x_left)
+    conservd[ijump]=statelft*frac+statergt*(1.0-frac)
+    @unsafe for ic=ijump+1:ilast
+        conservd[ic]=statergt
+    end
+end
+
+function bcmesh(fim::Int,lam::Int,ncells::Int, x::OffsetArray{Float64,1})
+#   ******************************************************************
+#   right side 
+#   ******************************************************************
+    if lam >= ncells
+        it=ncells
+        dx=x(it)-x(it-1)
+        @unsafe for ie=it+1:lam+1
+            x[ie]=x[it]+dx*(ie-it)
+        end
+    end
+#   ******************************************************************
+#   left side 
+#   ******************************************************************
+    if fim < 0
+        dx=x(1)-x(0)
+        @unsafe for ie=fim:-1
+            x[ie]=x[0]+dx*(ie)
+        end
+    end      
+end
+
+function bccells(fic::Int,lac::Int,ncells::Int, conservd::OffsetArray{Float64,1})
+#     ******************************************************************
+#     right side 
+#     ******************************************************************
+    if lac >= ncells
+        # println("treating right side")
+        @unsafe for ic=ncells:lac
+            # outgoing wave:
+            conservd[ic]=conservd[ncells-1]
+        end  
+    end
+#     ******************************************************************
+#     left side 
+#     ******************************************************************
+    if fic < 0
+        # println("treating left side")
+        @unsafe for ic=fic:-1
+            # outgoing wave:
+            # conservd[ic]=conservd[0]
+
+            # specified value
+            conservd[ic]=statelft
+        end  
+    end
+end
+

examples/scalar_law/PROGRAM1/upwind.jl

@@ -0,0 +1,35 @@
+#***********************************************************************
+#  Copyright 2006 John A. Trangenstein
+#
+#  This software is made available for research and instructional use 
+#  only. 
+#  You may copy and use this software without charge for these 
+#  non-commercial purposes, provided that the copyright notice and 
+#  associated text is reproduced on all copies.  
+#  For all other uses (including distribution of modified versions), 
+#  please contact the author at
+#    John A. Trangenstein
+#    Department of Mathematics
+#    Duke University
+#    Durham, NC 27708-0320
+#    USA
+#  or
+#    [email protected]
+#  
+#  This software is made available "as is" without any assurance that it
+#  is completely correct, or that it will work for your purposes.  
+#  Use the software at your own risk.
+#***********************************************************************
+
+function method(dt::Float64,fc::Int,lc::Int,fm::Int,lm::Int,fs::Int,ls::Int,ifirst::Int,ilast::Int,
+                conserved::OffsetArray{Float64,1},
+                flux::OffsetArray{Float64,1})
+#   ******************************************************************
+#   multiply fluxes by dt (ie, compute time integral over cell side)
+#   ******************************************************************
+    vdt=velocity*dt
+    @unsafe for ie=ifirst:ilast+1
+        flux[ie]=vdt*conserved[ie-1]
+    end
+end
+