Update MITgcm configs for PALEOboxes thread changes (#29)

* Update MITgcm configs for PALEOboxes thread changes Requires PALEOtoolkit/PALEOmodel.jl#105 PALEOtoolkit/PALEOboxes.jl#147 * Update Project.toml [compat] for PALEOmodel PALEOmodel v0.15.48 required for thread config changes
PALEOtoolkit · Dec 20, 2024 · 74fafaa · 74fafaa
1 parent f474b2d
commit 74fafaa
Show file tree

Hide file tree

Showing 12 changed files with 108 additions and 61 deletions.
diff --git a/examples/Project.toml b/examples/Project.toml
@@ -24,4 +24,4 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [compat]
-PALEOmodel = "0.15.8"
+PALEOmodel = "0.15.48"
diff --git a/examples/mitgcm/MITgcm_2deg8_COPDOM.yaml b/examples/mitgcm/MITgcm_2deg8_COPDOM.yaml
@@ -3,6 +3,7 @@ PO4MMbase:
     parameters:
         CIsotope: ScalarData 
         SIsotope: ScalarData  
+        threadsafe: false
     domains:
         global:
             # scalar domain
@@ -285,6 +286,7 @@ PO4MMcarbSCH4:
     parameters:
         CIsotope: IsotopeLinear # ScalarData
         SIsotope: IsotopeLinear
+        threadsafe: false
 
     domains:
         global:

diff --git a/examples/mitgcm/MITgcm_2deg8_PO4MMbase.jl b/examples/mitgcm/MITgcm_2deg8_PO4MMbase.jl
@@ -17,7 +17,8 @@ n_inner = 2
 # model = config_mitgcm_expts("PO4MMbase", ""); toutputs = [0.0, 1.0, 10.0, 100.0, 995,0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0] #, 1000.0, 1000.5]
 
 model = PB.create_model_from_config(
-    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase"
+    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase";
+    modelpars=Dict("threadsafe"=>use_threads),
 )
 
 toutputs = [0.0, 0.25, 0.5, 0.75, 1.0, 10.0] 
@@ -30,8 +31,17 @@ tstep_imp = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr
 output_filename = ""
 # output_filename = "MITgcm_PO4MMbase2deg8_3000yr_20210202"
 
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
+
 pickup_output = nothing
-initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output)  
+initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output)  
 pickup_output = nothing
 
 paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
@@ -74,7 +84,7 @@ else
 end
 
 
-isempty(output_filename) || PALEOmodel.OutputWriters.save_jld2(paleorun.output, output_filename)
+isempty(output_filename) || PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename)
 
 
 show(PB.show_variables(paleorun.model), allrows=true)

diff --git a/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl b/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl
@@ -13,7 +13,8 @@ include("plot_mitgcm.jl")
 use_threads = true
 
 model = PB.create_model_from_config(
-    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4"
+    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4";
+    modelpars=Dict("threadsafe"=>use_threads),
 )
 
 toutputs = [0, 1.0, 10.0] # , 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0]
@@ -31,8 +32,16 @@ tstep = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr
 # output_filename = "MITgcm_PO4MMcarbSCH42deg8FP64_2000yr_20230422"
 output_filename = ""
 
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
 pickup_output = nothing
-initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output)  
+initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output)  
 pickup_output = nothing
 
 paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
@@ -64,7 +73,7 @@ else
     @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep)
 end
 
-isempty(output_filename) || PALEOmodel.OutputWriters.save_jld2(paleorun.output, output_filename)
+isempty(output_filename) || PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename)
 
 show(PB.show_variables(paleorun.model), allrows=true)
 println()

diff --git a/examples/mitgcm/MITgcm_2deg8_abiotic.jl b/examples/mitgcm/MITgcm_2deg8_abiotic.jl
@@ -14,7 +14,8 @@ use_threads = false
 
 
 model = PB.create_model_from_config(
-    joinpath(@__DIR__, "MITgcm_2deg8_abiotic.yaml"), "abiotic_O2"
+    joinpath(@__DIR__, "MITgcm_2deg8_abiotic.yaml"), "abiotic_O2";
+    modelpars=Dict("threadsafe"=>use_threads),
 )
 
 toutputs_relative = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0]
@@ -31,11 +32,20 @@ outfile_root = ""
 
 toutputs = []
 
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
+
 for iseg in 1:num_segments
     if iseg > 1
         !isempty(outfile_root) || error("outfile_root is empty")
         pickup_filename = build_outfilename(outfile_root, iseg-1)
-        pickup_output = PALEOmodel.OutputWriters.load_jld2!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename)
+        pickup_output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename)
 
         tstart = PB.get_data(pickup_output, "ocean.tmodel")[end] 
 
