Merge pull request #31 from PALEOtoolkit/mitgcm_updates_fixes

Tidy up and bugfix mitgcm configs
PALEOtoolkit · Dec 22, 2024 · 00baead · 00baead
2 parents bb1fbf3 + 5739753
commit 00baead
Show file tree

Hide file tree

Showing 18 changed files with 853 additions and 1,317 deletions.
diff --git a/.gitignore b/.gitignore
@@ -6,5 +6,6 @@ LocalPreferences.toml
 *.jld2
 *.nc
 *.zip
+*_log.txt
 examples/shelf1D/S2P3_transport_20240614
 examples/romglb/romaniello2010_transport
diff --git a/examples/Project.toml b/examples/Project.toml
@@ -5,6 +5,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
@@ -14,16 +15,16 @@ PALEOcopse = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
 PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
 PALEOocean = "41de04b1-2efd-44ae-92ae-39d71a4fd99b"
 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"
+ThreadPinning = "811555cd-349b-4f26-b7bc-1f208b848042"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [compat]
-PALEOmodel = "0.15.48"
+PALEOaqchem = "0.3.18"
+PALEOboxes = "0.21.40"
+PALEOmodel = "0.15.49"
 PALEOsediment = "0.3.3"
-
diff --git a/examples/mitgcm/MITgcm_2deg8_PO4MMbase.jl b/examples/mitgcm/MITgcm_2deg8_PO4MMbase.jl
@@ -1,6 +1,6 @@
 using Logging
 
-using Plots; plotlyjs(size=(750, 565))
+using Plots
 
 import PALEOboxes as PB
 import PALEOmodel
@@ -10,22 +10,36 @@ global_logger(ConsoleLogger(stderr,Logging.Info))
 include("config_mitgcm_expts.jl")
 include("plot_mitgcm.jl")
 
-use_threads = false
+use_threads = true
 use_split   = false
-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]
+# transport matrices have deltaT=1200 s (timestep used to run the ocean model)
+tstep_explicit_s::Int = 86400 # s timestep to use for explicit transport
+n_inner = 8         # if use_split, number of implicit timesteps per explicit timestep
+# n_inner = 72        # max value, gives implicit timestep 1200 s
 
 model = PB.create_model_from_config(
-    joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase";
-    modelpars=Dict("threadsafe"=>use_threads),
+    joinpath(@__DIR__, "MITgcm_2deg8_PO4MMbase.yaml"), "PO4MMbase";
+    modelpars=Dict(
+        "threadsafe"=>use_threads,
+        "transport_pack_chunk_width" => 4, # use SIMD optimization for transport matrix
+    ),
 )
 
 toutputs = [0.0, 0.25, 0.5, 0.75, 1.0, 10.0] 
 # 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]
 
+# configure timestepping
+tstep_explicit_yr = tstep_explicit_s/PB.Constants.k_secpyr # yr
+if use_split
+    tstep_implicit_s::Int = tstep_explicit_s / n_inner
+else
+    tstep_implicit_s::Int = tstep_explicit_s 
+end
+@info "using timesteps tstep_implicit_s $tstep_implicit_s tstep_explicit_s $tstep_explicit_s tstep_explicit_yr $tstep_explicit_yr"
 transportMITgcm = PB.get_reaction(model, "ocean", "transportMITgcm")
-tstep_imp = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr
+PB.setvalue!(transportMITgcm.pars.Aimp_deltat, tstep_implicit_s)
+transportMITgcm = nothing # holds large transport matrix arrays !
 
 
 output_filename = ""
@@ -34,7 +48,8 @@ output_filename = ""
 if use_threads
     method_barrier = PB.reaction_method_thread_barrier(
         PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
-        PALEOmodel.ThreadBarriers.wait_barrier
+        PALEOmodel.ThreadBarriers.wait_barrier;
+        operatorID = [1], # if use_split = true, only operatorID 1 (explict transport matrix) has dependency between tiles
     )
 else
     method_barrier = nothing
@@ -48,38 +63,33 @@ paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMem
 
 if !use_threads
     if use_split       
-        @info "using tstep_outer=$n_inner x $tstep_imp yr"
+        @info "using tstep_outer=$tstep_explicit_yr (yr), n_inner $n_inner, tstep_inner $(tstep_explicit_yr/n_inner) yr"
         cellrange_outer = PB.create_default_cellrange(paleorun.model, operatorID=1)
         cellrange_inner = PB.create_default_cellrange(paleorun.model, operatorID=2)
 
         @time PALEOmodel.ODEfixed.integrateSplitEuler(
-            paleorun, initial_state, modeldata, toutputs, tstep_imp*n_inner, n_inner,
-            cellrange_outer=cellrange_outer,
-            cellrange_inner=cellrange_inner
+            paleorun, initial_state, modeldata, toutputs, tstep_explicit_yr, n_inner,
+            cellranges_outer=cellrange_outer,
+            cellranges_inner=cellrange_inner
         )
     else
-        @info "using tstep=$tstep_imp yr"
-        @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_imp)
+        @info "using tstep=$tstep_explicit_yr (yr)"
+        @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_explicit_yr)
     end
 else
-    # Threads.nthreads() == 4 || error("use_threads requires 4 threads, Threads.nthreads()=", Threads.nthreads())
-
-    # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread
-    # tiles = [(1:44, :, :), (45:72, :, :), (73:98, :, :), (99:128, :, :)]  # 4 threads
-
-    # cellranges = PB.Grids.get_tiled_cellranges(paleorun.model, tiles)  # vector of vectors, 1 per tile
-
     cellranges = PB.Grids.get_tiled_cellranges(paleorun.model, Threads.nthreads(), "ocean")
 
     if use_split
-        @info "using tstep_outer=$n_inner x $tstep_imp yr"
+        @info "using tstep_outer=$tstep_explicit_yr (yr), n_inner $n_inner, tstep_inner $(tstep_explicit_yr/n_inner) yr"
         cellranges_outer = PB.Grids.get_tiled_cellranges(paleorun.model, Threads.nthreads(), "ocean", operatorID=1)
         cellranges_inner = PB.Grids.get_tiled_cellranges(paleorun.model, Threads.nthreads(), "ocean", operatorID=2)
-        @time PALEOmodel.ODEfixed.integrateSplitEulerthreads(paleorun, initial_state, modeldata, toutputs , tstep_imp*n_inner, n_inner,
-                cellranges_outer=cellranges_outer, cellranges_inner=cellranges_inner)
+        @time PALEOmodel.ODEfixed.integrateSplitEulerthreads(
+            paleorun, initial_state, modeldata, toutputs , tstep_explicit_yr, n_inner;
+            cellranges_outer=cellranges_outer, cellranges_inner=cellranges_inner
+        )
     else
-        @info "using tstep=$tstep_imp yr"
-        @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep_imp)
+        @info "using tstep=$tstep_explicit_yr (yr)"
+        @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep_explicit_yr)
     end
 end