Update model

scripts_future_API/tm_stokes_wip.jl

@@ -4,6 +4,7 @@ using FastIce.Fields
 using FastIce.Utils
 using FastIce.BoundaryConditions
 using FastIce.Models.FullStokes.Isothermal
+using FastIce.Physics
 
 using KernelAbstractions
 
@@ -40,32 +41,41 @@ iter_params = (
     Δτ = ( Pr = r / θ_dτ, τ = (xx = dτ_r, yy = dτ_r, zz = dτ_r, xy = dτ_r, xz = dτ_r, yz = dτ_r), V = (x = nudτ, y = nudτ, z = nudτ)),
 )
 
-physics = (rheology = GlensLawRheology(1.0, 1))
+backend = CPU()
+
+physics = (rheology = GlensLawRheology(1), )
+other_fields = (
+    A = Field(backend, grid, Center()),
+)
+
+init_incl(x, y, z, x0, y0, z0, r, ηi, ηm) = ifelse((x-x0)^2 + (y-y0)^2 + (z-z0)^2 < r^2, ηi, ηm)
+set!(other_fields.A, grid, init_incl; continuous = true, parameters = (x0 = 0.0, y0 = 0.0, z0 = 0.5, r = 0.2, ηi = 1e-1, ηm = 1.0))
 
 model = IsothermalFullStokesModel(;
-    backend = CPU(),
+    backend,
     grid,
+    physics,
     boundary_conditions,
     iter_params,
+    other_fields
 )
 
 set!(model.fields.Pr, 0.0)
 foreach(x -> set!(x, 0.0), model.fields.τ)
 Isothermal.apply_bcs!(model.backend, model.grid, model.fields, model.boundary_conditions.stress)
 
-foreach(x -> set!(x, 0.0), model.fields.V)
+# foreach(x -> set!(x, 0.0), model.fields.V)
 
-# set!(model.fields.V.x, grid, (x, y, z, ebg) -> -x*ebg; continuous=true, parameters = (ebg, ))
-# set!(model.fields.V.y, grid, (x, y, z, ebg) ->  y*ebg; continuous=true, parameters = (ebg, ))
-# set!(model.fields.V.z, 0.0)
+set!(model.fields.V.x, grid, (x, y, z, ebg) -> -x*ebg; continuous=true, parameters = (ebg, ))
+set!(model.fields.V.y, grid, (x, y, z, ebg) ->  y*ebg; continuous=true, parameters = (ebg, ))
+set!(model.fields.V.z, 0.0)
 Isothermal.apply_bcs!(model.backend, model.grid, model.fields, model.boundary_conditions.velocity)
 
-init_incl(x, y, z, x0, y0, z0, r, ηi, ηm) = ifelse((x-x0)^2 + (y-y0)^2 + (z-z0)^2 < r^2, ηi, ηm)
-set!(model.fields.η, grid, init_incl; continuous = true, parameters = (x0 = 0.0, y0 = 0.0, z0 = 0.5, r = 0.2, ηi = 1e-1, ηm = 1.0))
+set!(model.fields.η, other_fields.A)
 extrapolate!(data(model.fields.η))
 
 
-for it in 1:10
+for it in 1:100
     advance_iteration!(model, 0.0, 1.0; async = false)
 end
 

src/FastIce.jl

@@ -3,6 +3,7 @@ module FastIce
 using KernelAbstractions
 
 include("utils.jl")
+include("macros.jl")
 include("logging.jl")
 include("grids.jl")
 include("fields.jl")

src/models/full_stokes/common/macros.jl → src/macros.jl

@@ -1,3 +1,9 @@
+module Macros
+
+export @all, @inn, @inn_x, @inn_y, @inn_z, @inn_xy, @inn_xz, @inn_yz
+export @∂_x, @∂_y, @∂_z, @∂_xi, @∂_yi, @∂_zi
+export @av_xy, @av_xz, @av_yz, @av_xyi, @av_xzi, @av_yzi
+
 macro all(A) esc(:($A[ix, iy, iz])) end
 
 macro inn(A)   esc(:($A[ix+1, iy+1, iz+1])) end
