added WIP (mostly works modulu bugs) 1d acoustics.

Author: alsam

examples/acoustics_1d_example1/ReadDataFile.jl

@@ -0,0 +1,96 @@
+module ReadDataFile
+
+export Reader, readstr, readbool, readint, readflt, readintarray, readfltarray, close_reader
+
+type Reader
+       
+    fname::String
+    istream::IOStream
+    line::String
+    lineno::Int
+
+    # a constructor: opens a file, sets lineno
+    function Reader(fname::String)
+        istream = open(fname, "r")
+        r = new(fname,istream, "", 0)
+        finalizer(r, close_reader)
+        return r
+    end
+end
+
+function readln(r::Reader)
+    const skip_rexp = r"^\s*(?:#|$)"
+    const rexp = r"^\s*(.*?)\s*(?:(#|=:)\s*(.*?)\s*$|$)"
+
+    while true
+        r.line = readline(r.istream)
+        r.lineno += 1
+        # skip empty lines and comments starting from #
+        if match(skip_rexp, r.line) == nothing
+            break
+        end
+    end
+    m = match(rexp, r.line)
+    r.line = m.captures[1]
+end
+
+function readstr(r::Reader)
+    readln(r)
+    return r.line
+end
+
+function readbool(r::Reader)
+    readln(r)
+    if r.line != "T" && r.line != "F"
+        local fname = r.fname, lineno = r.lineno
+        error("$fname($lineno): error parsing boolean")
+    end
+    return r.line == "T"
+end
+
+function readint(r::Reader)
+    readln(r)
+    return parse(Int,r.line)
+end
+
+function readflt(r::Reader)
+    readln(r)
+    return parse(Float64, r.line)
+end
+
+function readraw(r::Reader)
+    readln(r)
+    return split(r.line)
+end
+
+function readarray(r::Reader, parser, len)
+    buf = readraw(r)
+    buflen = length(buf)
+    @assert(buflen == len)
+    return [parser(buf[i]) for i = 1:buflen]
+end
+
+readintarray(r::Reader, len::Int) = readarray(r, x -> parse(Int, x), len)
+readfltarray(r::Reader, len::Int) = readarray(r, x -> parse(Float64, x), len)
+
+function close_reader(r::Reader)
+    Base.close(r.istream)
+end
+
+# ```jlcon
+# julia> include("ReadDataFile.jl")
+# 
+# julia> using ReadDataFile
+# 
+# julia> r = Reader("claw.data")
+# -I- claw.data 6 blank lines were skipped
+# Reader("claw.data",IOStream(<file claw.data>),"1",6,r"^\s*(.*?)\s*(?:#\s*(.*?)\s*$|$)")
+# 
+# julia> readint(r)
+# 1
+# 
+# julia> readflt(r)
+# -1.0
+# ```
+
+end # ReadDataFile

examples/acoustics_1d_example1/SetProb.jl

@@ -0,0 +1,62 @@
+
+module SetProb
+
+export
+    CParam, read_prob_data
+
+type CParam
+    rho   ::Float64
+    bulk  ::Float64
+    cc    ::Float64
+    zz    ::Float64
+
+    beta  ::Float64
+
+    # default constructor
+    function CParam()
+        new(0.0, 0.0, 0.0, 0.0, 0.0)
+    end
+end
+
+
+# Set the material parameters for the acoustic equations
+
+using ReadDataFile
+
+function read_prob_data(parms::CParam)
+
+    r = Reader("setprob.data")
+
+    # density:
+    parms.rho  = readflt(r)
+
+    # bulk modulus:
+    parms.bulk = readflt(r)
+
+    # sound speed:
+    parms.cc = sqrt(parms.bulk/parms.rho)
+
+    # impedance:
+    parms.zz = parms.cc*parms.rho
+
+    # beta for initial conditions:
+    parms.beta = readflt(r)
+
+    close_reader(r)
+end # read_prob_data
+
+end # SetProb
+
+using .SetProb
+using Base.Test.@test
+
+function test_read_data()
+    const parms = CParam()
+    read_prob_data(parms)
+    println("CParam = $parms")
+    return true
+end
+
+# @test test_read_data()
+
+

examples/acoustics_1d_example1/b4step1.jl

@@ -0,0 +1,15 @@
+function b4step1(mbc::Int, mx::Int, meqn::Int, q, xloweri::Float64, dxi::Float64, t::Float64, dt::Float64, mauxi::Int, aux)
+#
+#    ! Called before each call to step1.
+#    ! Use to set time-dependent aux arrays or perform other tasks.
+#    !
+#    ! This default version does nothing. 
+# 
+#    implicit none
+#    integer, intent(in) :: mbc,mx,meqn,maux
+#    real(kind=8), intent(in) :: xlower,dx,t,dt
+#    real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc)
+#    real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc)
+#
+end # subroutine b4step1
+

examples/acoustics_1d_example1/bc1.jl

