Skip to content

Commit 40735f4

Browse files
Test with dtmin 0
I'm running the test that was suggested in SciML/Sundials.jl#423 . > I just tried the same test case using some other solvers (this is a small stiff system of 19 equations), all at reltol=1e-4: How did you do that? I didn't see a loop over the solvers in the CI here at all. I'd like to see if this option is all that's needed.
1 parent 95ad37c commit 40735f4

File tree

4 files changed

+49
-44
lines changed

4 files changed

+49
-44
lines changed

examples/COPSE/COPSE_bergman2004_bergman2004.jl

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ include("compare_output.jl")
1717
include("copse_bergman2004_bergman2004_expts.jl")
1818

1919

20-
# load archived model output
20+
# load archived model output
2121
comparemodel = CompareOutput.copse_output_load("bergman2004","")
2222
# comparemodel = nothing
2323

@@ -42,24 +42,26 @@ if use_TEMP_DAE
4242
println("integrate, DAE")
4343
# first run includes JIT time
4444
@time PALEOmodel.ODE.integrateDAE(
45-
run, initial_state, modeldata, (-600e6, 0),
45+
run, initial_state, modeldata, (-600e6, 0),
4646
solvekwargs=(
47+
dtmin=0.0,
4748
saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
4849
)
4950
)
50-
51-
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
51+
52+
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
5253
# @time PALEOmodel.ODE.integrateDAEForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=IDA(linear_solver=:KLU))
5354
else
5455
println("integrate, ODE")
55-
# first run includes JIT time
56+
# first run includes JIT time
5657
@time PALEOmodel.ODE.integrate(
57-
run, initial_state, modeldata, (-600e6, 0),
58+
run, initial_state, modeldata, (-600e6, 0),
5859
solvekwargs=(
60+
dtmin=0.0,
5961
saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
6062
)
61-
)
62-
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
63+
)
64+
# sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
6365
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=CVODE_BDF(linear_solver=:KLU))
6466
end
6567

examples/COPSE/COPSE_reloaded_reloaded.ipynb

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -100,8 +100,8 @@
100100
"initial_deriv = similar(initial_state)\n",
101101
"\n",
102102
"PALEOmodel.SolverFunctions.ModelODE(modeldata)(\n",
103-
" initial_deriv, \n",
104-
" initial_state, \n",
103+
" initial_deriv,\n",
104+
" initial_state,\n",
105105
" nothing,\n",
106106
" 0.0\n",
107107
")\n",
@@ -139,12 +139,11 @@
139139
"\n",
140140
"# first run includes JIT time\n",
141141
"@time PALEOmodel.ODE.integrate(\n",
142-
" paleorun, initial_state, modeldata, (-1000e6, 0), \n",
142+
" paleorun, initial_state, modeldata, (-1000e6, 0),\n",
143143
" solvekwargs=(\n",
144144
" reltol=1e-4,\n",
145145
" )\n",
146-
"); \n",
147-
" "
146+
");\n"
148147
]
149148
},
150149
{
@@ -303,7 +302,7 @@
303302
"metadata": {},
304303
"outputs": [],
305304
"source": [
306-
"# model-model plot \n",
305+
"# model-model plot\n",
307306
"copse_reloaded_reloaded_plot([output_baseline_6C, output_baseline_3C]);\n"
308307
]
309308
},
@@ -379,7 +378,7 @@
379378
"println(\"state and sms vars:\")\n",
380379
"println()\n",
381380
"\n",
382-
"for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv)) \n",
381+
"for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv))\n",
383382
" 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",
384383
" 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",
385384
"end\n",

examples/COPSE/COPSE_reloaded_reloaded.jl

Lines changed: 11 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
using Logging
22
import DataFrames
33

4-
using Plots
4+
using Plots
55

66
import PALEOboxes as PB
77
import PALEOmodel
@@ -19,7 +19,7 @@ global_logger(ConsoleLogger(stderr,Logging.Info))
1919

2020
@info "Start $(@__FILE__)"
2121

22-
# load archived model output
22+
# load archived model output
2323
# comparemodel=nothing
2424

2525
# Baseline Phanerozoic configuration
@@ -43,13 +43,13 @@ tspan=(-1000e6, 0)
4343
# OOE oscillations with VCI and linear mocb
4444
# comparemodel = nothing
4545
# model = copse_reloaded_reloaded_expts(
46-
# "reloaded",
46+
# "reloaded",
4747
# [
4848
# "VCI",
4949
# "mocbProdLinear",
5050
# ("set_par", "global", "force_LIPs", "co2releasefield", "CO2max"),
5151
# ],
52-
# )
52+
# )
5353
# tspan=(-1000e6, 1e6)
5454

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

7171
@time PALEOmodel.ODE.integrate(
72-
run, initial_state, modeldata, tspan,
72+
run, initial_state, modeldata, tspan,
7373
solvekwargs=( # https://diffeq.sciml.ai/stable/basics/common_solver_opts/
74+
dtmin=0.0,
7475
reltol=1e-4,
7576
# reltol=1e-5,
7677
# saveat=1e6,
7778
# dtmax=1e4,
7879
)
79-
)
80+
)
8081

8182
# println("integrate, jacobian from autodifferentiation")
8283
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-5,))
83-
84+
8485
# PB.TestUtils.bench_model(run.model)
8586

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

101102
#####################################
102-
# Check diff against comparison model
103+
# Check diff against comparison model
103104
########################################
104105

105-
if !isnothing(comparemodel)
106+
if !isnothing(comparemodel)
106107
diff = CompareOutput.compare_copse_output(run.output, comparemodel)
107108
firstpoint=50
108109
show(sort(DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean), :std), allrows=true)
109110

110111
println()
111112

