diff --git a/.gitignore b/.gitignore index eb7dec6..d218c86 100644 --- 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 index 6405e9b..c591372 100644 --- 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 index dbbc3e8..0f904b5 100644 --- 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 diff --git a/examples/mitgcm/MITgcm_2deg8_PO4MMbase.yaml b/examples/mitgcm/MITgcm_2deg8_PO4MMbase.yaml new file mode 100644 index 0000000..f9cce76 --- /dev/null +++ b/examples/mitgcm/MITgcm_2deg8_PO4MMbase.yaml @@ -0,0 +1,283 @@ +# MITgcm 2.8deg P, O atmosphere-ocean test configuration +PO4MMbase: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + threadsafe: false + transport_pack_chunk_width: 0 # 0 (no SIMD optimization), 4, 8 + domains: + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2_total, -138.0*ocean.DOP_total] # DOP O2eq = 106.0 + 2 * 16.0 + variable_links: + sum: total_O2 + + total_P: + class: ReactionSum + parameters: + vars_to_add: [ocean.P_total, ocean.DOP_total] + variable_links: + sum: total_P + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxOceanfloor: + reactions: + particulatefluxtarget: + class: ReactionFluxTarget + operatorID: [1, 2] + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate + + solutefluxtarget: + class: ReactionFluxTarget + operatorID: [1, 2] + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["P", "O2"] #, "DIC::CIsotope", "TAlk", "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute + + atm: + + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + ocean: + reactions: + force_temperature: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_2.8deg/GCM/Theta_gcm.mat + data_var: Tgcm + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + constant_offset: 273.15 # convert C -> Kelvin + variable_links: + F: temp + + force_salinity: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_2.8deg/GCM/Salt_gcm.mat + data_var: Sgcm + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + variable_links: + F: sal + + reservoir_P: + class: ReactionReservoirTotal + operatorID: [1, 2] + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + reservoir_O2: + class: ReactionReservoirTotal + operatorID: [1, 2] + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + reservoir_DOP: + class: ReactionReservoirTotal + operatorID: [1, 2] + variable_links: + R*: DOP* + variable_attributes: + R:initial_value: 0.0e-3 # concentration m-3 + + + transportMITgcm: + class: ReactionOceanTransportTMM + operatorID: [1, 2] + parameters: + base_path: $TMMDir$/MITgcm_2.8deg + Aimp_deltat: 86400 # s 24 x 3600 + TMfpsize: 64 + pack_chunk_width: external%transport_pack_chunk_width + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.02 + + bioprod: + class: ReactionBioProdMMPop + operatorID: [2] + parameters: + depthlimit: -119.0 # first two layers only (top of layer 3 is -120.0m) + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.67 + + k_poptype: Constant + k_uPO4: 2.0e-3 # mol m-3 yr-1 + + k_nuttype: PO4MM + k_KPO4: 0.5e-3 # mol m-3 0.5 uM + + k_lightlim: MM + k_Ic: 30.0 + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + domprod_P: DOP_sms + + dopdecay: + class: ReactionParticleDecay + operatorID: [2] + parameters: + decay_timescale: 0.5 # yr + variable_links: + Particle*: DOP* + decayflux: DOP_decay + dopdecaycomponents: + class: ReactionFluxToComponents + operatorID: [2] + parameters: + outputflux_prefix: remin_ + outputflux_names: ["Corg", "N", "P"] + outputflux_stoich: [106.0, 16.0, 1.0] # must match bioprod stoich + variable_links: + inputflux: DOP_decay + + biopump: + class: ReactionExportDirectColumn + operatorID: [2] + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: Martin + martin_rovera: 0.858 + martin_depthmin: 120.0 # base of second layer + + reminocean: + class: ReactionReminO2 + operatorID: [2] + parameters: + + variable_links: + soluteflux_*: "*_sms" + + oceansurface: + reactions: + force_par: + class: ReactionForceInsolationModernEarth + + variable_links: + insolation: surface_insol + + force_wind_speed: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_2.8deg/BiogeochemData/wind_speed.mat + data_var: windspeed + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + variable_links: + F: wind_speed + + + force_open_area_fraction: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_2.8deg/BiogeochemData/ice_fraction.mat + data_var: Fice + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + scale: -1.0 # convert sea ice fraction to open area (0 - 1) + constant_offset: 1.0 + variable_links: + F: open_area_fraction + + + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston_fixed: false + + transfer_fluxAtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: ocean.oceansurface.$fluxname$_sms + variable_links: + # output_CO2: ocean.oceansurface.DIC_sms + oceanfloor: + reactions: + reminoceanfloor: + class: ReactionReminO2 + operatorID: [2] + parameters: + + variable_links: + remin*: particulateflux* + soluteflux_*: fluxOceanfloor.soluteflux_* + + transfer_particulatefluxOceanfloor: + class: ReactionFluxTransfer + operatorID: [1, 2] + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + variable_links: + + transfer_solutefluxOceanfloor: + class: ReactionFluxTransfer + operatorID: [1, 2] + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + variable_links: + diff --git a/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl b/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl index bcf448b..50fd9af 100644 --- a/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl +++ b/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.jl @@ -1,76 +1,113 @@ using Logging - -using Plots; plotlyjs(size=(750, 565)) +import LoggingExtras +using Plots import PALEOboxes as PB import PALEOmodel import PALEOocean global_logger(ConsoleLogger(stderr,Logging.Info)) + include("config_mitgcm_expts.jl") include("plot_mitgcm.jl") use_threads = true +use_split = false + +# 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 + +pickup_output = nothing + +# output_filename = "MITgcm_PO4MMcarbSCH42deg8_2000yr_20241222" +output_filename = "" + +if !isempty(output_filename) + logfile_name = output_filename*"_log.txt" + println("writing log output to ", logfile_name); flush(stdout) + logfile = open(logfile_name, "w") + global_logger(LoggingExtras.MinLevelLogger(LoggingExtras.FileLogger(logfile, always_flush=true), Logging.Info)) +else + logfile = nothing +end model = PB.create_model_from_config( - joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMcarbSCH4"; - modelpars=Dict("threadsafe"=>use_threads), + joinpath(@__DIR__, "MITgcm_2deg8_PO4MMcarbSCH4.yaml"), "PO4MMcarbSCH4"; + modelpars=Dict( + "threadsafe"=>use_threads, + "carbchem_simd_width"=>"FP32P8", # use SIMD optimisation for carbonate chemistry + "transport_pack_chunk_width" => 4, # use SIMD optimization for transport matrix + # "transport_pack_chunk_width" => 8, # use SIMD optimization for transport matrix + ), ) -toutputs = [0, 1.0, 10.0] # , 100.0, 1000.0, 1999.5, 2000.0, 2999.5, 3000.0] +toutputs = [0, 1.0, 10.0] # short run for testing # toutputs = [0, 1.0, 10.0, 100.0, 1000.0, 1999.5, 2000.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) +# 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 = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr -@info "using tstep=$tstep yr" +PB.setvalue!(transportMITgcm.pars.Aimp_deltat, tstep_implicit_s) +transportMITgcm = nothing # holds large transport matrix arrays ! -# output_filename = "MITgcm_PO4MMcarbSCH42deg8FP64_3000yr_20210210" -# 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 + PALEOmodel.ThreadBarriers.wait_barrier; + operatorID = [1], # if use_split = true, only operatorID 1 (explict transport matrix) has dependency between tiles ) else method_barrier = nothing end -pickup_output = nothing + + initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output) pickup_output = nothing paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) if !use_threads - @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep) + if use_split + @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_explicit_yr, n_inner, + cellranges_outer=cellrange_outer, + cellranges_inner=cellrange_inner + ) + else + @info "using tstep=$tstep_explicit_yr (yr)" + @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_explicit_yr) + end else - # if Threads.nthreads() == 4 - # # 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 - # elseif Threads.nthreads() == 8 - # # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - # tiles = [(1:23, :, :), (24:44, :, :), (45:61, :, :), (62:72, :, :), (73:83, :, :), (84:98, :, :), (99:116, :, :) , (117:128, :, :)] # 8 threads - # elseif Threads.nthreads() == 16 - # # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - # tiles = [(1:12, :, :), (13:22, :, :), (23:32, :, :), (33:44, :, :), - # (45:55, :, :), (56:61, :, :), (62:67, :, :), (68:72, :, :), - # (73:77, :, :), (78:83, :, :), (84:90, :, :), (91:98, :, :), - # (99:109, :, :) , (110:116, :, :) , (117:121, :, :), (122:128, :, :),] # 16 threads - - # else - # error("no config for nthreads=$(Threads.nthreads())") - # end - - # 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") - - @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep) + if use_split + @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_explicit_yr, n_inner; + cellranges_outer=cellranges_outer, cellranges_inner=cellranges_inner + ) + else + @info "using tstep=$tstep_explicit_yr (yr)" + cellranges = PB.Grids.get_tiled_cellranges(paleorun.model, Threads.nthreads(), "ocean") + @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep_explicit_yr) + end end isempty(output_filename) || PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename) @@ -111,3 +148,7 @@ pager(:newpage) +if !isnothing(logfile) + close(logfile); logfile = nothing + global_logger(ConsoleLogger(stderr, Logging.Info)) +end \ No newline at end of file diff --git a/examples/mitgcm/MITgcm_2deg8_COPDOM.yaml b/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.yaml similarity index 61% rename from examples/mitgcm/MITgcm_2deg8_COPDOM.yaml rename to examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.yaml index 06b62b4..8fce7a7 100644 --- a/examples/mitgcm/MITgcm_2deg8_COPDOM.yaml +++ b/examples/mitgcm/MITgcm_2deg8_PO4MMcarbSCH4.yaml @@ -1,284 +1,3 @@ -# MITgcm 2.8deg P, O atmosphere-ocean test configuration -PO4MMbase: - parameters: - CIsotope: ScalarData - SIsotope: ScalarData - threadsafe: false - domains: - global: - # scalar domain - - reactions: - total_O2: - class: ReactionSum - parameters: - vars_to_add: [atm.O2, ocean.O2_total, -138.0*ocean.DOP_total] # DOP O2eq = 106.0 + 2 * 16.0 - variable_links: - sum: total_O2 - - total_P: - class: ReactionSum - parameters: - vars_to_add: [ocean.P_total, ocean.DOP_total] - variable_links: - sum: total_P - - fluxAtmtoOceansurface: - reactions: - fluxtarget: - class: ReactionFluxTarget - - parameters: - flux_totals: true - fluxlist: ["O2"] - - fluxOceanfloor: - reactions: - particulatefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: particulateflux_ - fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate - - solutefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: soluteflux_ - fluxlist: ["P", "O2"] #, "DIC::CIsotope", "TAlk", "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute - - atm: - - - reactions: - reservoir_O2: - class: ReactionReservoirAtm - - variable_links: - R*: O2* - pRatm: pO2atm - pRnorm: pO2PAL - variable_attributes: - R:norm_value: 3.7e19 # present-day atmospheric level - R:initial_value: 3.7e19 - - transfer_AtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Distribute - transfer_multiplier: -1.0 - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: $fluxname$_sms - - ocean: - reactions: - force_temperature: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_2.8deg/GCM/Theta_gcm.mat - data_var: Tgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - constant_offset: 273.15 # convert C -> Kelvin - variable_links: - F: temp - - force_salinity: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_2.8deg/GCM/Salt_gcm.mat - data_var: Sgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: sal - - reservoir_P: - class: ReactionReservoirTotal - operatorID: [1, 2] - variable_links: - R*: P* - variable_attributes: - R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) - - reservoir_O2: - class: ReactionReservoirTotal - operatorID: [1, 2] - variable_links: - R*: O2* - variable_attributes: - R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) - - reservoir_DOP: - class: ReactionReservoirTotal - operatorID: [1, 2] - variable_links: - R*: DOP* - variable_attributes: - R:initial_value: 0.0e-3 # concentration m-3 - - - transportMITgcm: - class: ReactionOceanTransportTMM - operatorID: [1, 2] - parameters: - base_path: $TMMDir$/MITgcm_2.8deg - TMfpsize: 64 - pack_chunk_width: 4 - - light: - class: ReactionLightColumn - parameters: - background_opacity: 0.02 - - bioprod: - class: ReactionBioProdMMPop - operatorID: [2] - parameters: - depthlimit: -119.0 # first two layers only (top of layer 3 is -120.0m) - - rCorgPO4: 106.0 - rNPO4: 16.0 - rCcarbCorg: 0.25 - rCcarbCorg_fixed: true - - nuDOM: 0.67 - - k_poptype: Constant - k_uPO4: 2.0e-3 # mol m-3 yr-1 - - k_nuttype: PO4MM - k_KPO4: 0.5e-3 # mol m-3 0.5 uM - - k_lightlim: MM - k_Ic: 30.0 - k_Irel: 1.0 - - variable_links: - partprod_*: export_* - domprod_P: DOP_sms - - dopdecay: - class: ReactionParticleDecay - operatorID: [2] - parameters: - decay_timescale: 0.5 # yr - variable_links: - Particle*: DOP* - decayflux: DOP_decay - dopdecaycomponents: - class: ReactionFluxToComponents - operatorID: [2] - parameters: - outputflux_prefix: remin_ - outputflux_names: ["Corg", "N", "P"] - outputflux_stoich: [106.0, 16.0, 1.0] # must match bioprod stoich - variable_links: - inputflux: DOP_decay - - biopump: - class: ReactionExportDirectColumn - operatorID: [2] - parameters: - fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] - transportfloor: true - exportfunction: Martin - martin_rovera: 0.858 - martin_depthmin: 120.0 # base of second layer - - reminocean: - class: ReactionReminO2 - operatorID: [2] - parameters: - - variable_links: - soluteflux_*: "*_sms" - - oceansurface: - reactions: - force_par: - class: ReactionForceInsolationModernEarth - - variable_links: - insolation: surface_insol - - force_wind_speed: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_2.8deg/BiogeochemData/wind_speed.mat - data_var: windspeed - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: wind_speed - - - force_open_area_fraction: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_2.8deg/BiogeochemData/ice_fraction.mat - data_var: Fice - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - scale: -1.0 # convert sea ice fraction to open area (0 - 1) - constant_offset: 1.0 - variable_links: - F: open_area_fraction - - - airsea_O2: - class: ReactionAirSeaO2 - parameters: - piston_fixed: false - - transfer_fluxAtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Identity - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: ocean.oceansurface.$fluxname$_sms - variable_links: - # output_CO2: ocean.oceansurface.DIC_sms - oceanfloor: - reactions: - reminoceanfloor: - class: ReactionReminO2 - operatorID: [2] - parameters: - - variable_links: - remin*: particulateflux* - soluteflux_*: fluxOceanfloor.soluteflux_* - - transfer_particulatefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ - output_fluxes: particulateflux_$fluxname$ - variable_links: - - transfer_solutefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ - output_fluxes: ocean.oceanfloor.$fluxname$_sms - variable_links: - ############################################################ # MITgcm 2.8deg P, O with Cinorg, S, CH4 ############################################################ @@ -287,6 +6,8 @@ PO4MMcarbSCH4: CIsotope: IsotopeLinear # ScalarData SIsotope: IsotopeLinear threadsafe: false + carbchem_simd_width: "1" # "1" (no SIMD optimization) FP32P8 FP64P4 + transport_pack_chunk_width: 0 # 0 (no SIMD optimization), 4, 8 domains: global: @@ -415,14 +136,16 @@ PO4MMcarbSCH4: reservoir_P: - class: ReactionReservoirTotal + class: ReactionReservoirTotal + operatorID: [1, 2] variable_links: R*: P* variable_attributes: R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) reservoir_O2: - class: ReactionReservoirTotal + class: ReactionReservoirTotal + operatorID: [1, 2] variable_links: R*: O2* variable_attributes: @@ -430,6 +153,7 @@ PO4MMcarbSCH4: reservoir_DIC: class: ReactionReservoirTotal + operatorID: [1, 2] parameters: field_data: external%CIsotope variable_links: @@ -441,6 +165,7 @@ PO4MMcarbSCH4: reservoir_DOC: class: ReactionReservoirTotal + operatorID: [1, 2] parameters: field_data: external%CIsotope variable_links: @@ -451,7 +176,8 @@ PO4MMcarbSCH4: R:norm_value: 1.0e-3 # for scaling only reservoir_TAlk: - class: ReactionReservoirTotal + class: ReactionReservoirTotal + operatorID: [1, 2] variable_links: R*: TAlk* variable_attributes: @@ -460,6 +186,7 @@ PO4MMcarbSCH4: reservoir_SO4: class: ReactionReservoirTotal + operatorID: [1, 2] parameters: field_data: external%SIsotope variable_links: @@ -471,6 +198,7 @@ PO4MMcarbSCH4: reservoir_H2S: class: ReactionReservoirTotal + operatorID: [1, 2] parameters: field_data: external%SIsotope limit_delta_conc: 1e-6 # mol m-3 to limit delta **EXPERIMENTAL** @@ -483,6 +211,7 @@ PO4MMcarbSCH4: reservoir_CH4: class: ReactionReservoirTotal + operatorID: [1, 2] parameters: field_data: external%CIsotope limit_delta_conc: 1e-6 # mol m-3 to limit delta **EXPERIMENTAL** @@ -495,12 +224,13 @@ PO4MMcarbSCH4: carbchem: class: ReactionCO2SYS + operatorID: [1, 2] parameters: components: ["Ci", "B", "S", "F", "Omega", "H2S"] defaultconcs: ["TF", "TB", "Ca"] solve_pH: solve # solve for Hfree iteratively outputs: ["pCO2", "xCO2dryinp", "CO2", "CO3", "OmegaCA", "OmegaAR"] - simd_width: FP32P8 # FP64P4 + simd_width: external%carbchem_simd_width # SIMD optimization pHtol: 1.2e-5 # 100*eps(Float32) variable_links: TCi_conc: DIC_conc @@ -512,10 +242,12 @@ PO4MMcarbSCH4: transportMITgcm: class: ReactionOceanTransportTMM + operatorID: [1, 2] parameters: base_path: $TMMDir$/MITgcm_2.8deg + Aimp_deltat: 86400 # s 24 x 3600 TMfpsize: 64 # 32 - pack_chunk_width: 4 # 8 + pack_chunk_width: external%transport_pack_chunk_width light: class: ReactionLightColumn @@ -532,6 +264,7 @@ PO4MMcarbSCH4: bioprod: class: ReactionBioProdMMPop + operatorID: [2] parameters: depthlimit: -60.0 # surface cells only @@ -560,6 +293,7 @@ PO4MMcarbSCH4: docdecay: class: ReactionParticleDecay + operatorID: [2] parameters: decay_timescale: 0.5 # yr field_data: external%CIsotope @@ -568,6 +302,7 @@ PO4MMcarbSCH4: decayflux: DOC_decay docdecaycomponents: class: ReactionFluxToComponents + operatorID: [2] parameters: outputflux_prefix: remin_ outputflux_names: ["Corg::Isotope", "N", "P"] @@ -578,7 +313,7 @@ PO4MMcarbSCH4: biopumporg: class: ReactionExportDirectColumn - + operatorID: [2] parameters: fluxlist: ["P", "N", "Corg::CIsotope"] transportfloor: true @@ -588,7 +323,7 @@ PO4MMcarbSCH4: biopumpcarb: class: ReactionExportDirectColumn - + operatorID: [2] parameters: fluxlist: ["Ccarb::CIsotope"] transportfloor: true @@ -598,7 +333,7 @@ PO4MMcarbSCH4: reminocean: class: ReactionReminO2_SO4_CH4 - + operatorID: [2] parameters: SO4reminlimit: 1000.0e-3 @@ -607,13 +342,13 @@ PO4MMcarbSCH4: redox_H2S_O2: class: ReactionRedoxH2S_O2 - + operatorID: [2] parameters: R_H2S_O2: 3.65e2 # 3.65e3 # (mol m-3) yr-1 redox_CH4_O2: class: ReactionRedoxCH4_O2 - + operatorID: [2] parameters: R_CH4_O2: 1.0e2 # 10.0e3 # (mol m-3) yr-1 @@ -681,7 +416,7 @@ PO4MMcarbSCH4: reactions: reminoceanfloor: class: ReactionReminO2_SO4_CH4 - + operatorID: [2] parameters: SO4reminlimit: 1000.0e-3 diff --git a/examples/mitgcm/MITgcm_2deg8_abiotic.jl b/examples/mitgcm/MITgcm_2deg8_abiotic.jl index 66d09cd..216f77d 100644 --- a/examples/mitgcm/MITgcm_2deg8_abiotic.jl +++ b/examples/mitgcm/MITgcm_2deg8_abiotic.jl @@ -15,22 +15,23 @@ use_threads = false model = PB.create_model_from_config( joinpath(@__DIR__, "MITgcm_2deg8_abiotic.yaml"), "abiotic_O2"; - modelpars=Dict("threadsafe"=>use_threads), + modelpars=Dict( + "threadsafe"=>use_threads, + "transport_pack_chunk_width" => 4, # use SIMD optimization for transport matrix + ), ) -toutputs_relative = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0] +toutputs = [0, 0.25, 0.5, 0.75, 1.0] #, 10.0] +tstep_explicit_s::Int = 86400 # s timestep to use for explicit transport +tstep_implicit_s::Int = tstep_explicit_s +tstep_explicit_yr = tstep_explicit_s/PB.Constants.k_secpyr # yr +@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 = transportMITgcm.pars.Aimp_deltat[]/PB.Constants.k_secpyr +PB.setvalue!(transportMITgcm.pars.Aimp_deltat, tstep_implicit_s) +transportMITgcm = nothing # holds large transport matrix arrays ! -@info "using tstep=$tstep yr" - -num_segments = 1 -outfile_root = "" -# num_segments = 20 -# outfile_root = "MITgcm_PO4MMbase2deg8_100yr_20210201" - -toutputs = [] +output_filename = "" if use_threads method_barrier = PB.reaction_method_thread_barrier( @@ -41,51 +42,23 @@ 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_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename) - - tstart = PB.get_data(pickup_output, "ocean.tmodel")[end] - - else - pickup_output = nothing - tstart = 0.0 - end - - global initial_state - global modeldata - global toutputs - - initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output) - pickup_output = nothing - - toutputs = toutputs_relative .+ tstart - - global paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) - - if !use_threads - @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep) - 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(run.model, Threads.nthreads(), "ocean") - - @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep) - end - - if !isempty(outfile_root) - output_filename = build_outfilename(outfile_root, iseg) - PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename) - end +initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier) + +paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) + +if !use_threads + @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep_explicit_yr) +else + cellranges = PB.Grids.get_tiled_cellranges(run.model, Threads.nthreads(), "ocean") + + @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep_explicit_yr) +end + +if !isempty(output_filename) + PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename) end + show(PB.show_variables(paleorun.model), allrows=true) println() diff --git a/examples/mitgcm/MITgcm_2deg8_abiotic.yaml b/examples/mitgcm/MITgcm_2deg8_abiotic.yaml index c04b031..638565f 100644 --- a/examples/mitgcm/MITgcm_2deg8_abiotic.yaml +++ b/examples/mitgcm/MITgcm_2deg8_abiotic.yaml @@ -2,7 +2,8 @@ abiotic_O2: parameters: CIsotope: ScalarData - threadsafe: false + threadsafe: false + transport_pack_chunk_width: 0 # 0 (no SIMD optimization), 4, 8 domains: global: # scalar domain @@ -95,7 +96,7 @@ abiotic_O2: parameters: base_path: $TMMDir$/MITgcm_2.8deg TMfpsize: 64 - pack_chunk_width: 0 # 4 + pack_chunk_width: external%transport_pack_chunk_width light: class: ReactionLightColumn diff --git a/examples/mitgcm/MITgcm_2deg8_test.jl b/examples/mitgcm/MITgcm_2deg8_test.jl deleted file mode 100644 index 3864d12..0000000 --- a/examples/mitgcm/MITgcm_2deg8_test.jl +++ /dev/null @@ -1,85 +0,0 @@ -using Logging - -using Plots; plotlyjs(size=(750, 565)) - -import PALEOboxes as PB -import PALEOmodel -import PALEOocean - -global_logger(ConsoleLogger(stderr,Logging.Info)) -include("config_mitgcm_expts.jl") - -use_threads = true - -# 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 = 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) - - -transportMITgcm = PB.get_reaction(model, "ocean", "transportMITgcm") -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; method_barrier) - -run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) - -if !use_threads - @time PALEOmodel.ODEfixed.integrateEuler(run, initial_state, modeldata, toutputs , tstep) -else - #= - if Threads.nthreads() == 4 - # 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 - elseif Threads.nthreads() == 8 - # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - tiles = [(1:23, :, :), (24:44, :, :), (45:61, :, :), (62:72, :, :), (73:83, :, :), (84:98, :, :), (99:116, :, :) , (117:128, :, :)] # 8 threads - elseif Threads.nthreads() == 16 - # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - tiles = [(1:12, :, :), (13:22, :, :), (23:32, :, :), (33:44, :, :), - (45:55, :, :), (56:61, :, :), (62:67, :, :), (68:72, :, :), - (73:77, :, :), (78:83, :, :), (84:90, :, :), (91:98, :, :), - (99:109, :, :) , (110:116, :, :) , (117:121, :, :), (122:128, :, :),] # 16 threads - - else - error("no config for nthreads=$(Threads.nthreads())") - end - - cellranges = PB.Grids.get_tiled_cellranges(run.model, tiles) # vector of vectors, 1 per tile -=# - cellranges = PB.Grids.get_tiled_cellranges(run.model, Threads.nthreads(), "ocean") - - @time PALEOmodel.ODEfixed.integrateEulerthreads(run, initial_state, modeldata, cellranges, toutputs , tstep) -end - - -show(PB.show_variables(run.model), allrows=true) -println() - diff --git a/examples/mitgcm/MITgcm_ECCO_COPDOM.yaml b/examples/mitgcm/MITgcm_ECCO_COPDOM.yaml deleted file mode 100644 index a1cf6c9..0000000 --- a/examples/mitgcm/MITgcm_ECCO_COPDOM.yaml +++ /dev/null @@ -1,644 +0,0 @@ -######################################################### -# MITgcm ECCO P, O atmosphere-ocean test configuration -####################################################### -PO4MMbase: - parameters: - CIsotope: ScalarData - SIsotope: ScalarData - threadsafe: false - domains: - global: - # scalar domain - - reactions: - total_O2: - class: ReactionSum - parameters: - vars_to_add: [atm.O2, ocean.O2_total] - variable_links: - sum: total_O2 - - fluxAtmtoOceansurface: - reactions: - fluxtarget: - class: ReactionFluxTarget - - parameters: - flux_totals: true - fluxlist: ["O2"] - - fluxOceanfloor: - reactions: - particulatefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: particulateflux_ - fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate - - solutefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: soluteflux_ - fluxlist: ["P", "O2"] #, "DIC::CIsotope", "TAlk", "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute - - atm: - - - reactions: - reservoir_O2: - class: ReactionReservoirAtm - - variable_links: - R*: O2* - pRatm: pO2atm - pRnorm: pO2PAL - variable_attributes: - R:norm_value: 3.7e19 # present-day atmospheric level - R:initial_value: 3.7e19 - - transfer_AtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Distribute - transfer_multiplier: -1.0 - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: $fluxname$_sms - - ocean: - reactions: - force_temperature: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Theta_gcm.mat - data_var: Tgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - constant_offset: 273.15 # convert C -> Kelvin - variable_links: - F: temp - - force_salinity: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Salt_gcm.mat - data_var: Sgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: sal - - reservoir_P: - class: ReactionReservoirTotal - variable_links: - R*: P* - variable_attributes: - R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) - - reservoir_O2: - class: ReactionReservoirTotal - variable_links: - R*: O2* - variable_attributes: - R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) - - - - transportMITgcm: - class: ReactionOceanTransportTMM - parameters: - base_path: $TMMDir$/MITgcm_ECCO - Aimp_deltat: 43200 # s 12 x 3600 - TMfpsize: 32 - pack_chunk_width: 4 # SIMD transport - - - light: - class: ReactionLightColumn - parameters: - background_opacity: 0.02 - - bioprod: - class: ReactionBioProdMMPop - parameters: - depthlimit: -119.0 # equivalent of top two layers in 2.8 deg grid (top of layer 3 in 2.8deg grid is -120.0m) - - rCorgPO4: 106.0 - rNPO4: 16.0 - rCcarbCorg: 0.25 - rCcarbCorg_fixed: true - - nuDOM: 0.66 - - k_poptype: Constant - k_uPO4: 3.0e-3 # mol m-3 yr-1 3 uM yr-1 - - k_nuttype: PO4MM - k_KPO4: 0.5e-3 # mol m-3 0.5 uM - - k_lightlim: MM - k_Ic: 30.0 - k_Irel: 1.0 - - variable_links: - partprod_*: export_* - domprod_*: remin_* # no explicit DOM pool, so reroute directly to remin - - biopump: - class: ReactionExportDirectColumn - - parameters: - fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] - transportfloor: true - exportfunction: Martin - martin_rovera: 0.858 - martin_depthmin: 120.0 # base of second layer - - reminocean: - class: ReactionReminO2 - - parameters: - - variable_links: - soluteflux_*: "*_sms" - - oceansurface: - reactions: - force_par: - class: ReactionForceInsolationModernEarth - - variable_links: - insolation: surface_insol - - force_wind_speed: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/wind_speed.mat - data_var: windspeed - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: wind_speed - - - force_open_area_fraction: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/ice_fraction.mat - data_var: Fice - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - scale: -1.0 # convert sea ice fraction to open area (0 - 1) - constant_offset: 1.0 - variable_links: - F: open_area_fraction - - - airsea_O2: - class: ReactionAirSeaO2 - parameters: - piston_fixed: false - - transfer_fluxAtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Identity - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: ocean.oceansurface.$fluxname$_sms - variable_links: - # output_CO2: ocean.oceansurface.DIC_sms - - oceanfloor: - reactions: - reminoceanfloor: - class: ReactionReminO2 - - parameters: - - variable_links: - remin*: particulateflux* - soluteflux_*: fluxOceanfloor.soluteflux_* - - transfer_particulatefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ - output_fluxes: particulateflux_$fluxname$ - variable_links: - - transfer_solutefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ - output_fluxes: ocean.oceanfloor.$fluxname$_sms - variable_links: - -######################################### -# MITgcm ECCO P, O with Cinorg, S, CH4 -########################################## -PO4MMcarbSCH4: - parameters: - CIsotope: IsotopeLinear # ScalarData - SIsotope: IsotopeLinear - - domains: - global: - # scalar domain - - reactions: - total_O2eq: - class: ReactionSum - parameters: - vars_to_add: [atm.O2, ocean.O2_total, -2*ocean.CH4_total, -2*ocean.H2S_total] - component_to_add: 1 # we just want the first component (total) from isotope variables - variable_links: - sum: total_O2eq - - total_C: - class: ReactionSum - parameters: - vars_to_add: [atm.CO2, ocean.DIC_total, ocean.CH4_total] - variable_links: - sum: total_C - - total_S: - class: ReactionSum - parameters: - vars_to_add: [ocean.SO4_total, ocean.H2S_total] - variable_links: - sum: total_S - - fluxAtmtoOceansurface: - reactions: - fluxtarget: - class: ReactionFluxTarget - - parameters: - flux_totals: true - fluxlist: ["O2", "CO2::CIsotope"] - - fluxOceanfloor: - reactions: - particulatefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: particulateflux_ - fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate - - solutefluxtarget: - class: ReactionFluxTarget - operatorID: [1, 2] - parameters: - flux_totals: true - target_prefix: soluteflux_ - fluxlist: ["P", "O2", "DIC::CIsotope", "TAlk", "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute - - atm: - - - reactions: - reservoir_O2: - class: ReactionReservoirAtm - - variable_links: - R*: O2* - pRatm: pO2atm - pRnorm: pO2PAL - variable_attributes: - R:norm_value: 3.7e19 # present-day atmospheric level - R:initial_value: 3.7e19 - - reservoir_CO2: - class: ReactionReservoirAtm - parameters: - field_data: external%CIsotope - - variable_links: - R*: CO2* - pRatm: pCO2atm - pRnorm: pCO2PAL - variable_attributes: - R:norm_value: 4.956e16 # pre ind 280e-6 - R:initial_value: 4.956e16 # pre ind 280e-6 - - transfer_AtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Distribute - transfer_multiplier: -1.0 - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: $fluxname$_sms - - ocean: - reactions: - force_temperature: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Theta_gcm.mat - data_var: Tgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - constant_offset: 273.15 # convert C -> Kelvin - variable_links: - F: temp - - force_salinity: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Salt_gcm.mat - data_var: Sgcm - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: sal - - - reservoir_P: - class: ReactionReservoirTotal - variable_links: - R*: P* - variable_attributes: - R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) - - - reservoir_O2: - class: ReactionReservoirTotal - variable_links: - R*: O2* - variable_attributes: - R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) - - - reservoir_DIC: - class: ReactionReservoirTotal - parameters: - field_data: external%CIsotope - variable_links: - R*: DIC* - variable_attributes: - R:initial_value: 2291.74e-3 # 2231.486e-6 mol kg-1 * 1027 kg m-3 - R:initial_delta: -1.0 - R:norm_value: 1000.0e-3 # for scaling only - - reservoir_TAlk: - class: ReactionReservoirTotal - variable_links: - R*: TAlk* - variable_attributes: - R:initial_value: 2426.8e-3 # 2363e-6 mol kg-1 * 1027 kg m-3 - R:norm_value: 1000.0e-3 # for scaling only - - reservoir_SO4: - class: ReactionReservoirTotal - parameters: - field_data: external%SIsotope - variable_links: - R*: SO4* - variable_attributes: - R:initial_value: 28756.0e-3 # concentration mol m-3 ~ 28e-3 mol/kg * 1027 kg m-3 - R:initial_delta: 1.0 - R:norm_value: 1000.0e-3 # for scaling only - - reservoir_H2S: - class: ReactionReservoirTotal - parameters: - field_data: external%SIsotope - limit_delta_conc: 1e-6 # mol m-3 to limit delta **EXPERIMENTAL** - variable_links: - R*: H2S* - variable_attributes: - R:initial_value: 1e-6 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 - R:initial_delta: 0.0 - R:norm_value: 1.0e-3 # for scaling only - - reservoir_CH4: - class: ReactionReservoirTotal - parameters: - field_data: external%CIsotope - limit_delta_conc: 1e-6 # mol m-3 to limit delta **EXPERIMENTAL** - variable_links: - R*: CH4* - variable_attributes: - R:initial_value: 1e-6 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 - R:initial_delta: 0.0 - R:norm_value: 1.0e-3 # for scaling only - - carbchem: - class: ReactionCO2SYS - parameters: - components: ["Ci", "B", "S", "F", "Omega", "H2S"] - defaultconcs: ["TF", "TB", "Ca"] - solve_pH: solve # solve for Hfree iteratively - outputs: ["pCO2", "xCO2dryinp", "CO2", "CO3", "OmegaCA", "OmegaAR"] - simd_width: FP32P8 # FP64P4 - pHtol: 1.2e-5 # 100*eps(Float32) - variable_links: - TCi_conc: DIC_conc - TS_conc: SO4_conc - CO2: CO2_conc - pCO2: pCO2 - OmegaCA: OmegaCA - TH2S_conc: H2S_conc - - transportMITgcm: - class: ReactionOceanTransportTMM - parameters: - base_path: $TMMDir$/MITgcm_ECCO - - light: - class: ReactionLightColumn - parameters: - background_opacity: 0.02 - - const_cisotopes: - class: ReactionScalarConst - parameters: - constnames: ["D_mccb_DIC", "D_B_mccb_mocb"] - variable_attributes: - D_mccb_DIC:initial_value: 0.0 - D_B_mccb_mocb:initial_value: 25.0 - - bioprod: - class: ReactionBioProdMMPop - parameters: - depthlimit: -60.0 # surface cells only - - rCorgPO4: 106.0 - rNPO4: 16.0 - # rCcarbCorg: 0.25 - rCcarbCorg_fixed: false - k_r0: 0.0485 - k_eta: 0.7440 - - nuDOM: 0.66 - - k_poptype: Constant - k_uPO4: 3.0e-3 # mol m-3 yr-1 3 uM yr-1 - - k_nuttype: PO4MM - k_KPO4: 0.5e-3 # mol m-3 0.5 uM - - k_lightlim: MM - k_Ic: 30.0 - k_Irel: 1.0 - - variable_links: - partprod_*: export_* - domprod_*: remin_* # no explicit DOM pool, so reroute directly to remin - - biopumporg: - class: ReactionExportDirectColumn - - parameters: - fluxlist: ["P", "N", "Corg::CIsotope"] - transportfloor: true - exportfunction: Martin - martin_rovera: 0.858 - martin_depthmin: 120.0 # base of second layer - - biopumpcarb: - class: ReactionExportDirectColumn - - parameters: - fluxlist: ["Ccarb::CIsotope"] - transportfloor: true - exportfunction: SumExp - input_frac: [0.55, 0.45] # 2-G model of Ridgwell & Hargreaves (2007) - sumexp_scale: [1890.5, 1e6] - - reminocean: - class: ReactionReminO2_SO4_CH4 - - parameters: - SO4reminlimit: 1000.0e-3 - - variable_links: - soluteflux_*: "*_sms" - - redox_H2S_O2: - class: ReactionRedoxH2S_O2 - - parameters: - R_H2S_O2: 3.65e2 # 3.65e3 # (mol m-3) yr-1 - - redox_CH4_O2: - class: ReactionRedoxCH4_O2 - - parameters: - R_CH4_O2: 1.0e2 # 10.0e3 # (mol m-3) yr-1 - - oceansurface: - reactions: - force_par: - class: ReactionForceInsolationModernEarth - - variable_links: - insolation: surface_insol - - force_wind_speed: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/wind_speed.mat - data_var: windspeed - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - variable_links: - F: wind_speed - - - force_open_area_fraction: - class: ReactionForceGrid - parameters: - matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/ice_fraction.mat - data_var: Fice - time_var: "" - tidx_start: 1 - tidx_end: 12 - cycle_time: 1.0 - scale: -1.0 # convert sea ice fraction to open area (0 - 1) - constant_offset: 1.0 - variable_links: - F: open_area_fraction - - airsea_O2: - class: ReactionAirSeaO2 - parameters: - piston_fixed: false - - airsea_CO2: - class: ReactionAirSeaCO2 - parameters: - - piston_fixed: false - moistair: false # no sat H2O correction - - variable_links: - Xatm_delta: atm.CO2_delta - Xocean_delta: ocean.oceansurface.DIC_delta - - transfer_fluxAtmtoOceansurface: - class: ReactionFluxTransfer - parameters: - transfer_matrix: Identity - input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ - output_fluxes: ocean.oceansurface.$fluxname$_sms - variable_links: - output_CO2: ocean.oceansurface.DIC_sms - - oceanfloor: - reactions: - reminoceanfloor: - class: ReactionReminO2_SO4_CH4 - - parameters: - SO4reminlimit: 1000.0e-3 - - variable_links: - remin*: particulateflux* - O2_conc: ocean.oceanfloor.O2_conc - SO4_*: ocean.oceanfloor.SO4_* # conc and delta - - soluteflux_*: fluxOceanfloor.soluteflux_* - - transfer_particulatefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ - output_fluxes: particulateflux_$fluxname$ - variable_links: - - transfer_solutefluxOceanfloor: - class: ReactionFluxTransfer - operatorID: [1, 2] - parameters: - transfer_matrix: Identity - input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ - output_fluxes: ocean.oceanfloor.$fluxname$_sms - variable_links: - - diff --git a/examples/mitgcm/MITgcm_ECCO_abiotic.jl b/examples/mitgcm/MITgcm_ECCO_PO4MMbase.jl similarity index 50% rename from examples/mitgcm/MITgcm_ECCO_abiotic.jl rename to examples/mitgcm/MITgcm_ECCO_PO4MMbase.jl index dfec64e..78139b6 100644 --- a/examples/mitgcm/MITgcm_ECCO_abiotic.jl +++ b/examples/mitgcm/MITgcm_ECCO_PO4MMbase.jl @@ -1,5 +1,5 @@ using Logging - +import LoggingExtras using Plots import PALEOboxes as PB @@ -11,22 +11,36 @@ include("config_mitgcm_expts.jl") include("plot_mitgcm.jl") use_threads = true -do_benchmarks = false + +outfile_root = "MITgcm_PO4MMbaseECCO_100yr_20241220" + +logfile_name = outfile_root*"_1_log.txt" +println("writing log output to ", logfile_name); flush(stdout) +logfile = open(logfile_name, "w") +global_logger(LoggingExtras.MinLevelLogger(LoggingExtras.FileLogger(logfile, always_flush=true), Logging.Info)) model = PB.create_model_from_config( - joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase"; - modelpars=Dict("threadsafe"=>use_threads), + joinpath(@__DIR__, "MITgcm_ECCO_PO4MMbase.yaml"), "PO4MMbase"; + modelpars=Dict( + "threadsafe"=>use_threads, + "transport_pack_chunk_width" => 4, # use SIMD optimization for transport matrix + ), ) toutputs_relative = [0.0, 1.0, 10.0, 100.0] +# toutputs_relative = [0.0, 1.0] + +tstep_explicit_s::Int = 43200 # s timestep to use for explicit transport +tstep_implicit_s::Int = tstep_explicit_s +tstep_explicit_yr = tstep_explicit_s/PB.Constants.k_secpyr # yr +@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 = transportMITgcm.par_Aimp_deltat[]/PB.Constants.k_secpyr -transportMITgcm = nothing +PB.setvalue!(transportMITgcm.pars.Aimp_deltat, tstep_implicit_s) +transportMITgcm = nothing # holds large transport matrix arrays ! -@info "using tstep=$tstep yr" num_segments = 20 -outfile_root = "MITgcm_PO4MMbaseECCO_100yr_20210130" +# num_segments = 2 toutputs=[] @@ -39,40 +53,53 @@ else method_barrier = nothing end +initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier) + for iseg in 1:num_segments -# for iseg in 13:num_segments +# for iseg in 9:num_segments +# iseg = 1 +# if true + @info """ + + ============================================================================== + start iseg $iseg + =============================================================================== + """ + if iseg > 1 pickup_filename = build_outfilename(outfile_root, iseg-1) pickup_output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), pickup_filename) + + PALEOmodel.set_statevar_from_output!(modeldata, pickup_output) + global initial_state = PALEOmodel.get_statevar(modeldata.solver_view_all) tstart = PB.get_data(pickup_output, "ocean.tmodel")[end] + pickup_output = nothing else - pickup_output = nothing tstart = 0.0 end - initial_state, modeldata = PALEOmodel.initialize!(model; method_barrier, pickup_output) - pickup_output = nothing toutputs = toutputs_relative .+ tstart + @info "toutputs: $toutputs" global paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) if !use_threads - @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep) + @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep_explicit_yr) 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 - # tiles = [(1:124, :, :), (125:202, :, :), (203:275, :, :), (276:360, :, :)] # 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") - @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep) + @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep_explicit_yr) end output_filename = build_outfilename(outfile_root, iseg) + @info """ + ============================================================================== + end iseg $iseg + toutputs $toutputs + output_filename $output_filename + =============================================================================== + + """ PALEOmodel.OutputWriters.save_netcdf(paleorun.output, output_filename) end @@ -105,3 +132,7 @@ plot_PO4MMbase( pager(:newpage) +if !isnothing(logfile) + close(logfile); logfile = nothing + global_logger(ConsoleLogger(stderr, Logging.Info)) +end diff --git a/examples/mitgcm/MITgcm_ECCO_PO4MMbase.yaml b/examples/mitgcm/MITgcm_ECCO_PO4MMbase.yaml new file mode 100644 index 0000000..f789bce --- /dev/null +++ b/examples/mitgcm/MITgcm_ECCO_PO4MMbase.yaml @@ -0,0 +1,251 @@ +######################################################### +# MITgcm ECCO P, O atmosphere-ocean test configuration +####################################################### +PO4MMbase: + parameters: + CIsotope: ScalarData + SIsotope: ScalarData + threadsafe: false + transport_pack_chunk_width: 0 # 0 (no SIMD optimization), 4, 8 + domains: + global: + # scalar domain + + reactions: + total_O2: + class: ReactionSum + parameters: + vars_to_add: [atm.O2, ocean.O2_total] + variable_links: + sum: total_O2 + + fluxAtmtoOceansurface: + reactions: + fluxtarget: + class: ReactionFluxTarget + + parameters: + flux_totals: true + fluxlist: ["O2"] + + fluxOceanfloor: + reactions: + particulatefluxtarget: + class: ReactionFluxTarget + operatorID: [1, 2] + parameters: + flux_totals: true + target_prefix: particulateflux_ + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] # fluxlist_BioParticulate + + solutefluxtarget: + class: ReactionFluxTarget + operatorID: [1, 2] + parameters: + flux_totals: true + target_prefix: soluteflux_ + fluxlist: ["P", "O2"] #, "DIC::CIsotope", "TAlk", "SO4::SIsotope", "H2S::SIsotope", "CH4::CIsotope"] # fluxlist_Solute + + atm: + + + reactions: + reservoir_O2: + class: ReactionReservoirAtm + + variable_links: + R*: O2* + pRatm: pO2atm + pRnorm: pO2PAL + variable_attributes: + R:norm_value: 3.7e19 # present-day atmospheric level + R:initial_value: 3.7e19 + + transfer_AtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Distribute + transfer_multiplier: -1.0 + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: $fluxname$_sms + + ocean: + reactions: + force_temperature: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Theta_gcm.mat + data_var: Tgcm + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + constant_offset: 273.15 # convert C -> Kelvin + variable_links: + F: temp + + force_salinity: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_ECCO/GCM/Salt_gcm.mat + data_var: Sgcm + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + variable_links: + F: sal + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + variable_attributes: + R:initial_value: 2.208e-3 # concentration m-3 (1027 kg m-3 * 2.15e-6 mol/kg-sw) + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + variable_attributes: + R:initial_value: 0.2054 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + + + + transportMITgcm: + class: ReactionOceanTransportTMM + parameters: + base_path: $TMMDir$/MITgcm_ECCO + Aimp_deltat: 43200 # s 12 x 3600 + # TMfpsize: 32 + TMfpsize: 64 + pack_chunk_width: external%transport_pack_chunk_width + + light: + class: ReactionLightColumn + parameters: + background_opacity: 0.02 + + bioprod: + class: ReactionBioProdMMPop + parameters: + depthlimit: -119.0 # equivalent of top two layers in 2.8 deg grid (top of layer 3 in 2.8deg grid is -120.0m) + + rCorgPO4: 106.0 + rNPO4: 16.0 + rCcarbCorg: 0.25 + rCcarbCorg_fixed: true + + nuDOM: 0.66 + + k_poptype: Constant + k_uPO4: 3.0e-3 # mol m-3 yr-1 3 uM yr-1 + + k_nuttype: PO4MM + k_KPO4: 0.5e-3 # mol m-3 0.5 uM + + k_lightlim: MM + k_Ic: 30.0 + k_Irel: 1.0 + + variable_links: + partprod_*: export_* + domprod_*: remin_* # no explicit DOM pool, so reroute directly to remin + + biopump: + class: ReactionExportDirectColumn + + parameters: + fluxlist: ["P", "N", "Corg::CIsotope", "Ccarb::CIsotope"] + transportfloor: true + exportfunction: Martin + martin_rovera: 0.858 + martin_depthmin: 120.0 # base of second layer + + reminocean: + class: ReactionReminO2 + + parameters: + + variable_links: + soluteflux_*: "*_sms" + + oceansurface: + reactions: + force_par: + class: ReactionForceInsolationModernEarth + + variable_links: + insolation: surface_insol + + force_wind_speed: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/wind_speed.mat + data_var: windspeed + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + variable_links: + F: wind_speed + + + force_open_area_fraction: + class: ReactionForceGrid + parameters: + matlab_file: $TMMDir$/MITgcm_ECCO/BiogeochemData/ice_fraction.mat + data_var: Fice + time_var: "" + tidx_start: 1 + tidx_end: 12 + cycle_time: 1.0 + scale: -1.0 # convert sea ice fraction to open area (0 - 1) + constant_offset: 1.0 + variable_links: + F: open_area_fraction + + + airsea_O2: + class: ReactionAirSeaO2 + parameters: + piston_fixed: false + + transfer_fluxAtmtoOceansurface: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$ + output_fluxes: ocean.oceansurface.$fluxname$_sms + variable_links: + # output_CO2: ocean.oceansurface.DIC_sms + + oceanfloor: + reactions: + reminoceanfloor: + class: ReactionReminO2 + + parameters: + + variable_links: + remin*: particulateflux* + soluteflux_*: fluxOceanfloor.soluteflux_* + + transfer_particulatefluxOceanfloor: + class: ReactionFluxTransfer + operatorID: [1, 2] + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: particulateflux_$fluxname$ + variable_links: + + transfer_solutefluxOceanfloor: + class: ReactionFluxTransfer + operatorID: [1, 2] + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.soluteflux_$fluxname$ + output_fluxes: ocean.oceanfloor.$fluxname$_sms + variable_links: + diff --git a/examples/mitgcm/MITgcm_ECCO_abiotic_test.jl b/examples/mitgcm/MITgcm_ECCO_abiotic_test.jl deleted file mode 100644 index d81b3b3..0000000 --- a/examples/mitgcm/MITgcm_ECCO_abiotic_test.jl +++ /dev/null @@ -1,90 +0,0 @@ -using Logging - -using Plots - -import PALEOboxes as PB -import PALEOmodel -import PALEOocean - -global_logger(ConsoleLogger(stderr,Logging.Info)) -include("config_mitgcm_expts.jl") -include("plot_mitgcm.jl") - -use_threads = true -do_benchmarks = false - -model = PB.create_model_from_config( - joinpath(@__DIR__, "MITgcm_ECCO_COPDOM.yaml"), "PO4MMbase"; - modelpars=Dict("threadsafe"=>use_threads), -) - -toutputs = [0.0, 1.0] # , 10.0, 100.0] - -transportMITgcm = PB.get_reaction(model, "ocean", "transportMITgcm") -tstep = transportMITgcm.par_Aimp_deltat[]/PB.Constants.k_secpyr -transportMITgcm = nothing - -@info "using tstep=$tstep yr" - -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; method_barrier, pickup_output) -pickup_output = nothing - -paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) - -if !use_threads - @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs , tstep) -else - # if Threads.nthreads() == 4 - # # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - # tiles = [(1:124, :, :), (125:202, :, :), (203:275, :, :), (276:360, :, :)] # 4 threads - # elseif Threads.nthreads() == 8 - # # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - # tiles = [(1:62, :, :), (63:120, :, :), - # (121:167, :, :), (168:200, :, :), - # (201:232, :, :), (233:275, :, :), - # (276:325, :, :) , (326:360, :, :)] # 8 threads - # elseif Threads.nthreads() == 16 - # # indices are slightly uneven to equalize the number of active (mostly ocean) cells per thread - # tiles = [(1:31, :, :), (32:62, :, :), (63:88, :, :), (89:120, :, :), - # (121:148, :, :), (149:167, :, :), (168:184, :, :), (185:199, :, :), - # (200:215, :, :), (216:230, :, :), (231:250, :, :), (251:274, :, :), - # (275:304, :, :), (305:325, :, :), (326:340, :, :), (341:360, :, :),] # 16 threads - - # else - # error("no config for nthreads=$(Threads.nthreads())") - # end - - # 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") - - for (t,tc) in enumerate(cellranges) - numcells = 0 - for c in tc - numcells += length(c.indices) - end - println("thread $t numcells $numcells") - end - @time PALEOmodel.ODEfixed.integrateEulerthreads(paleorun, initial_state, modeldata, cellranges, toutputs , tstep) -end - -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/README.md b/examples/mitgcm/README.md index e930315..3432267 100644 --- a/examples/mitgcm/README.md +++ b/examples/mitgcm/README.md @@ -16,7 +16,7 @@ Tracers: atmosphere O2, ocean O2 julia> include("MITgcm_2deg8_abiotic.jl") -Wallclock time: 10.3 s core-1 (model yr)-1 with default timestep = 86400 s (1 day) +Wallclock time: 2.3 s (model yr)-1 using 1 core, with default timestep = 86400 s (1 day) ## 2.8 degree P, O2 @@ -26,7 +26,7 @@ Tracers: atmosphere O2, ocean O2, P, DOP julia> include("MITgcm_2deg8_PO4MMbase.jl") -Wallclock time: 10.0 s core-1 (model yr)-1 with default timestep = 86400 s (1 day) +Wallclock time: 1.4 s (model yr)-1 using 4 cores, with default timestep = 86400 s (1 day) ## 2.8 degree P, O2, S, DIC/TAlk @@ -38,7 +38,7 @@ Tracers: atmosphere O2, CO2(x2), ocean O2, P, DOC(x2), H2S(x2), SO4(x2), CH4(x2) julia> include("MITgcm_2deg8_PO4MMcarbSCH4.jl") -Wallclock time: 27.2 s core-1 (model yr)-1 with default timestep = 86400 s (1 day) +Wallclock time: 2.8 s (model yr)-1 using 4 cores, with default timestep = 86400 s (1 day) ![Surface P](images/surface_P_2deg8_PO4MMcarbSCH4_2000yr.svg) @@ -50,3 +50,12 @@ Wallclock time: 27.2 s core-1 (model yr)-1 with default timestep = 86400 s (1 d ###### Figure 2 *H2S concentration at 2000yr, sections at longitudes corresponding to red lines in Figure 1* +## ECCO (~1 degree) P, O2 + +Minimal test case for biotic ocean. + +Tracers: atmosphere O2, ocean O2, P + + julia> include("MITgcm_ECCO_PO4MMbase.jl") + +Wallclock time: 50.0 s (model yr)-1 using 4 cores, with default timestep = 43200 s (1 day) \ No newline at end of file diff --git a/examples/mitgcm/plot_mitgcm.jl b/examples/mitgcm/plot_mitgcm.jl index b404fb9..fb8f340 100644 --- a/examples/mitgcm/plot_mitgcm.jl +++ b/examples/mitgcm/plot_mitgcm.jl @@ -71,9 +71,15 @@ function plot_PO4MMbase( ) end - pager( - plot(title="P total", output, ["global.total_P"], ylabel="total (mol)",), - ) + if PB.has_variable(output, "global.total_P") + pager( + plot(title="P total", output, ["global.total_P"], ylabel="total (mol)",), + ) + elseif PB.has_variable(output, "ocean.P_total") + pager( + plot(title="P total", output, ["ocean.P_total"], ylabel="total (mol)",), + ) + end for t in tbioprod pager( diff --git a/examples/mitgcm/runtests.jl b/examples/mitgcm/runtests.jl index 2ffe0ac..4c16975 100644 --- a/examples/mitgcm/runtests.jl +++ b/examples/mitgcm/runtests.jl @@ -20,19 +20,24 @@ skipped_testsets = [ include("config_mitgcm_expts.jl") model = PB.create_model_from_config( - joinpath(@__DIR__, "MITgcm_2deg8_COPDOM.yaml"), "PO4MMbase" + joinpath(@__DIR__, "MITgcm_2deg8_PO4MMbase.yaml"), "PO4MMbase" ) toutputs = [0.0, 0.25, 0.5, 0.75, 1.0] + tstep_explicit_s::Int = 86400 # s timestep to use for explicit transport + tstep_implicit_s::Int = tstep_explicit_s + tstep_explicit_yr = tstep_explicit_s/PB.Constants.k_secpyr # yr + @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 ! + initial_state, modeldata = PALEOmodel.initialize!(model) paleorun = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory()) - @info "using tstep=$tstep_imp yr" - @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_imp) + @time PALEOmodel.ODEfixed.integrateEuler(paleorun, initial_state, modeldata, toutputs, tstep_explicit_yr) println("conservation checks:") conschecks = [ diff --git a/src/ocean/OceanTransportMatrix.jl b/src/ocean/OceanTransportMatrix.jl index 9681ca5..b78b671 100644 --- a/src/ocean/OceanTransportMatrix.jl +++ b/src/ocean/OceanTransportMatrix.jl @@ -406,7 +406,9 @@ struct PackedBuffer{pack_eltype, pack_datatype, pack_chunk_width, num_components packed_conc_array = zeros(pack_eltype, packed_conc_pad_width, ncells) packed_conc_vec = vec(packed_conc_array) - buffer = [Vector{pack_datatype}(undef, num_chunks) for t in 1:Threads.nthreads()] + # buffer = [Vector{pack_datatype}(undef, num_chunks) for t in 1:Threads.nthreads()] + # julia 1.11.2 workaround: enforce 64-byte alignment by creating large arrays + buffer = [Vector{pack_datatype}(undef, num_chunks + 256) for t in 1:Threads.nthreads()] return new{pack_eltype, pack_datatype, pack_chunk_width, num_components, num_chunks, packed_conc_pad_width}( packed_conc_array, @@ -480,8 +482,13 @@ function _kernel_transport_tr_packed( tm = (nzval1[idx]*wt1 + nzval2[idx]*wt2)*volume[jcol] irow = rowval[idx] conc_vec_idx = packed_conc_pad_width*(irow-1) + 1 - for cidx in 1:num_chunks - pack_sms_buffer[cidx] += SIMD.vload(pack_datatype, packed_conc_vec, conc_vec_idx)*tm + for cidx in 1:num_chunks + # slow on julia 1.11 + # pack_sms_buffer[cidx] += SIMD.vload(pack_datatype, packed_conc_vec, conc_vec_idx)*tm + # needs buffer alignment fix on julia 1.11.2 + x = SIMD.vload(pack_datatype, packed_conc_vec, conc_vec_idx)*tm + y = pack_sms_buffer[cidx] + pack_sms_buffer[cidx] = x + y conc_vec_idx += pack_chunk_width end end diff --git a/src/ocean/OceanTransportTMM.jl b/src/ocean/OceanTransportTMM.jl index 4967b84..fd617de 100644 --- a/src/ocean/OceanTransportTMM.jl +++ b/src/ocean/OceanTransportTMM.jl @@ -173,6 +173,7 @@ function _read_grids(rj::ReactionOceanTransportTMM) (rj.index_perm, rj.index_perm_inverse) = PB.Grids.linear_kji_order(rj.grid_ocean, v_i, v_j, v_k) v_i = v_i[rj.index_perm]; v_j = v_j[rj.index_perm]; v_k = v_k[rj.index_perm] end + PB.setfrozen!(rj.pars.kji_order) PB.Grids.set_linear_index(rj.grid_ocean, v_i, v_j, v_k) @info " set ocean linear <--> cartesian mapping for $(rj.grid_ocean.ncells) cells" @@ -219,6 +220,12 @@ function PB.register_methods!(rj::ReactionOceanTransportTMM) (PB.VarList_namedtuple(PALEOocean.Ocean.grid_vars_all),) ) + PB.add_method_setup!( + rj, + setup_transport_TMM, # reads transport matrices + (), + ) + return nothing end @@ -301,7 +308,7 @@ function PB.register_dynamic_methods!(rj::ReactionOceanTransportTMM) PB.VarList_components(transport_sms_vars), PB.VarList_nothing(), # not using input_vars ), - preparefn=prepare_do_transport_TMM, # read matrices, add buffer + preparefn=PALEOocean.Ocean.prepare_transport, # add buffer operatorID=[rj.operatorID[1]], p=:trspt_dAexp_tr # field name in rj to use ) @@ -342,8 +349,7 @@ function PB.register_dynamic_methods!(rj::ReactionOceanTransportTMM) ( PB.VarList_components(transport_conc_vars), PB.VarList_single(sequencer_var), - ), - preparefn=prepare_do_transport_TMM_packed, # read matrices + ), p=packed_buffer, ) # Aexp operatorID[1] @@ -373,29 +379,27 @@ function PB.register_dynamic_methods!(rj::ReactionOceanTransportTMM) end + PB.setfrozen!(rj.pars.TMfpsize, rj.pars.pack_chunk_width) + return nothing end -function prepare_do_transport_TMM(m::PB.ReactionMethod, vardata) - rj = m.reaction - - # read matrix data in prepare so it is deferred until initialize - _read_matrix_data(rj) - - return PALEOocean.Ocean.prepare_transport(m, vardata) -end +function setup_transport_TMM( + m::PB.ReactionMethod, + (), + cellrange::PB.AbstractCellRange, + attribute_name, +) + attribute_name == :setup || return -function prepare_do_transport_TMM_packed(m::PB.ReactionMethod, vardata) rj = m.reaction - # read matrix data in prepare so it is deferred until initialize _read_matrix_data(rj) - return vardata + return nothing end - function _read_matrix_data(rj::ReactionOceanTransportTMM) # calculate model time for list of matrices @@ -415,9 +419,7 @@ function _read_matrix_data(rj::ReactionOceanTransportTMM) ) rj.matrix_tinterp = PB.LinInterp(rj.matrix_times, 1.0) end - - tmp_Aexp_tr= [] - tmp_Aimp_tr = [] + PB.setfrozen!(rj.pars.use_annualmean, rj.pars.num_seasonal) # transpose, convert datatype, optionally permute function permute_indices_transpose(A) @@ -450,39 +452,27 @@ function _read_matrix_data(rj::ReactionOceanTransportTMM) ] Aimp_name = "Aimp" end - - tmp_Aexp_tr = Vector{Any}(nothing, rj.pars.num_seasonal[]) - tmp_Aimp_tr = Vector{Any}(nothing, rj.pars.num_seasonal[]) - + # work out upscaling factor for implicit matrix deltaT = rj.grid["deltaT"] rj.pars.Aimp_deltat[] % deltaT == 0 || error("requested Aimp_deltat $(rj.pars.Aimp_deltat[]) is not a multiple of matrix deltaT = $deltaT") rj.Aimp_mult = rj.pars.Aimp_deltat[]/deltaT Aimp_deltat_yr = rj.pars.Aimp_deltat[]/PB.Constants.k_secpyr + PB.setfrozen!(rj.pars.Aimp_deltat) # For large (1 deg) matrices, attempt to minimise memory use # NB: the underlying hdf5 library is not thread safe so we can't use threads without # adding a lot of complexity / using a lot of memory - # tlock = ReentrantLock() - # Threads.@threads - for i in 1:num_matrices - matAexp = nothing - matAimp = nothing - # @Base.lock tlock begin # serialize disk reads to avoid hammering the disk - Aexp_path = Aexp_paths[i] - @info " reading explicit matrix data from $Aexp_path" - file = MAT.matopen(Aexp_path) - matAexp = MAT.read(file, Aexp_name) - close(file) - - Aimp_path = Aimp_paths[i] - @info " reading implicit matrix data from $Aimp_path" - file = MAT.matopen(Aimp_path) - matAimp = MAT.read(file, Aimp_name) - close(file) - # end + tmp_Aexp_tr = Vector{Any}(nothing, rj.pars.num_seasonal[]) + for i in 1:num_matrices + Aexp_path = Aexp_paths[i] + @info " reading explicit matrix data from $Aexp_path" + file = MAT.matopen(Aexp_path) + matAexp = MAT.read(file, Aexp_name) + close(file) + tmp_Aexp_tr[i] = permute_indices_transpose(matAexp) matAexp = nothing if rj.pars.sal_norm[] @@ -491,6 +481,22 @@ function _read_matrix_data(rj::ReactionOceanTransportTMM) end # Aexp is already in differential form, convert s-1 to yr-1 tmp_Aexp_tr[i] .= PB.Constants.k_secpyr.*tmp_Aexp_tr[i] + end + # convert to CSR format with common sparsity pattern for fast multiply x vector + rj.trspt_dAexp_tr = PALEOocean.Ocean.create_common_sparsity_tr!( + tmp_Aexp_tr, + do_transpose=false, + TMeltype=rj.TMeltype + ) + tmp_Aexp_tr = [] + + tmp_Aimp_tr = Vector{Any}(nothing, rj.pars.num_seasonal[]) + for i in 1:num_matrices + Aimp_path = Aimp_paths[i] + @info " reading implicit matrix data from $Aimp_path" + file = MAT.matopen(Aimp_path) + matAimp = MAT.read(file, Aimp_name) + close(file) # convert implicit matrix tmp_Aimp_tr[i] = permute_indices_transpose(matAimp) @@ -507,15 +513,6 @@ function _read_matrix_data(rj::ReactionOceanTransportTMM) tmp_Aimp_tr[i] = (tmp_Aimp_tr[i]^rj.Aimp_mult - LinearAlgebra.I)./Aimp_deltat_yr end - - # convert to CSR format with common sparsity pattern for fast multiply x vector - rj.trspt_dAexp_tr = PALEOocean.Ocean.create_common_sparsity_tr!( - tmp_Aexp_tr, - do_transpose=false, - TMeltype=rj.TMeltype - ) - tmp_Aexp_tr = [] - # convert to CSR format with common sparsity pattern for fast multiply x vector rj.trspt_dAimp_tr = PALEOocean.Ocean.create_common_sparsity_tr!( tmp_Aimp_tr, @@ -523,10 +520,13 @@ function _read_matrix_data(rj::ReactionOceanTransportTMM) TMeltype=rj.TMeltype ) tmp_Aimp_tr = [] - + + PB.setfrozen!(rj.pars.sal_norm) + return nothing end + function do_transport_TMM( m::PB.ReactionMethod, (grid_tforce_vars, conc_components, sms_components, _, buffer), @@ -569,6 +569,7 @@ function do_transport_TMM_pack_conc( return nothing end + function do_transport_TMM_packed( m::PB.ReactionMethod, (grid_tforce_vars, sms_components, dummy_var),