Changeset 26522


Ignore:
Timestamp:
10/29/21 08:05:53 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: changing the structure and done with ModelProcessor

Location:
issm/trunk-jpl/src/jl
Files:
2 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r26518 r26522  
    1 using Printf
     1# Julia version of the Ice-sheet and Sea-level System Model
     2#
     3# Author: Mathieu Morlighem
     4# email:  mathieu.morlighem@dartmouth.edu
    25
    3 #Define classes used by model struct
    4 abstract type AbstractMesh end
    5 mutable struct Mesh2dTriangle <: AbstractMesh
    6         numberofvertices::Int64
    7         numberofelements::Int64
    8         x::Vector{Float64}
    9         y::Vector{Float64}
    10         elements::Matrix{Int64}
     6module ISSM
     7
     8include("classes.jl")
     9include("utils.jl")
     10include("triangle.jl")
     11include("parameterization.jl")
     12include("solve.jl")
     13#include("plotmodel.jl")
     14
    1115end
    12 function Mesh2dTriangle() #{{{
    13         return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(undef, 0, 0))
    14 end# }}}
    15 function Base.show(io::IO, this::Mesh2dTriangle)# {{{
    16 
    17         println(io,typeof(this),":")
    18         for name in fieldnames(typeof(this))
    19                 a=getfield(this,name)
    20                 print(io,"   $(name) = ")
    21                 if !isempty(a)
    22                         println(io, typeof(a), " of size ", size(a))
    23                 else
    24                         println(io,"empty")
    25                 end
    26         end
    27 end# }}}
    28 mutable struct Mesh3dPrism{T} <: AbstractMesh
    29         numberofvertices::Int64
    30         numberofelements::Int64
    31         numberoflayers::Int64
    32         x::Vector{Float64}
    33         y::Vector{Float64}
    34         z::Vector{Float64}
    35         elements::Matrix{Int64}
    36 end
    37 function Mesh3dPrism() #{{{
    38         return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0))
    39 end# }}}
    40 mutable struct Geometry
    41         surface::Vector{Float64}
    42         base::Vector{Float64}
    43         thickness::Vector{Float64}
    44         bed::Vector{Float64}
    45 end
    46 function Geometry() #{{{
    47         return Geometry( Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0))
    48 end# }}}
    49 mutable struct Mask
    50         ocean_levelset::Vector{Float64}
    51         ice_levelset::Vector{Float64}
    52 end
    53 function Mask() #{{{
    54         return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
    55 end# }}}
    56 mutable struct Initialization
    57         vx::Vector{Float64}
    58         vy::Vector{Float64}
    59 end
    60 function Initialization() #{{{
    61         return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
    62 end# }}}
    63 mutable struct Stressbalance
    64         spcvx::Vector{Float64}
    65         spcvy::Vector{Float64}
    66         restol::Float64
    67         reltol::Float64
    68         abstol::Float64
    69         maxiter::Int64
    70 end
    71 function Stressbalance() #{{{
    72         return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0., 0., 0., 0)
    73 end# }}}
    74 mutable struct Constants
    75         g::Float64
    76         yts::Float64
    77 end
    78 function Constants() #{{{
    79         return Constants( 9.81,  365*24*3600.)
    80 end# }}}
    81 mutable struct Materials
    82         rho_ice::Float64
    83         rho_water::Float64
    84         rho_freshwater::Float64
    85         mu_water::Float64
    86         heatcapacity::Float64
    87         latentheat::Float64
    88         thermalconductivity::Float64
    89         temperateiceconductivity::Float64
    90         effectiveconductivity_averaging::Int64
    91         meltingpoint::Float64
    92         beta::Float64
    93         mixed_layer_capacity::Float64
    94         thermal_exchange_velocity::Float64
    95         rheology_B::Vector{Float64}
    96         rheology_n::Vector{Float64}
    97         rheology_law::String
    98 end
    99 function Materials() #{{{
    100         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")
    101 end# }}}
    102 
    103 mutable struct model
    104         mesh::AbstractMesh
    105         geometry::Geometry
    106         mask::Mask
    107         materials::Materials
    108         initialization::Initialization
    109         stressbalance::Stressbalance
    110         constants::Constants
    111 end
    112 function model() #{{{
    113         return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
    114 end#}}}
    115 function Base.show(io::IO, md::model)# {{{
    116 
    117         compact = get(io, :compact, false)
    118 
    119         println(io,"Model:")
    120         @printf "%19s: %-22s -- %s\n" "mesh" typeof(md.mesh) "mesh properties"
    121         @printf "%19s: %-22s -- %s\n" "geometry" typeof(md.geometry) "surface elevation, bedrock topography, ice thickness,..."
    122         @printf "%19s: %-22s -- %s\n" "mask" typeof(md.mask) "defines grounded and floating regions"
    123         @printf "%19s: %-22s -- %s\n" "materials" typeof(md.materials) "material properties"
    124         @printf "%19s: %-22s -- %s\n" "initialization" typeof(md.initialization) "initial state"
    125         @printf "%19s: %-22s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
    126         @printf "%19s: %-22s -- %s\n" "constants" typeof(md.constants) "physical constants"
    127 
    128 end# }}}
    129 
    130 #utils
    131 function issmdir() #{{{
    132         issmdir = ENV["ISSM_DIR"]
    133 
    134         if isempty(issmdir)
    135                 error("Could not determine the location of ISSM")
    136         else
    137                 return issmdir
    138         end
    139 end#}}}
    140 function archread(filename::String,variablename::String) #{{{
    141 
    142         #initialize variables
    143         found = false
    144 
    145         #open file
    146         output = open(filename, "r") do f
    147 
    148                 while !eof(f)
    149                         reclen  = bswap(read(f, Int32))
    150                         rectype = bswap(read(f, Int32))
    151                         if rectype!=1
    152                                 error("Expected variable of type string")
    153                         else
    154                                 fieldname_length = bswap(read(f, Int32))
    155                                 field_name = String(read(f, fieldname_length))
    156                         end
    157                         rec_length = bswap(read(f, Int32))
    158                         field_type = bswap(read(f, Int32))
    159                         if field_type==2
    160                                 data = bswap(read(f, Float64))
    161                         elseif field_type==3
    162                                 rows = bswap(read(f, Int32))
    163                                 cols = bswap(read(f, Int32))
    164                                 data = reinterpret(Float64, read(f, sizeof(Float64)*rows*cols))
    165                                 data .= ntoh.(data)
    166                                 data = reshape(data, (rows,cols))
    167                                 data = collect(data)
    168                                 if cols == 1
    169                                         data = vec(data)
    170                                 end
    171                         else
    172                                 error("Error: Encountered invalid field type when reading data.")
    173                         end
    174 
    175                         if field_name == variablename
    176                                 found = true
    177                                 return data
    178                         end
    179                 end
    180         end
    181 
    182         return output
    183 end# }}}
    184 function InterpFromMeshToMesh2d(index::Array,x::Vector,y::Vector,data::Vector,xout::Vector,yout::Vector) #{{{
    185         #prepare input arrays
    186         nods = Cint(length(x))
    187         nels = Cint(size(index,1))
    188         nods_interp = Cint(length(xout))
    189         Cindex=Array{Cint,1}(undef,length(index))
    190         for i in 1:size(index,1)
    191                 for j in 1:3
    192                         Cindex[(i-1)*3+j] = Int32(index[i,j])
    193                 end
    194         end
    195         Cx    = Array{Cdouble,1}(undef,nods)
    196         Cy    = Array{Cdouble,1}(undef,nods)
    197         Cdata = Array{Cdouble,1}(undef,nods)
    198         for i in 1:nods
    199                 Cx[i]    = x[i]
    200                 Cy[i]    = y[i]
    201                 Cdata[i] = data[i]
    202         end
    203         Cxout = Array{Cdouble,1}(undef,nods_interp)
    204         Cyout = Array{Cdouble,1}(undef,nods_interp)
    205         for i in 1:nods_interp
    206                 Cxout[i] = xout[i]
    207                 Cyout[i] = yout[i]
    208         end
    209 
    210         Cdataout = Vector{Float64}(undef,nods_interp)
    211 
    212         #This is not working....
    213         rc=ccall( (:InterpFromMeshToMesh2dx,"libISSMCore"),
    214                                 Cint, (Ptr{Ptr{Cdouble}},Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint),
    215                                 Ref(Ref(Cdataout)), Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels,
    216                                 Ref(Cdata), nods, 1, Ref(Cxout), Ref(Cyout), nods_interp)
    217 
    218         #Process output
    219         dataout = Vector{Float64}(undef,nods_interp)
    220         for i in 1:nods_interp
    221                 dataout[i] = Cdataout[i]
    222         end
    223 
    224         return dataout
    225 end #}}}
    226 function solve(md::model,solution::String) #{{{
    227 
    228         if solution=="sb" || solution=="Stressbalance"
    229                 solutionstring = "StressbalanceSolution"
    230         else
    231                 error("solutionstring "*solution*" not supported!");
    232         end
    233 
    234         IssmCore(md)
    235 
    236         return md
    237 end #}}}

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

    r26500 r26522  
    1 include("issm.jl")
     1include("classes.jl")
    22include("exp.jl")
    33

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

    r26492 r26522  
    22
    33using GLMakie
     4include("classes.jl")
     5using .ISSM
    46
    5 function plotmodel(md::model,data)
     7function plotmodel(md::ISSM.model,data::Vector)
    68
    79        vertexcolor  = :black

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

    r26518 r26522  
    1 include("issm.jl")
     1include("classes.jl")
    22include("issmenums.jl")
    33
    … …  
    209209
    210210        #Solve
    211         error("not there yet")
    212211        Stressbalance(femmodel)
    213212
    … …  
    317316end#}}}
    318317function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{
    319         error("stop")
    320 end#}}}
     318
     319        #load constraints from model
     320        spcvx = md.stressbalance.spcvx
     321        spcvy = md.stressbalance.spcvy
     322
     323        count = 1
     324        for i in 1:md.mesh.numberofvertices
     325                if ~isnan(spcvx[i])
     326                        push!(constraints,Constraint(count,i,1,spcvx[i]))
     327                        count+=1
     328                end
     329                if ~isnan(spcvy[i])
     330                        push!(constraints,Constraint(count,i,2,spcvy[i]))
     331                        count+=1
     332                end
     333        end
     334
     335end#}}}
     336function Stressbalance(femmodel::FemModel)
     337        error("STOP")
     338end
    321339
    322340#Element functions

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

    r26500 r26522  
    11#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia
    2 include("triangle.jl")
    3 include("parameterization.jl")
    4 include("solve.jl")
     2include("issm.jl")
     3using .ISSM
    54
    6 md = model()
    7 md = triangle(md,"../../test/Exp/Square.exp",180000.)
    8 md = setmask(md,"all","")
     5md = ISSM.model()
     6md = ISSM.triangle(md,"../../test/Exp/Square.exp",180000.)
     7md = ISSM.setmask(md,"all","")
    98
    109#Geometry
    … …  
    2120
    2221#Initial velocity
    23 x     = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x")
    24 y     = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
    25 vx    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
    26 vy    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
    27 index = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index")
     22x     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","x")
     23y     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
     24vx    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
     25vy    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
     26index = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","index")
    2827md.initialization.vx=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
    2928md.initialization.vy=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
    … …  
    3635md.stressbalance.abstol=NaN
    3736
    38 md=solve(md,"Stressbalance")
     37#Boundary conditions
     38md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices)
     39md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices)
     40pos = findall(md.mesh.vertexonboundary)
     41md.stressbalance.spcvx[pos] .= 0.0
     42md.stressbalance.spcvy[pos] .= 0.0
     43
     44md=ISSM.solve(md,"Stressbalance")

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

    r26500 r26522  
    5151
    5252using Printf #needed for sprintf
    53 include("issm.jl")
     53include("classes.jl")
    5454include("exp.jl")
    5555
    … …  
    162162        points    = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio_out.pointlist,    (2,Int(ctio_out.numberofpoints)), own=true))'
    163163        triangles = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.trianglelist, (3,Int(ctio_out.numberoftriangles)), own=true))' .+1
     164        segments  = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.segmentlist,  (2,Int(ctio_out.numberofsegments)), own=true))' .+1
     165       
     166        #assign output
    164167        md.mesh = Mesh2dTriangle()
    165168        md.mesh.numberofvertices = ctio_out.numberofpoints
    … …  
    168171        md.mesh.y                = points[:,2]
    169172        md.mesh.elements         = triangles
     173        md.mesh.segments         = segments
     174
     175        #post processing
     176        md.mesh.vertexonboundary = zeros(Bool,md.mesh.numberofvertices)
     177        md.mesh.vertexonboundary[md.mesh.segments] .= true
    170178
    171179   return md
Note: See TracChangeset for help on using the changeset viewer.