@@ -48,7 +58,7 @@ for iseg in 1:num_segments
     global modeldata
     global toutputs
 
-    initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output) 
+    initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output) 
     pickup_output = nothing
 
     toutputs = toutputs_relative .+ tstart
@@ -72,7 +82,7 @@ for iseg in 1:num_segments
 
     if !isempty(outfile_root)
         output_filename = build_outfilename(outfile_root, iseg)
-        PALEOmodel.OutputWriters.save_jld2(paleorun.output, output_filename)
+        PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename)
     end
 end
 

diff --git a/examples/mitgcm/MITgcm_2deg8_abiotic.yaml b/examples/mitgcm/MITgcm_2deg8_abiotic.yaml
@@ -1,7 +1,8 @@
 # MITgcm 2.8deg O2-only atmosphere-ocean test configuration
 abiotic_O2:
     parameters:
-        CIsotope: ScalarData  
+        CIsotope: ScalarData
+        threadsafe: false 
     domains:
         global:
             # scalar domain

diff --git a/examples/mitgcm/MITgcm_2deg8_test.jl b/examples/mitgcm/MITgcm_2deg8_test.jl
@@ -11,11 +11,23 @@ include("config_mitgcm_expts.jl")
 
 use_threads = true
 
-# model = config_mitgcm_expts("abiotic_O2", ""); toutputs = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0]
+# model = PB.create_model_from_config(
+#     joinpath(@__DIR__, "MITgcm_2deg8_abiotic.yaml"), "abiotic_O2"
+# )
+# toutputs = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0]
 
-# model = config_mitgcm_expts("PO4MMbase", ""); toutputs = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0, 99.5, 100.0] #, 1000.0, 1000.5]
 
-model = config_mitgcm_expts("PO4MMcarbSCH4", ""); toutputs = [0, 1.0] # , 10.0, 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0]
+# model = PB.create_model_from_config(
+#     joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase"
+# )
+# toutputs = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0, 99.5, 100.0] #, 1000.0, 1000.5]
+
+
+model = PB.create_model_from_config(
+    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4";
+    modelpars=Dict("threadsafe"=>use_threads),
+)
+toutputs = [0, 1.0] # , 10.0, 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0]
 # start with low oxygen to test marine sulphur system
 PB.set_variable_attribute!(model, "atm", "O2", :initial_value, 0.1*3.71e19)
 PB.set_variable_attribute!(model, "ocean", "O2", :initial_value, 0.1*0.2054)
@@ -26,8 +38,16 @@ tstep = transportMITgcm.par_Aimp_deltat[]/PB.Constants.k_secpyr
 
 @info "using tstep=$tstep yr"
 
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
 
-initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads)  
+initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier)  
 
 run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
 

diff --git a/examples/mitgcm/MITgcm_ECCO_COPDOM.yaml b/examples/mitgcm/MITgcm_ECCO_COPDOM.yaml
@@ -4,7 +4,8 @@
 PO4MMbase:
     parameters:
         CIsotope: ScalarData 
-        SIsotope: ScalarData  
+        SIsotope: ScalarData
+        threadsafe: false 
     domains:
         global:
             # scalar domain

diff --git a/examples/mitgcm/MITgcm_ECCO_abiotic.jl b/examples/mitgcm/MITgcm_ECCO_abiotic.jl
@@ -13,7 +13,11 @@ include("plot_mitgcm.jl")
 use_threads = true
 do_benchmarks = false
 
-model = config_mitgcm_expts("PO4MMbaseECCO", ""); toutputs_relative = [0.0, 1.0, 10.0, 100.0]
+model = PB.create_model_from_config(
+    joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase";
+    modelpars=Dict("threadsafe"=>use_threads),
+)
+toutputs_relative = [0.0, 1.0, 10.0, 100.0]
 
 transportMITgcm = PB.get_reaction(model, "ocean", "transportMITgcm")
 tstep = transportMITgcm.par_Aimp_deltat[]/PB.Constants.k_secpyr
@@ -25,17 +29,27 @@ num_segments = 20
 outfile_root = "MITgcm_PO4MMbaseECCO_100yr_20210130"
 
 toutputs=[]
+
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
+
 for iseg in 1:num_segments
 # for iseg in 13:num_segments
     if iseg > 1
         pickup_filename = build_outfilename(outfile_root, iseg-1)
-        pickup_output = PALEOmodel.OutputWriters.load_jld2!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename)
+        pickup_output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename)
         tstart = PB.get_data(pickup_output, "ocean.tmodel")[end]
     else
         pickup_output = nothing
         tstart = 0.0
     end
-    initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output)  
+    initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output)  
     pickup_output = nothing
 
     toutputs = toutputs_relative .+ tstart
@@ -59,7 +73,7 @@ for iseg in 1:num_segments
     end
 
     output_filename = build_outfilename(outfile_root, iseg)
