Add error check WIP

scripts_future_API/tm_stokes_wip.jl

@@ -47,7 +47,10 @@ function main()
                            y = (VBC(x_bc, y_bc, 0.0), VBC(x_bc, y_bc, 0.0)),
                            z = (free_slip, free_surface))
 
-    gravity = (x=0.0, y=0.0, z=0.0)
+    ρg(x, y, z) = 0.0
+    gravity = (x=FunctionField(ρg, grid, (Vertex(), Center(), Center())),
+               y=FunctionField(ρg, grid, (Center(), Vertex(), Center())),
+               z=FunctionField(ρg, grid, (Center(), Center(), Vertex())))
 
     # numerics
     niter   = 10maximum(size(grid))
@@ -111,8 +114,9 @@ function main()
     for iter in 1:niter
         advance_iteration!(model, 0.0, 1.0; async=false)
         if (iter % ncheck == 0)
+            evaluate_error!(model; async=false)
             println("iter/nx = $(iter/maximum(size(grid)))")
-            evaluate_error(model; async=false)
+            println("  err = [Pr $(maximum(abs.(model.fields.r_Pr))), V.x $(maximum(abs.(model.fields.r_V.x))), V.y $(maximum(abs.(model.fields.r_V.y))), V.z $(maximum(abs.(model.fields.r_V.z)))]")
         end
         # if iter % ncheck == 0
         #     plt.Pr[3][] = interior(model.fields.Pr)[:, size(grid, 2)÷2, :]

src/Models/full_stokes/isothermal/boundary_conditions.jl

@@ -44,13 +44,13 @@ for (dim, val) in _COORDINATES
 
     ex_res_vel = ntuple(length(_COORDINATES)) do I
         kind = I == val ? :(FullCell) : :(HalfCell)
-        Expr(:(=), Symbol(:V, _COORDINATES[I][1]), :(DirichletBC{$kind}(0.0)))
+        Expr(:(=), Symbol(:r_V, _COORDINATES[I][1]), :(DirichletBC{$kind}(0.0)))
     end
 
     ex_res_vel = Expr(:tuple, ex_res_vel...)
 
-    ex_res_slip_vel = Expr(:(=), Symbol(:V, dim), :(DirichletBC{FullCell}(0.0)))
-    ex_res_slip_vel = Expr(:tuple, ex_slip_vel)
+    ex_res_slip_vel = Expr(:(=), Symbol(:r_V, dim), :(DirichletBC{FullCell}(0.0)))
+    ex_res_slip_vel = Expr(:tuple, ex_res_slip_vel)
 
     @eval begin
         extract_stress_bc(::Val{$val}, bc::BoundaryCondition{Traction}) = $ex_tr

src/Models/full_stokes/isothermal/isothermal.jl

@@ -1,7 +1,7 @@
 module Isothermal
 
 export BoundaryCondition, Traction, Velocity, Slip
-export IsothermalFullStokesModel, advance_iteration!, advance_timestep!
+export IsothermalFullStokesModel, advance_iteration!, advance_timestep!, evaluate_error!
 
 using FastIce.Architectures
 using FastIce.Physics
@@ -64,9 +64,9 @@ function make_fields_mechanics(backend, grid::CartesianGrid{3})
                z=Field(backend, grid, (Center(), Center(), Vertex()); halo=1)),
             # residual
             r_Pr=Field(backend, grid, Center(); halo=0),
-            r_V=(x=Field(backend, grid, (Vertex(), Center(), Center()); halo=0),
-                 y=Field(backend, grid, (Center(), Vertex(), Center()); halo=0),
-                 z=Field(backend, grid, (Center(), Center(), Vertex()); halo=0)))
+            r_V=(x=Field(backend, grid, (Vertex(), Center(), Center()); halo=1),
+                 y=Field(backend, grid, (Center(), Vertex(), Center()); halo=1),
+                 z=Field(backend, grid, (Center(), Center(), Vertex()); halo=1)))
 end
 
 function IsothermalFullStokesModel(;
@@ -126,14 +126,16 @@ function advance_iteration!(model::IsothermalFullStokesModel, t, Δt; async=true
     return
 end
 
-function evaluate_error(model::IsothermalFullStokesModel; async=true)
+function evaluate_error!(model::IsothermalFullStokesModel; async=true)
     (; Pr, τ, V, r_Pr, r_V) = model.fields
     ρg              = model.gravity
     hide_boundaries = model.hide_boundaries
     outer_width     = model.outer_width
 
+    Δ = NamedTuple{(:x, :y, :z)}(spacing(model.grid))
+
     launch!(model.arch, model.grid, compute_residuals! => (r_V, r_Pr, Pr, τ, V, ρg, model.grid, Δ);
-            location=Vertex(), expand=1, boundary_conditions=model.boundary_conditions.residuals_vel, hide_boundaries, outer_width)
+            location=Vertex(), boundary_conditions=model.boundary_conditions.residuals_vel, hide_boundaries, outer_width)
 
     async || synchronize(backend(model.arch))
     return