Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Test with dtmin 0 #33

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Test with dtmin 0 #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 10 additions & 8 deletions examples/COPSE/COPSE_bergman2004_bergman2004.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ include("compare_output.jl")
include("copse_bergman2004_bergman2004_expts.jl")


# load archived model output
# load archived model output
comparemodel = CompareOutput.copse_output_load("bergman2004","")
# comparemodel = nothing

Expand All @@ -42,24 +42,26 @@ if use_TEMP_DAE
println("integrate, DAE")
# first run includes JIT time
@time PALEOmodel.ODE.integrateDAE(
run, initial_state, modeldata, (-600e6, 0),
run, initial_state, modeldata, (-600e6, 0),
solvekwargs=(
dtmin=0.0,
saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
)
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?

# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
# @time PALEOmodel.ODE.integrateDAEForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=IDA(linear_solver=:KLU))
else
println("integrate, ODE")
# first run includes JIT time
# first run includes JIT time
@time PALEOmodel.ODE.integrate(
run, initial_state, modeldata, (-600e6, 0),
run, initial_state, modeldata, (-600e6, 0),
solvekwargs=(
dtmin=0.0,
saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
)
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
)
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=CVODE_BDF(linear_solver=:KLU))
end

Expand Down
13 changes: 6 additions & 7 deletions examples/COPSE/COPSE_reloaded_reloaded.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,8 +100,8 @@
"initial_deriv = similar(initial_state)\n",
"\n",
"PALEOmodel.SolverFunctions.ModelODE(modeldata)(\n",
" initial_deriv, \n",
" initial_state, \n",
" initial_deriv,\n",
" initial_state,\n",
" nothing,\n",
" 0.0\n",
")\n",
Expand Down Expand Up @@ -139,12 +139,11 @@
"\n",
"# first run includes JIT time\n",
"@time PALEOmodel.ODE.integrate(\n",
" paleorun, initial_state, modeldata, (-1000e6, 0), \n",
" paleorun, initial_state, modeldata, (-1000e6, 0),\n",
" solvekwargs=(\n",
" reltol=1e-4,\n",
" )\n",
"); \n",
" "
");\n"
]
},
{
Expand Down Expand Up @@ -303,7 +302,7 @@
"metadata": {},
"outputs": [],
"source": [
"# model-model plot \n",
"# model-model plot\n",
"copse_reloaded_reloaded_plot([output_baseline_6C, output_baseline_3C]);\n"
]
},
Expand Down Expand Up @@ -379,7 +378,7 @@
"println(\"state and sms vars:\")\n",
"println()\n",
"\n",
"for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv)) \n",
"for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv))\n",
" println(var.name, \"\\t\\t\", PB.get_var_type(var), \"\\t\\t\", PB.get_attribute(var, :units), \"\\t\\t\", PB.get_attribute(var, :description))\n",
" println(var_sms.name, \"\\t\\t\", PB.get_var_type(var_sms), \"\\t\\t\", PB.get_attribute(var_sms, :units), \"\\t\\t\", PB.get_attribute(var_sms, :description))\n",
"end\n",
Expand Down
21 changes: 11 additions & 10 deletions examples/COPSE/COPSE_reloaded_reloaded.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
using Logging
import DataFrames

using Plots
using Plots

import PALEOboxes as PB
import PALEOmodel
Expand All @@ -19,7 +19,7 @@ global_logger(ConsoleLogger(stderr,Logging.Info))

@info "Start $(@__FILE__)"

# load archived model output
# load archived model output
# comparemodel=nothing

# Baseline Phanerozoic configuration
Expand All @@ -43,13 +43,13 @@ tspan=(-1000e6, 0)
# OOE oscillations with VCI and linear mocb
# comparemodel = nothing
# model = copse_reloaded_reloaded_expts(
# "reloaded",
# "reloaded",
# [
# "VCI",
# "mocbProdLinear",
# ("set_par", "global", "force_LIPs", "co2releasefield", "CO2max"),
# ],
# )
# )
# tspan=(-1000e6, 1e6)

######################################
Expand All @@ -69,18 +69,19 @@ println("integrate, no jacobian")
run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())

@time PALEOmodel.ODE.integrate(
run, initial_state, modeldata, tspan,
run, initial_state, modeldata, tspan,
solvekwargs=( # https://diffeq.sciml.ai/stable/basics/common_solver_opts/
dtmin=0.0,
reltol=1e-4,
# reltol=1e-5,
# saveat=1e6,
# dtmax=1e4,
)
)
)

# println("integrate, jacobian from autodifferentiation")
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-5,))

# PB.TestUtils.bench_model(run.model)

##############################
Expand All @@ -99,18 +100,18 @@ pager=PALEOmodel.PlotPager(9, (legend_background_color=nothing, ))
copse_reloaded_reloaded_plot(run.output, pager=pager)

#####################################
# Check diff against comparison model
# Check diff against comparison model
########################################

if !isnothing(comparemodel)
if !isnothing(comparemodel)
diff = CompareOutput.compare_copse_output(run.output, comparemodel)
firstpoint=50
show(sort(DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean), :std), allrows=true)