112113
# Overlay plots for comparison
113-
CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager)
114+
CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager)
114115
end
115116

116117
####################################################################################################

examples/COPSE/runtests.jl

Lines changed: 22 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ include("compare_output.jl")
1313

1414
@testset "COPSE examples" begin
1515
skipped_testsets = [
16-
# "COPSE_reloaded_reloaded",
16+
# "COPSE_reloaded_reloaded",
1717
# "COPSE_reloaded_bergman2004",
1818
# "COPSE_bergman2004_bergman2004",
1919
# "COPSE_reloaded_reloaded.ipynb",
@@ -24,7 +24,7 @@ skipped_testsets = [
2424

2525
include("copse_reloaded_reloaded_expts.jl")
2626

27-
# load archived model output
27+
# load archived model output
2828
comparemodel = CompareOutput.copse_output_load("reloaded", "reloaded_baseline")
2929

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

3535
@time PALEOmodel.ODE.integrate(
36-
run, initial_state, modeldata, (-1000e6, 0),
36+
run, initial_state, modeldata, (-1000e6, 0),
3737
solvekwargs=(
38+
dtmin=0.0,
3839
reltol=1e-4,
3940
)
40-
)
41+
)
4142
# @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-4,))
4243

4344
# conservation checks
44-
conschecks = [
45+
conschecks = [
4546
("total_C", :v, 1e-6 ),
4647
("total_C", :v_moldelta, 1e-6),
4748
("total_S", :v, 1e-6),
4849
("total_S", :v_moldelta, 1e-6),
4950
("total_redox", nothing, 1e-6)
50-
]
51+
]
5152
for (varname, propertyname, rtol) in conschecks
5253
startval, endval = PB.get_property(
53-
PB.get_data(run.output, "global."*varname),
54+
PB.get_data(run.output, "global."*varname),
5455
propertyname=propertyname,
5556
)[[1, end]]
5657
println("check $varname $startval $endval $rtol")
@@ -62,7 +63,7 @@ skipped_testsets = [
6263
firstpoint=50
6364
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
6465

65-
stdlimits = [
66+
stdlimits = [
6667
(:mccb_delta_diff, 1.5e-3),
6768
(:Sr_delta_diff, 4.0e-7),
6869
(:pCO2PAL_reldiff, 6.0e-4),
@@ -80,7 +81,7 @@ skipped_testsets = [
8081

8182
end
8283

83-
!("COPSE_reloaded_bergman2004" in skipped_testsets) &&
84+
!("COPSE_reloaded_bergman2004" in skipped_testsets) &&
8485
@testset "COPSE_reloaded_bergman2004" begin
8586
include("copse_reloaded_bergman2004_expts.jl")
8687

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

9596
@time PALEOmodel.ODE.integrate(
96-
run, initial_state, modeldata, (-1000e6, 0),
97+
run, initial_state, modeldata, (-1000e6, 0),
9798
solvekwargs=(
9899
reltol=1e-4,
100+
dtmin=0.0,
99101
# saveat=1e6,
100102
),
101103
)
102104

103105
# conservation checks
104-
conschecks = [
106+
conschecks = [
105107
("total_C", :v, 1e-6 ),
106108
("total_C", :v_moldelta, 1e-6),
107109
("total_S", :v, 1e-6),
108110
("total_S", :v_moldelta, 1e-6),
109111
("total_redox", nothing, 1e-6)
110-
]
112+
]
111113
for (varname, propertyname, rtol) in conschecks
112114
startval, endval = PB.get_property(
113115
PB.get_data(run.output, "global."*varname),
@@ -122,7 +124,7 @@ end
122124
firstpoint = 100
123125
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
124126

125-
stdlimits = [
127+
stdlimits = [
126128
(:mccb_delta_diff, 4.0e-3),
127129
(:pCO2PAL_reldiff, 3.0e-4),
128130
(:pO2PAL_reldiff, 4.0e-4),
@@ -139,7 +141,7 @@ end
139141

140142
end
141143

142-
!("COPSE_bergman2004_bergman2004" in skipped_testsets) &&
144+
!("COPSE_bergman2004_bergman2004" in skipped_testsets) &&
143145
@testset "COPSE_bergman2004_bergman2004" begin
144146
include("copse_bergman2004_bergman2004_expts.jl")
145147

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

152154
run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
153-
155+
154156
@time PALEOmodel.ODE.integrate(
155157
run, initial_state, modeldata, (-600e6, 0),
156158
solvekwargs=(
159+
dtmin=0.0,
157160
saveat=1e6,
158161
)
159162
)
160163

161164
# conservation checks
162-
conschecks = [
165+
conschecks = [
163166
("total_C", :v, 1e-6 ),
164167
("total_C", :v_moldelta, 1e-6),
165168
("total_S", :v, 1e-6),
166169
("total_S", :v_moldelta, 1e-6),
167170
("total_redox", nothing, 1e-6)
168-
]
171+
]
169172
for (varname, propertyname, rtol) in conschecks
170173
startval, endval = PB.get_property(
171174
PB.get_data(run.output, "global."*varname),
@@ -180,7 +183,7 @@ end
180183
firstpoint = 1
181184
diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
182185

183-
stdlimits = [
186+
stdlimits = [
184187
(:mccb_delta_diff, 3.0e-2),
185188
(:pCO2PAL_reldiff, 9.0e-3),
186189
(:pO2PAL_reldiff, 6.0e-3),
@@ -197,7 +200,7 @@ end
197200

198201
end
199202

200-
!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) &&
203+
!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) &&
201204
@testset "COPSE_reloaded_reloaded.ipynb" begin
202205
# unicodeplots()
203206
gr() # headless environment will require ENV["GKSwstype"] = "100"

0 commit comments

Comments
 (0)