-    PALEOmodel.OutputWriters.save_jld2(paleorun.output, output_filename)
+    PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename)
 end
 
 

diff --git a/examples/mitgcm/MITgcm_ECCO_abiotic_test.jl b/examples/mitgcm/MITgcm_ECCO_abiotic_test.jl
@@ -14,7 +14,8 @@ use_threads = true
 do_benchmarks = false
 
 model = PB.create_model_from_config(
-    joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase"
+    joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase";
+    modelpars=Dict("threadsafe"=>use_threads),
 )
 
 toutputs = [0.0, 1.0] # , 10.0, 100.0]
@@ -27,8 +28,17 @@ transportMITgcm = nothing
 
 output_filename = ""
 
+if use_threads
+    method_barrier = PB.reaction_method_thread_barrier(
+        PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
+        PALEOmodel.ThreadBarriers.wait_barrier
+    )
+else
+    method_barrier = nothing
+end
+
 pickup_output = nothing
-initial_state, modeldata = PALEOmodel.initialize!(model, threadsafe=use_threads, pickup_output=pickup_output)  
+initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output)  
 pickup_output = nothing
 
 paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
@@ -70,7 +80,7 @@ else
     @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep)
 end
 
-isempty(output_filename) || PALEOmodel.OutputWriters.save_jld2(paleorun.output, output_filename)
+isempty(output_filename) || PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename)
 
 
 show(PB.show_variables(paleorun.model), allrows=true)

diff --git a/examples/mitgcm/config_mitgcm_expts.jl b/examples/mitgcm/config_mitgcm_expts.jl
@@ -8,38 +8,6 @@ using Printf
 include("../atmreservoirreaction.jl") # temporary solution to make ReactionReservoirAtm available
 
 
-function config_mitgcm_expts(baseconfig, expt, extrapars=Dict())
-
-
-    if baseconfig == "abiotic_O2"
-
-        model = PB.create_model_from_config(
-            joinpath(@__DIR__, "MITgcm_2deg8_abiotic.yaml"), "abiotic_O2", modelpars=extrapars)
-
-    elseif baseconfig == "PO4MMbase"
-
-            model = PB.create_model_from_config(
-                joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase", modelpars=extrapars)
-
-    elseif baseconfig == "PO4MMbaseECCO"
-
-            model = PB.create_model_from_config(
-                joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase", modelpars=extrapars)
-
-    elseif baseconfig == "PO4MMcarbSCH4"
-            model = PB.create_model_from_config(
-                    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4", modelpars=extrapars)
-
-    else
-        error("unrecognized baseconfig='$(baseconfig)'")
-    end
-
-    return model
-end
-
-
-
-
 function build_outfilename(outfileroot, segment_number)
     return outfileroot*@sprintf("_%04i", segment_number)
 end

diff --git a/examples/mitgcm/runtests.jl b/examples/mitgcm/runtests.jl
@@ -19,26 +19,28 @@ skipped_testsets = [
 
     include("config_mitgcm_expts.jl")
 
-    model = config_mitgcm_expts("PO4MMbase", "")
+    model = PB.create_model_from_config(
+        joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase"
+    )
     toutputs = [0.0, 0.25, 0.5, 0.75, 1.0]
 
     transportMITgcm = PB.get_reaction(model, "ocean", "transportMITgcm")
     tstep_imp = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr
 
     initial_state, modeldata = PALEOmodel.initialize!(model)
 
-    run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
+    paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
 
     @info "using tstep=$tstep_imp yr"
-    @time PALEOmodel.ODEfixed.integrateEuler(run, initial_state, modeldata, toutputs, tstep_imp)
+    @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_imp)
 
     println("conservation checks:")
     conschecks = [
         ("global", "total_O2",          1e-6),
         ("global",  "total_P",          1e-6),
     ]    
     for (domname, varname, rtol) in conschecks
-        startval, endval = PB.get_data(run.output, domname*"."*varname)[[1, end]]
+        startval, endval = PB.get_data(paleorun.output, domname*"."*varname)[[1, end]]
         println("  check $domname.$varname $startval $endval $rtol")
         @test isapprox(startval, endval, rtol=rtol)
     end
@@ -48,7 +50,7 @@ skipped_testsets = [
         ("ocean", "Prod_Corg_total",  3.6975e15,    1e-4),
     ]    
     for (domname, varname, checkval, rtol) in checkvals
-        outputval = PB.get_data(run.output, domname*"."*varname)[end]
+        outputval = PB.get_data(paleorun.output, domname*"."*varname)[end]
         println("  check $domname.$varname $outputval $checkval $rtol")
         @test isapprox(outputval, checkval, rtol=rtol)
     end