Merge pull request #27 from PALEOtoolkit/romglb_test_update

add some 'romglb' tests that got missed

Project.toml

@@ -37,6 +37,7 @@ Plots = "1.38"
 Preferences = "1.3"
 SIMD = "3.4"
 SnoopPrecompile = "1.0"
+SparseDiffTools = "2.0"
 SpecialFunctions = "1.0, 2.0"
 Sundials = "4.0"
 TestEnv = "1.0"
@@ -52,9 +53,10 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PALEOcopse = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
 PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [targets]
-test = ["DataFrames", "DiffEqBase", "Documenter", "Downloads", "Logging", "PALEOcopse", "PALEOmodel", "Plots", "Sundials", "Test", "ZipFile"]
+test = ["DataFrames", "DiffEqBase", "Documenter", "Downloads", "Logging", "PALEOcopse", "PALEOmodel", "Plots", "Sundials", "SparseDiffTools", "Test", "ZipFile"]

examples/Project.toml

@@ -17,6 +17,7 @@ PALEOsediment = "e0a37952-6f01-4236-91ff-62fdc855f67b"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

examples/romglb/runtests.jl

@@ -3,6 +3,8 @@ using Logging
 using DiffEqBase
 using Sundials
 import DataFrames
+import SparseArrays
+import SparseDiffTools
 
 import PALEOboxes as PB
 
@@ -14,6 +16,7 @@ import PALEOcopse
 @testset "romglb examples" begin
 
 skipped_testsets = [
+    # "airsea_O2",
     # "O2_only", 
     # "P_O2",
     # "P_O2_S_Carb_open",
@@ -32,6 +35,81 @@ include("SedimentationRate_dev.jl")
 
 include("config_ocean_romglb_expts.jl")
 
+!("airsea_O2" in skipped_testsets) && @testset "airsea_O2" begin
+
+    model = PB.create_model_from_config(
+        joinpath(@__DIR__, "PALEO_examples_romglb_cfg.yaml"), "romglb_abiotic_O2";
+        modelpars=Dict(
+            "matdir"=>matdir,
+        )
+    )
+
+    # Test OceanBase domain configuration
+    @test PB.get_num_domains(model) == 6
+
+    global_domain = PB.get_domain(model, "global")
+    @test PB.get_length(global_domain) == 1
+
+    ocean_domain = PB.get_domain(model, "ocean")
+    @test PB.get_length(ocean_domain) == 79
+
+
+    # test OceanBase variables
+
+    initial_state, modeldata = PALEOmodel.initialize!(model)
+
+    ocean_modelcreated_vars_dict = Dict([(var.name, var) for var in PB.get_variables(ocean_domain, hostdep=false)])
+
+    println("ocean  model created variables after initialize!:")
+    for (name, var ) in ocean_modelcreated_vars_dict
+        println("\t", PB.fullname(var), " = ", PB.get_data(var, modeldata))
+    end
+
+    # bodge a test for ocean tracer with single non-zero cell
+    ocean_T_data = PB.get_data(PB.get_variable(ocean_domain,"T"), modeldata)
+    ocean_T_data .= 0.0
+    ocean_T_data[1] = 1.0
+    # update initial_state with our bodged values
+    initial_state = PALEOmodel.get_statevar(modeldata.solver_view_all)
+
+    # Check model derivative
+
+    PB.do_deriv(modeldata.dispatchlists_all)
+
+    println("state, sms variables after check model derivative:")
+    for (state_var, sms_var) in PB.IteratorUtils.zipstrict(PB.get_vars(modeldata.solver_view_all.stateexplicit), 
+                                    PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv))
+        println(PB.fullname(state_var), " ", PB.get_data(state_var, modeldata))
+        println(PB.fullname(sms_var), " ", PB.get_data(sms_var, modeldata))
+    end
+
+    # check conservation
+    ocean_T_sms_data = PB.get_data(PB.get_variable(ocean_domain,"T_sms"), modeldata)
+    sum_T_sms = sum(ocean_T_sms_data)
+    println("check conservation: sum(Tracer_sms)=",sum_T_sms)
+    @test abs(sum_T_sms) < 1e-15
+
+    # check jacobian sparsity calculation
+    jac_prototype = PALEOmodel.JacobianAD.calcJacobianSparsitySparsityTracing!(model, modeldata, initial_state, 0.0)
+    @test SparseArrays.nnz(jac_prototype) == 819
+    jac_proto_fill = PALEOmodel.SparseUtils.fill_sparse_jac(jac_prototype)
+    colors = SparseDiffTools.matrix_colors(copy(jac_proto_fill))
+    @test maximum(colors) == 23  
+
+    # integrate to approx steady state
+    paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
+    PALEOmodel.ODE.integrateForwardDiff(paleorun, initial_state, modeldata, (0, 1e5), solvekwargs=(reltol=1e-5,))  # first run includes JIT time
+
+    # check conservation
+    T_total = PB.get_data(paleorun.output, "ocean.T_total")
+    @test abs(T_total[1] - 1.0) < 1e-16
+    @test abs(T_total[end] - T_total[1]) < 1e-4
+
+    total_O2 = PB.get_data(paleorun.output, "global.total_O2")
+    @test abs(total_O2[end] - total_O2[1]) < 1e-7*total_O2[1]
+
+end
+
 !("O2_only" in skipped_testsets) && @testset "O2_only" begin
 
     model = PB.create_model_from_config(
@@ -216,3 +294,5 @@ end
 
 
 end
+
+nothing # so no output printed when run from REPL