Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 3 additions & 1 deletion .typos.toml
Original file line number Diff line number Diff line change
Expand Up @@ -12,4 +12,6 @@ myu = "myu"
ND = "ND"
nd = "nd"
ned = "ned"
sitl = "sitl"
sitl = "sitl"
Ser = "Ser"
inh = "inh"
161 changes: 161 additions & 0 deletions benchmarks/Bio/BioPreDynB4.jmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
---
title: BioPreDyn-bench B4 Work-Precision Diagrams
author: Chris Rackauckas
---

The following benchmark is of 35 ODEs modeling the metabolism of Chinese
Hamster Ovary (CHO) cells from the
[BioPreDyn-bench](https://doi.org/10.1186/s12918-015-0144-4) benchmark suite
(Problem B4). The model describes a batch fermentation process with 32
reactions including glycolysis, TCA cycle, amino acid metabolism, and product
protein formation. The kinetics use log-linear elasticity coefficients with
117 parameters.

The model is manually translated from the original BioPreDyn-bench Matlab/C
files since the SBML representation uses rate rules rather than reactions.

```julia
using OrdinaryDiffEq, DiffEqDevTools, Sundials, Plots,
LSODA, LinearAlgebra, BenchmarkTools,
SciMLBenchmarks

gr()
```

## Load the model

```julia
include(joinpath(@__DIR__, "Models", "BioPreDyn", "b4_model.jl"))

prob = ODEProblem{true, SciMLBase.FullSpecialize}(
biopredyn_b4!, B4_U0, B4_TSPAN, B4_NOMINAL_PARAMETERS)
@show length(prob.u0)
```

## Picture of the solution

```julia
sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
plot(sol; legend = false, fmt = :png,
title = "BioPreDyn B4: CHO Cell Metabolism (35 ODEs)")
```

## Generate Reference Solution

```julia
@time sol = solve(prob, Rodas5P(), abstol = 1/10^14, reltol = 1/10^14)
test_sol = TestSolution(sol);
```

## Setups

```julia
default(legendfontsize = 7, framestyle = :box, gridalpha = 0.3, gridlinewidth = 2.5)
```

## Work-Precision Diagrams (High Tolerances)

```julia
abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4)

setups = [
Dict(:alg => Rosenbrock23()),
Dict(:alg => TRBDF2()),
Dict(:alg => FBDF()),
Dict(:alg => QNDF()),
Dict(:alg => Rodas4()),
Dict(:alg => Rodas5()),
Dict(:alg => RadauIIA5()),
Dict(:alg => lsoda()),
Dict(:alg => CVODE_BDF()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["Rosenbrock23" "TRBDF2" "FBDF" "QNDF" "Rodas4" "Rodas5" "RadauIIA5" "lsoda" "CVODE_BDF"]
plot(wp; label = names, title = "BioPreDyn B4: High Tolerances")
```

## Work-Precision Diagrams (Low Tolerances)

```julia
abstols = 1.0 ./ 10.0 .^ (7:12)
reltols = 1.0 ./ 10.0 .^ (4:9)

setups = [
Dict(:alg => FBDF()),
Dict(:alg => QNDF()),
Dict(:alg => Rodas4()),
Dict(:alg => Rodas5P()),
Dict(:alg => KenCarp4()),
Dict(:alg => KenCarp47()),
Dict(:alg => Kvaerno5()),
Dict(:alg => CVODE_BDF()),
Dict(:alg => lsoda()),
Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["FBDF" "QNDF" "Rodas4" "Rodas5P" "KenCarp4" "KenCarp47" "Kvaerno5" "CVODE_BDF" "lsoda" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: Low Tolerances")
```

## Work-Precision Diagrams (SDIRK Methods)

```julia
abstols = 1.0 ./ 10.0 .^ (5:10)
reltols = 1.0 ./ 10.0 .^ (2:7)

setups = [
Dict(:alg => KenCarp3()),
Dict(:alg => KenCarp4()),
Dict(:alg => KenCarp5()),
Dict(:alg => Kvaerno3()),
Dict(:alg => Kvaerno4()),
Dict(:alg => Kvaerno5()),
Dict(:alg => TRBDF2()),
Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["KenCarp3" "KenCarp4" "KenCarp5" "Kvaerno3" "Kvaerno4" "Kvaerno5" "TRBDF2" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: SDIRK Methods")
```

## Summary of Results

```julia
abstols = 1.0 ./ 10.0 .^ (5:12)
reltols = 1.0 ./ 10.0 .^ (2:9)

setups = [
Dict(:alg => CVODE_BDF()),
Dict(:alg => QNDF()),
Dict(:alg => FBDF()),
Dict(:alg => Rodas4()),
Dict(:alg => Rodas5P()),
Dict(:alg => KenCarp4()),
Dict(:alg => lsoda()),
Dict(:alg => RadauIIA5()),
]

wp = WorkPrecisionSet(prob, abstols, reltols, setups; error_estimate = :l2,
saveat = B4_TSPAN[2]/1000.0, appxsol = test_sol, maxiters = Int(1e5), numruns = 10)

names = ["CVODE_BDF" "QNDF" "FBDF" "Rodas4" "Rodas5P" "KenCarp4" "lsoda" "RadauIIA5"]
plot(wp; label = names, title = "BioPreDyn B4: Summary",
left_margin = 10Plots.mm, right_margin = 10Plots.mm,
legendfontsize = 10, tickfontsize = 10, guidefontsize = 10,
legend = :topright, lw = 2, markersize = 6, size = (800, 600))
```

```julia, echo = false
using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
```
Loading
Loading