println()

# Overlay plots for comparison
CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager)
CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager)
end

####################################################################################################
Expand Down
41 changes: 22 additions & 19 deletions examples/COPSE/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ include("compare_output.jl")

@testset "COPSE examples" begin
skipped_testsets = [
# "COPSE_reloaded_reloaded",
# "COPSE_reloaded_reloaded",
# "COPSE_reloaded_bergman2004",
# "COPSE_bergman2004_bergman2004",
# "COPSE_reloaded_reloaded.ipynb",
Expand All @@ -24,7 +24,7 @@ skipped_testsets = [

include("copse_reloaded_reloaded_expts.jl")

# load archived model output
# load archived model output
comparemodel = CompareOutput.copse_output_load("reloaded", "reloaded_baseline")

model = copse_reloaded_reloaded_expts("reloaded", ["baseline"])
Expand All @@ -33,24 +33,25 @@ skipped_testsets = [
run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())

@time PALEOmodel.ODE.integrate(
run, initial_state, modeldata, (-1000e6, 0),
run, initial_state, modeldata, (-1000e6, 0),
solvekwargs=(
dtmin=0.0,
reltol=1e-4,
)
)
)
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-4,))

# conservation checks
conschecks = [
conschecks = [
("total_C", :v, 1e-6 ),
("total_C", :v_moldelta, 1e-6),
("total_S", :v, 1e-6),
("total_S", :v_moldelta, 1e-6),
("total_redox", nothing, 1e-6)
]
]
for (varname, propertyname, rtol) in conschecks
startval, endval = PB.get_property(
PB.get_data(run.output, "global."*varname),
PB.get_data(run.output, "global."*varname),
propertyname=propertyname,
)[[1, end]]
println("check $varname $startval $endval $rtol")
Expand All @@ -62,7 +63,7 @@ skipped_testsets = [
firstpoint=50
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)

stdlimits = [
stdlimits = [
(:mccb_delta_diff, 1.5e-3),
(:Sr_delta_diff, 4.0e-7),
(:pCO2PAL_reldiff, 6.0e-4),
Expand All @@ -80,7 +81,7 @@ skipped_testsets = [

end

!("COPSE_reloaded_bergman2004" in skipped_testsets) &&
!("COPSE_reloaded_bergman2004" in skipped_testsets) &&
@testset "COPSE_reloaded_bergman2004" begin
include("copse_reloaded_bergman2004_expts.jl")

Expand All @@ -93,21 +94,22 @@ end
run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())

@time PALEOmodel.ODE.integrate(
run, initial_state, modeldata, (-1000e6, 0),
run, initial_state, modeldata, (-1000e6, 0),
solvekwargs=(
reltol=1e-4,
dtmin=0.0,
# saveat=1e6,
),
)

# conservation checks
conschecks = [
conschecks = [
("total_C", :v, 1e-6 ),
("total_C", :v_moldelta, 1e-6),
("total_S", :v, 1e-6),
("total_S", :v_moldelta, 1e-6),
("total_redox", nothing, 1e-6)
]
]
for (varname, propertyname, rtol) in conschecks
startval, endval = PB.get_property(
PB.get_data(run.output, "global."*varname),
Expand All @@ -122,7 +124,7 @@ end
firstpoint = 100
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)

stdlimits = [
stdlimits = [
(:mccb_delta_diff, 4.0e-3),
(:pCO2PAL_reldiff, 3.0e-4),
(:pO2PAL_reldiff, 4.0e-4),
Expand All @@ -139,7 +141,7 @@ end

end

!("COPSE_bergman2004_bergman2004" in skipped_testsets) &&
!("COPSE_bergman2004_bergman2004" in skipped_testsets) &&
@testset "COPSE_bergman2004_bergman2004" begin
include("copse_bergman2004_bergman2004_expts.jl")

Expand All @@ -150,22 +152,23 @@ end
initial_state, modeldata = PALEOmodel.initialize!(model)

run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())

@time PALEOmodel.ODE.integrate(
run, initial_state, modeldata, (-600e6, 0),
solvekwargs=(
dtmin=0.0,
saveat=1e6,
)
)

# conservation checks
conschecks = [
conschecks = [
("total_C", :v, 1e-6 ),
("total_C", :v_moldelta, 1e-6),
("total_S", :v, 1e-6),
("total_S", :v_moldelta, 1e-6),
("total_redox", nothing, 1e-6)
]
]
for (varname, propertyname, rtol) in conschecks
startval, endval = PB.get_property(
PB.get_data(run.output, "global."*varname),
Expand All @@ -180,7 +183,7 @@ end
firstpoint = 1
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)

stdlimits = [
stdlimits = [
(:mccb_delta_diff, 3.0e-2),
(:pCO2PAL_reldiff, 9.0e-3),
(:pO2PAL_reldiff, 6.0e-3),
Expand All @@ -197,7 +200,7 @@ end

end

!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) &&
!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) &&
@testset "COPSE_reloaded_reloaded.ipynb" begin
# unicodeplots()
gr() # headless environment will require ENV["GKSwstype"] = "100"
Expand Down