@@ -0,0 +1,130 @@
+using OffsetArrays
+
+#
+#
+#=========================================================================#
+function bc1(meqn::Int, mbc::Int, mx::Int, xlower::Float64, dx::Float64,
+             q::OffsetArray{Float64,2}, maux::Int, aux::OffsetArray{Float64,2},
+             t::Float64, dt::Float64, mthbc::Array{Int,1})
+#=========================================================================#
+#
+#     # Standard boundary condition choices for claw2
+#
+#     # At each boundary  k = 1 (left),  2 (right):
+#     #   mthbc(k) =  0  for user-supplied BC's (must be inserted!)
+#     #            =  1  for zero-order extrapolation
+#     #            =  2  for periodic boundary coniditions
+#     #            =  3  for solid walls, assuming this can be implemented
+#     #                  by reflecting the data about the boundary and then
+#     #                  negating the 2'nd component of q.
+#     ------------------------------------------------
+#
+#     # Extend the data from the computational region
+#     #      i = 1, 2, ..., mx2
+#     # to the virtual cells outside the region, with
+#     #      i = 1-ibc  and   i = mx+ibc   for ibc=1,...,mbc
+#
+##      implicit double precision (a-h,o-z)
+##      dimension q(meqn,1-mbc:mx+mbc)
+##      dimension aux(maux,1-mbc:mx+mbc)
+##
+##      dimension mthbc(2)
+
+#
+#
+#-------------------------------------------------------
+#     # left boundary:
+#-------------------------------------------------------
+#      go to (100,110,120,130) mthbc(1)+1
+#
+    if mthbc[1] == 0
+#  100 continue
+
+        # user-specified boundary conditions go here in place of error output
+        println("*** ERROR *** mthbc(1)=0 and no BCs specified in bc1")
+
+    elseif mthbc[1] == 1
+
+#  110 continue
+        # zero-order extrapolation:
+        for ibc=1:mbc
+            for m=1:meqn
+                q[m,1-ibc] = q[m,1]
+            end
+        end
+#      go to 199
+
+    elseif mthbc[1] == 2 
+#  120 continue
+        # periodic:
+        for ibc=1:mbc
+            for m=1:meqn
+                q[m,1-ibc] = q[m,mx+1-ibc]
+            end
+        end
+#      go to 199
+
+    elseif mthbc[1] == 3 
+#  130 continue
+#     # solid wall (assumes 2'nd component is velocity or momentum in x):
+      for ibc=1:mbc
+         for m=1:meqn
+            q[m,1-ibc] = q[m,ibc]
+         end
+         q[2,1-ibc] = -q[2,ibc]
+      end
+#     # negate the normal velocity:
+#      go to 199
+
+#  199 continue
+    end # if mthbc[1] == ...
+
+#
+#-------------------------------------------------------
+#     # right boundary:
+#-------------------------------------------------------
+#      go to (200,210,220,230) mthbc(2)+1
+#
+    if mthbc[2] == 0
+#  200 continue
+        # user-specified boundary conditions go here in place of error output
+        println("*** ERROR *** mthbc(2)=0 and no BCs specified in bc2")
+        # go to 299
+
+    elseif mthbc[2] == 1
+#  210 continue
+        # zero-order extrapolation:
+        for ibc=1:mbc
+            for m=1:meqn
+                q[m,mx+ibc] = q[m,mx]
+            end
+        end
+#      go to 299
+
+    elseif mthbc[2] == 2 
+#  220 continue
+#     # periodic:  
+        for ibc=1:mbc
+            for m=1:meqn
+                q[m,mx+ibc] = q[m,ibc]
+            end
+        end
+#      go to 299
+
+    elseif mthbc[2] == 3 
+#  230 continue
+#     # solid wall (assumes 2'nd component is velocity or momentum in x):
+      for ibc=1:mbc
+         for m=1:meqn
+            q[m,mx+ibc] = q[m,mx+1-ibc]
+         end
+         q[2,mx+ibc] = -q[2,mx+1-ibc]
+      end
+#      go to 299
+
+    end # if mthbc[2] == ...
+#  299 continue
+#
+
+end
+

examples/acoustics_1d_example1/claw.data

@@ -0,0 +1,48 @@
+########################################################
+### DO NOT EDIT THIS FILE:  GENERATED AUTOMATICALLY ####
+### To modify data, edit  setrun.py                 ####
+###    and then "make .data"                        ####
+########################################################
+
+1                    =: num_dim             
+-1.0                 =: lower               
+1.0                  =: upper               
+400                  =: num_cells           
+
+2                    =: num_eqn             
+2                    =: num_waves           
+0                    =: num_aux             
+
+0.0                  =: t0                  
+
+1                    =: output_style        
+16                   =: num_output_times    
+0.8                  =: tfinal              
+T                    =: output_t0           
+
+1                    =: output_format       
+1 1                  =: iout_q              
+
+0.1                  =: dt_initial          
+1e+99                =: dt_max              
+1.0                  =: cfl_max             
+0.9                  =: cfl_desired         
+500                  =: steps_max           
+
+T                    =: dt_variable         
+2                    =: order               
+0                    =: verbosity           
+0                    =: source_split        
+0                    =: capa_index          
+F                    =: use_fwaves          
+
+4 4                  =: limiter             
+
+2                    =: num_ghost           
+1                    =: bc_lower            
+1                    =: bc_upper            
+
+F                    =: restart             
+'fort.q0006'         =: restart_file        
+0                    =: checkpt_style       
+