code runs, makes NaN

CliMA · Apr 4, 2024 · 8536daf · 8536daf
1 parent ad4f2a0
commit 8536daf
Show file tree

Hide file tree

Showing 7 changed files with 61 additions and 143 deletions.
diff --git a/docs/src/APIs/shared_utilities.md b/docs/src/APIs/shared_utilities.md
@@ -60,7 +60,6 @@ ClimaLand.boundary_var_types
 ClimaLand.make_tendency_jacobian
 ClimaLand.make_update_jacobian
 ClimaLand.set_dfluxBCdY!
-ClimaLand.AbstractTridiagonalW
 ClimaLand.get_drivers
 ```
 

diff --git a/experiments/Manifest.toml b/experiments/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "4616031e611fded1c2ac1e5271a5018186edc260"
+project_hash = "0bed51ae7dc5ef45124afa7b76d700cdefd4f573"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "016833eb52ba2d6bea9fcb50ca295980e728ee24"
@@ -358,6 +358,12 @@ git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1"
 uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
 version = "0.1.12"
 
+[[deps.CodeTracking]]
+deps = ["InteractiveUtils", "UUIDs"]
+git-tree-sha1 = "c0216e792f518b39b22212127d4a84dc31e4e386"
+uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
+version = "1.3.5"
+
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
 git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73"
@@ -514,6 +520,12 @@ version = "1.0.0"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.Debugger]]
+deps = ["CodeTracking", "Crayons", "Highlights", "InteractiveUtils", "JuliaInterpreter", "Markdown", "REPL"]
+git-tree-sha1 = "6798d210d8277b96cee6a6af387ed249b6e744ab"
+uuid = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
+version = "0.7.8"
+
 [[deps.DelaunayTriangulation]]
 deps = ["DataStructures", "EnumX", "ExactPredicates", "Random", "SimpleGraphs"]
 git-tree-sha1 = "d4e9dc4c6106b8d44e40cd4faf8261a678552c7c"
@@ -986,6 +998,12 @@ git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
 uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
 version = "2.8.1+1"
 
+[[deps.Highlights]]
+deps = ["DocStringExtensions", "InteractiveUtils", "REPL"]
+git-tree-sha1 = "0341077e8a6b9fc1c2ea5edc1e93a956d2aec0c7"
+uuid = "eafb193a-b7ab-5a9e-9068-77385905fa72"
+version = "0.5.2"
+
 [[deps.Hwloc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "ca0f6bf568b4bfc807e7537f081c81e35ceca114"
@@ -1173,6 +1191,12 @@ git-tree-sha1 = "3336abae9a713d2210bb57ab484b1e065edd7d23"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
 version = "3.0.2+0"
 
+[[deps.JuliaInterpreter]]
+deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
+git-tree-sha1 = "e9648d90370e2d0317f9518c9c6e0841db54a90b"
+uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
+version = "0.9.31"
+
 [[deps.JuliaNVTXCallbacks_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "af433a10f3942e882d3c671aacb203e006a5808f"

diff --git a/experiments/Project.toml b/experiments/Project.toml
@@ -8,6 +8,7 @@ ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"

diff --git a/src/shared_utilities/implicit_functions.jl b/src/shared_utilities/implicit_functions.jl
diff --git a/src/shared_utilities/implicit_tendencies.jl b/src/shared_utilities/implicit_tendencies.jl
@@ -1,5 +1,5 @@
 export make_tendency_jacobian,
-    make_update_jacobian, AbstractTridiagonalW, set_dfluxBCdY!
+    make_update_jacobian, set_dfluxBCdY!
 
 
 """
@@ -63,12 +63,3 @@ function set_dfluxBCdY!(
     ::AbstractBoundary,
     _...,
 )::Union{ClimaCore.Fields.FieldVector, Nothing} end
-
-"""
-    AbstractTridiagonalW
-
-An abstract type for tridiagonal Jacobian matrices.
-"""
-abstract type AbstractTridiagonalW end
-
-Base.similar(w::AbstractTridiagonalW) = w
diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
@@ -489,12 +489,30 @@ function ClimaLand.set_dfluxBCdY!(
     θ_rN = Domains.top_center_to_surface(θ_r)
     S_sN = Domains.top_center_to_surface(S_s)
 
+    # Extract the local geometry of the top face space
+    local_geometryN = Fields.local_geometry_field(KN)
+
     # Update dfluxBCdY at the top boundary in place
+    # Calculate the value and convert it to a Covariant3Vector
     @. p.soil.dfluxBCdY =
-        Geometry.WVector(KN * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
+        covariant3_unit_vector(local_geometryN) * (KN * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
     return nothing
 end
 
+"""
+    covariant3_unit_vector(local_geometry)
+
+A function to compute the unit vector in the direction of the normal
+to the surface.
+
+Adapted from ClimaAtmos.jl's unit_basis_vector_data function.
+"""
+function covariant3_unit_vector(local_geometry)
+    FT = Geometry.undertype(typeof(local_geometry))
+    data = FT(1) / Geometry._norm(Geometry.Covariant3Vector(FT(1)), local_geometry)
+    return Geometry.Covariant3Vector(data)
+end
+
 """
     ClimaLand.set_dfluxBCdY!(
         ::RichardsModel,
@@ -991,4 +1009,4 @@ boundary_var_types(
     model::RichardsModel{FT},
     bc::MoistureStateBC,
     ::ClimaLand.TopBoundary,
-) where {FT} = (FT, ClimaCore.Geometry.WVector{FT})
+) where {FT} = (FT, ClimaCore.Geometry.Covariant3Vector{FT})
diff --git a/src/standalone/Soil/rre.jl b/src/standalone/Soil/rre.jl
@@ -347,7 +347,7 @@ end
 
 
 """
-    ImplicitEquationJacobian{M, S} <: ClimaLand.AbstractTridiagonalW
+    ImplicitEquationJacobian{M, S}
 
 A struct containing the necessary information for constructing a tridiagonal
 Jacobian matrix (`W`) for solving Richards equation, treating only the vertical
@@ -361,7 +361,7 @@ block diagonal.
 Note that the diagonal, upper diagonal, and lower diagonal entry values
 are stored in this struct and updated in place.
 """
-struct ImplicitEquationJacobian{M, S} <: ClimaLand.AbstractTridiagonalW
+struct ImplicitEquationJacobian{M, S}
     "Jacobian matrix stored as a MatrixFields.FieldMatrix"
     matrix::M
     "Solver to use for solving the tridiagonal system"
@@ -409,6 +409,15 @@ function ImplicitEquationJacobian(Y::ClimaCore.Fields.FieldVector)
     return ImplicitEquationJacobian(matrix, solver)
 end
 
+Base.similar(w::ImplicitEquationJacobian) = w
+
+function LinearAlgebra.ldiv!(
+    x::Fields.FieldVector,
+    A::ImplicitEquationJacobian,
+    b::Fields.FieldVector,
+)
+    MatrixFields.field_matrix_solve!(A.solver, x, A.matrix, b)
+end
 
 """
     ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
@@ -449,7 +458,7 @@ function ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
             dfluxBCdY = p.soil.dfluxBCdY
             topBC_op = Operators.SetBoundaryOperator(
                 top = Operators.SetValue(dfluxBCdY),
-                bottom = Operators.SetValue(Geometry.WVector(zero(FT))),
+                bottom = Operators.SetValue(Geometry.Covariant3Vector(zero(FT))),
             )
             # Add term from top boundary condition before applying divergence
             # Note: need to pass 3D field on faces to `topBC_op`. Interpolating `K` to faces
@@ -465,9 +474,9 @@ function ClimaLand.make_update_jacobian(model::RichardsModel{FT}) where {FT}
                             θ_r,
                             S_s,
                         ),
-                    ) + MatrixFields.UpperBidiagonalMatrixRow(
+                    ) + MatrixFields.UpperDiagonalMatrixRow(
                         topBC_op(
-                            Geometry.WVector(zero(interpc2f_op(p.soil.K))),
+                            Geometry.Covariant3Vector(zero(interpc2f_op(p.soil.K))),
                         ),
                     )
                 )