Update volfrac script to mask API

PTsolvers · Apr 25, 2024 · c2228e3 · c2228e3
1 parent 4bf65aa
commit c2228e3
Showing 1 changed file with 15 additions and 14 deletions.
diff --git a/...API/test_levelsets_volumefractions_mpi.jl → scripts_future_API/test_levelsets_volfrac.jl b/...API/test_levelsets_volumefractions_mpi.jl → scripts_future_API/test_levelsets_volfrac.jl
@@ -1,6 +1,7 @@
 using Chmy.Architectures
 using Chmy.Grids
 using Chmy.Fields
+using Chmy.GridOperators
 using Chmy.KernelLaunch
 
 using FastIce.Writers
@@ -11,13 +12,13 @@ using JLD2
 using KernelAbstractions
 using CairoMakie
 # using CUDA
-using AMDGPU
-AMDGPU.allowscalar(false)
+# using AMDGPU
+# AMDGPU.allowscalar(false)
 
 # Select backend
-# backend = CPU()
+backend = CPU()
 # backend = CUDABackend()
-backend = ROCBackend()
+# backend = ROCBackend()
 
 do_h5_save = true
 
@@ -97,8 +98,8 @@ function main_vavilov(backend=CPU())
 
     data_elevation = extract_dem(arch, data_path)
 
-    nx, ny = 254, 254
-    nz     = max(conv(ceil(Int, data_elevation.lz / data_elevation.lx * nx), 30), 62)
+    nx, ny = 126, 126
+    nz     = max(conv(ceil(Int, data_elevation.lz / data_elevation.lx * nx), 30), 30)
     resol  = (nx, ny, nz)
 
     run_simulation(data_elevation..., resol...)
@@ -118,24 +119,24 @@ end
     # init fields
     Ψ = (na=Field(arch, grid, Vertex()),
          ns=Field(arch, grid, Vertex()))
-    wt = (na=volfrac_field(arch, grid),
-          ns=volfrac_field(arch, grid))
+    wt = (na=FieldMask(arch, grid),
+          ns=FieldMask(arch, grid))
     # comput level set
     compute_levelset_from_dem!(arch, launch, Ψ.na, surf, grid_2d, grid)
     compute_levelset_from_dem!(arch, launch, Ψ.ns, bed, grid_2d, grid)
     invert_levelset!(arch, launch, Ψ.ns, grid) # invert level set to set what's below the DEM surface as inside
     # volume fractions
     for phase in eachindex(Ψ)
-        compute_volfrac_from_levelset!(arch, launch, wt[phase], Ψ[phase], grid)
+        ω_from_ψ!(arch, launch, wt[phase], Ψ[phase], grid)
     end
 
     # compute physics or else
 
     # postprocessing
-    wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing
-    wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing
-    gather!(arch, wt_na_c, wt.na.c)
-    gather!(arch, wt_ns_c, wt.ns.c)
+    wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.ccc)) .* dims(topo)) : nothing
+    wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.ccc)) .* dims(topo)) : nothing
+    gather!(arch, wt_na_c, wt.na.ccc)
+    gather!(arch, wt_ns_c, wt.ns.ccc)
     # visualise
     if me == 0
         dem_bed  = Array(interior(bed))
@@ -169,7 +170,7 @@ end
 
     if do_h5_save
         h5names = String[]
-        fields = Dict("wt_na" => wt.na.c, "wt_ns" => wt.ns.c)
+        fields = Dict("wt_na" => wt.na.ccc, "wt_ns" => wt.ns.ccc)
         outdir = "out_visu_mpi"
         (me == 0) && mkpath(outdir)