Changeset 26492 for issm/trunk-jpl/src/jl/issm.jl

Timestamp:

10/20/21 10:01:31 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with archread, phew

File:

• issm/trunk-jpl/src/jl/issm.jl (modified) (4 diffs)

issm/trunk-jpl/src/jl/issm.jl

@@ -12 +12 @@
 function Mesh2dTriangle() #{{{
         return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int32}(undef, 0, 0))
+end# }}}
+function Base.show(io::IO, this::Mesh2dTriangle)# {{{
+
+        println(io,typeof(this),":")
+        for name in fieldnames(typeof(this))
+                a=getfield(this,name)
+                print(io,"   $(name) = ")
+                if !isempty(a)
+                        println(io, typeof(a), " of size ", size(a))
+                else
+                        println(io,"empty")
+                end
+        end
 end# }}}
 mutable struct Mesh3dPrism{T} <: AbstractMesh
@@ -41 +54 @@
         return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
 end# }}}
+mutable struct Initialization
+        vx::Vector{Float64}
+        vy::Vector{Float64}
+end
+function Initialization() #{{{
+        return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
+mutable struct Materials
+        rho_ice::Float64
+        rho_water::Float64
+        rho_freshwater::Float64
+        mu_water::Float64
+        heatcapacity::Float64
+        latentheat::Float64
+        thermalconductivity::Float64
+        temperateiceconductivity::Float64
+        effectiveconductivity_averaging::Int32
+        meltingpoint::Float64
+        beta::Float64
+        mixed_layer_capacity::Float64
+        thermal_exchange_velocity::Float64
+        rheology_B::Vector{Float64}
+        rheology_n::Vector{Float64}
+        rheology_law::String
+end
+function Materials() #{{{
+        return Materials(917., 1023., 1000., 0.001787, 2093., 3.34*10^5, 2.4, .24, 1, 273.15, 9.8*10^-8, 3974., 1.00*10^-4, Vector{Float64}(undef,0), Vector{Float64}(undef,0), "Cuffey")
+end# }}}
 
 mutable struct model
@@ -46 +87 @@
         geometry::Geometry
         mask::Mask
+        materials::Materials
+        initialization::Initialization
 end
 function model() #{{{
-        return model( Mesh2dTriangle(), Geometry(), Mask())
+        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -58 +101 @@
         @printf "%19s: %-22s -- %s\n" "geometry" typeof(md.geometry) "surface elevation, bedrock topography, ice thickness,..."
         @printf "%19s: %-22s -- %s\n" "mask" typeof(md.mask) "defines grounded and floating regions"
+        @printf "%19s: %-22s -- %s\n" "materials" typeof(md.materials) "material properties"
+        @printf "%19s: %-22s -- %s\n" "initialization" typeof(md.initialization) "initial state"
 
 end# }}}
 
-# Probably actually want something more like:
-# function solve(mesh::AbstractMesh, data::AbstractDataset)
-#     do_something_with_mesh(mesh) # This will dispatch
-#     do_something_with_data(data) # This will dispatch
-#     # do something
-# end
+#utils
+function issmdir() #{{{
+        issmdir = ENV["ISSM_DIR"]
+
+        if isempty(issmdir)
+                error("Could not determine the location of ISSM")
+        else
+                return issmdir
+        end
+end#}}}
+function archread(filename::String,variablename::String) #{{{
+
+        #initialize variables
+        found = false
+
+        #open file
+        output = open(filename, "r") do f
+
+                while !eof(f)
+                        reclen  = bswap(read(f, Int32))
+                        rectype = bswap(read(f, Int32))
+                        if rectype!=1
+                                error("Expected variable of type string")
+                        else
+                                fieldname_length = bswap(read(f, Int32))
+                                field_name = String(read(f, fieldname_length))
+                        end
+                        rec_length = bswap(read(f, Int32))
+                        field_type = bswap(read(f, Int32))
+                        if field_type==2
+                                data = bswap(read(f, Float64))
+                        elseif field_type==3
+                                rows = bswap(read(f, Int32))
+                                cols = bswap(read(f, Int32))
+                                data = reinterpret(Float64, read(f, sizeof(Float64)*rows*cols))
+                                data .= ntoh.(data)
+                                data = reshape(data, (rows,cols))
+                        else
+                                error("Error: Encountered invalid field type when reading data.")
+                        end
+
+                        if field_name == variablename
+                                found = true
+                                return data
+                        end
+                end
+        end
+
+        return output
+end# }}}