@@ -24,4 +30,6 @@ macro av_yzi(A) esc(:(0.25 * ($A[ix+1, iy, iz] + $A[ix+1, iy+1, iz] + $A[ix+1, i
 
 macro av_xy(A) esc(:(0.25 * ($A[ix, iy, iz] + $A[ix+1, iy, iz] + $A[ix+1, iy+1, iz] + $A[ix, iy+1, iz]))) end
 macro av_xz(A) esc(:(0.25 * ($A[ix, iy, iz] + $A[ix+1, iy, iz] + $A[ix+1, iy, iz+1] + $A[ix, iy, iz+1]))) end
-macro av_yz(A) esc(:(0.25 * ($A[ix, iy, iz] + $A[ix, iy+1, iz] + $A[ix, iy+1, iz+1] + $A[ix, iy, iz+1]))) end
+macro av_yz(A) esc(:(0.25 * ($A[ix, iy, iz] + $A[ix, iy+1, iz] + $A[ix, iy+1, iz+1] + $A[ix, iy, iz+1]))) end
+
+end

src/models/full_stokes/isothermal/isothermal.jl

@@ -15,7 +15,7 @@ function default_physics(::Type{T}) where T
     return (
         equation_of_state=default(IncompressibleIceEOS{T}),
         thermal_properties=default(IceThermalProperties{T}),
-        rheology=default(GlensLawRheology{T, Int64})
+        rheology=default(GlensLawRheology{Int64})
     )
 end
 
@@ -47,13 +47,17 @@ function make_fields_mechanics(backend, grid)
     )
 end
 
-function IsothermalFullStokesModel(; backend, grid, boundary_conditions, physics=nothing, iter_params, fields=nothing)
+function IsothermalFullStokesModel(; backend, grid, boundary_conditions, physics=nothing, iter_params, fields=nothing, other_fields=nothing)
     if isnothing(fields)
         mechanic_fields = make_fields_mechanics(backend, grid)
         rheology_fields = (η = Field(backend, grid, Center(), (1, 1, 1)),)
         fields = merge(mechanic_fields, rheology_fields)
     end
 
+    if !isnothing(other_fields)
+        fields = merge(fields, other_fields)
+    end
+
     if isnothing(physics)
         physics = default_physics(eltype(grid))
     end
@@ -237,13 +241,13 @@ function advance_iteration!(model::IsothermalFullStokesModel, t, Δt; async = tr
     set_bcs!(bcs) = apply_bcs!(model.backend, model.grid, model.fields, bcs)
 
     # stress
-    update_σ!(backend, 256, (nx + 2, ny + 2, nz + 2))(interior(Pr), interior(τ), V, η, Δτ, Δ)
+    update_σ!(backend, 256, (nx, ny, nz))(interior(Pr), interior(τ), V, η, Δτ, Δ)
     set_bcs!(model.boundary_conditions.stress)
     # velocity
     update_V!(backend, 256, (nx + 1, ny + 1, nz + 1))(interior(V), Pr, τ, η, Δτ, Δ)
     set_bcs!(model.boundary_conditions.velocity)
     # rheology
-    update_η!(backend, 256, (nx, ny, nz))(interior(η), τ, η_rh, η_rel)
+    update_η!(backend, 256, (nx, ny, nz))(interior(η), η_rh, η_rel, model.grid, model.fields)
     extrapolate!(data(η))
 
     async || synchronize(backend)

src/models/full_stokes/isothermal/kernels.jl

@@ -1,16 +1,16 @@
 using KernelAbstractions
 
-include("../common/macros.jl")
+using FastIce.Macros
 
-@kernel function update_η!(η, τ, η_rh, χ)
+@kernel function update_η!(η, η_rh, χ, grid, fields, args...)
     ix, iy, iz = @index(Global, NTuple)
-    @inbounds τII = sqrt(0.5 * (@inn(τ.xx)^2 + @inn(τ.yy)^2 + @inn(τ.zz)^2) + @av_xy(τ.xy)^2 + @av_xz(τ.xz)^2 + @av_yz(τ.yz)^2)
-    @inbounds @inn(η) = exp(log(@inn(η)) * (1 - χ) + log(η_rh(τII)) * χ)
+    ηt = η_rh(grid, ix, iy, iz, fields, args...)
+    @inbounds @all(η) = exp(log(@all(η)) * (1 - χ) + log(ηt) * χ)
 end
 
 @kernel function update_σ!(Pr, τ, V, η, Δτ, Δ)
     ix, iy, iz = @index(Global, NTuple)
-    @inbounds if ix <= size(Pr, 1) - 2 && iy <= size(Pr, 2) - 2 && iz <= size(Pr, 3) - 2
+    @inbounds if checkbounds(Bool, Pr, ix, iy, iz)
         ε̇xx = @∂_x(V.x) / Δ.x
         ε̇yy = @∂_y(V.y) / Δ.y
         ε̇zz = @∂_z(V.z) / Δ.z
@@ -35,15 +35,15 @@ end
 
 @kernel function update_V!(V, Pr, τ, η, Δτ, Δ)
     ix, iy, iz = @index(Global, NTuple)
-    @inbounds if ix <= size(V.x, 1) && iy <= size(V.x, 2) - 2 && iz <= size(V.x, 3) - 2
+    @inbounds if checkbounds(Bool, V.x, ix, iy, iz)
         ηx = max(η[ix, iy+1, iz+1], η[ix+1, iy+1, iz+1])
         @all(V.x) += ((-@∂_xi(Pr) + @∂_xi(τ.xx)) / Δ.x + @∂_y(τ.xy) / Δ.y + @∂_z(τ.xz) / Δ.z) / ηx * Δτ.V.x
     end
-    @inbounds if ix <= size(V.y, 1) - 2 && iy <= size(V.y, 2) && iz <= size(V.y, 3) - 2
+    @inbounds if checkbounds(Bool, V.y, ix, iy, iz)
         ηy = max(η[ix+1, iy, iz+1], η[ix+1, iy+1, iz+1])
         @all(V.y) += ((-@∂_yi(Pr) + @∂_yi(τ.yy)) / Δ.y + @∂_x(τ.xy) / Δ.x + @∂_z(τ.yz) / Δ.z) / ηy * Δτ.V.y
     end
-    @inbounds if ix <= size(V.z, 1) - 2 && iy <= size(V.z, 2) - 2 && iz <= size(V.z, 3)
+    @inbounds if checkbounds(Bool, V.z, ix, iy, iz)
         ηz = max(η[ix+1, iy+1, iz], η[ix+1, iy+1, iz+1])
         @all(V.z) += ((-@∂_zi(Pr) + @∂_zi(τ.zz)) / Δ.z + @∂_x(τ.xz) / Δ.x + @∂_y(τ.yz) / Δ.y) / ηz * Δτ.V.z
     end

src/physics.jl

@@ -5,6 +5,8 @@ export IncompressibleIceEOS, IceThermalProperties
 export IceRheology, GlensLawRheology, τ
 export default
 
+using FastIce.Macros
+
 abstract type AbstractPhysics end
 
 struct IncompressibleIceEOS{T}
@@ -23,15 +25,16 @@ default(::Type{IceThermalProperties{T}}) where {T} = IceThermalProperties(conver
 
 abstract type IceRheology end
 
-struct GlensLawRheology{T,I}
-    consistency::T
+struct GlensLawRheology{I}
     exponent::I
 end
 
-default(::Type{GlensLawRheology{T,I}}) where {T,I} = GlensLawRheology(convert(T, 2.4e-24), convert(I, 3))
+default(::Type{GlensLawRheology{I}}) where {I} = GlensLawRheology(convert(I, 3))
 
-@inline function (rh::GlensLawRheology{T})(τII::T) where {T}
-    return 0.5 / (rh.consistency * τII^(rh.exponent - 1))
+@inline function (rh::GlensLawRheology{T})(grid, ix, iy, iz, fields) where {T}
+    (; τ, A) = fields
+    @inbounds τII = sqrt(0.5 * (@inn(τ.xx)^2 + @inn(τ.yy)^2 + @inn(τ.zz)^2) + @av_xy(τ.xy)^2 + @av_xz(τ.xz)^2 + @av_yz(τ.yz)^2)
+    return 0.5 / (A[ix, iy, iz] * τII^(rh.exponent - 1))
 end
 
 end