Fix regression in DiffBioEq driver

.gitignore

@@ -0,0 +1 @@
+models/*

KPP.jl

@@ -1,5 +1,5 @@
 # Configure here:
-const mechanism = "saprc99"
+const mechanism = "small_strato"
 const mechanism_safe = mechanism
 
 # Choose the driver to use: DiffBioEq or Native (ODE constructed from scratch)
@@ -8,7 +8,7 @@ const driver = "Native"
 # Enable adaptive optimizer? This is an implementation of the adaptive reduce
 # algorithm described by Santillana et al., 2010.
 # Only compatible with the "Native" driver
-const incl_optimize_adapt = true
+const incl_optimize_adapt = false
 
 # NO USER CONFIGURABLE PARAMETERS BELOW
 include("parser.jl")

generator.jl

@@ -49,6 +49,7 @@ function generatemodel()
  # TIMESTEP_PARAM (in jlkpp_Timestep)
  # TIMESTEP_AFTER_PARAM***
  # TIMESTEP_AFTER_SOLVE
+ # TIMESTEP_SPC_LOOP_COND_LEFT, TIMESTEP_SPC_LOOP_COND_RIGHT (usually end)
  # KPP LAYER: Headers/registry.jl
  # REGISTRY_SPCLIST
  # REGISTRY_DEFAULT_IC (this is passed as jlkpp_ic to INITIAL_CONDITIONS)
@@ -60,6 +61,8 @@ function generatemodel()
  initialize_kpp_spc = ""
  initialize_kpp_eqn = ""
  timestep_param = ""
+ timestep_spc_loop_cond_left = ""
+ timestep_spc_loop_cond_right = ""
 
  registry_spclist = ""
  registry_default_ic = ""
@@ -91,13 +94,20 @@ addodes!(jlkpp_mechanism)
  # Adaptive reduction (Santillana et al., 2010)
  # Add extra code to zero out changes for the "fixed" slow species
  if incl_optimize_adapt == true
- initialize_kpp_eqn *= """
- # Adaptive reduction: slow species set rate of change to 0; concentrations
- # will be manually updated later.
- for spc in 1:length(u)
- jlkpp_adapt_slow[spc] && (du[spc] = 0)
- end
-"""
+# initialize_kpp_eqn *= """
+# # Adaptive reduction: slow species set rate of change to 0; concentrations
+# # will be manually updated later.
+# for spc in 1:length(u)
+# jlkpp_adapt_slow[spc] && (du[spc] = 0)
+# end
+# """
+
+ # Avoid the compute in du by using ternary operators -
+ # this is experimental for the sake of benchmarking
+ initialize_kpp_eqn = replace(initialize_kpp_eqn, r"du\[(?<id>\d+)\] = (?<expr>.*)" => s"du[\g<id>] = jlkpp_adapt_slow[\g<id>] ? 0 : \g<expr>")
+
+ timestep_spc_loop_cond_left = "if jlkpp_adapt_slow[spc] == false"
+ timestep_spc_loop_cond_right = "end"
  end
 
  # Build the final KPP block code
@@ -122,7 +132,7 @@ end
  # For DiffBioEq its just an include.
  # For Native usually there are optimization options, etc. you can use
  setup_extra_declares = ""
- setup_extra_declares_kpp = ""
+ setup_extra_declares_kpp = "const jlkpp_nspecies = " * string(modelparams["species"]) * "\n"
  generate_oprob = ""
  if driver == "DiffBioEq"
  setup_extra_declares *= "using DiffEqBiological\n"
@@ -131,7 +141,7 @@ end
 
  # generate_oprob *= "jlkpp_oprob = jlkpp_Compile(jlkpp_mechanism, @view(chem_species[:,1,1,1]))"
  # @view is incompatble with SUNDIALS solver
- generate_oprob *= "jlkpp_oprob = jlkpp_Compile(jlkpp_mechanism_t!, chem_species[1,1,1])"
+ generate_oprob *= "jlkpp_oprob = jlkpp_Compile(jlkpp_mechanism, chem_species[1,1,1])"
  elseif driver == "Native"
  setup_extra_declares *= "# KPP.jl Native Driver: Use precomputed Jacobian?\n"
  setup_extra_declares *= "const jlkpp_prejac = true\n"
@@ -143,7 +153,7 @@ end
  setup_extra_declares *= "using Sundials\n"
 
  setup_extra_declares *= "# KPP.jl: Time stepping method. 'remake' remakes ODE problem, 'interfaced' uses the DiffEq Integrator interface.\n"
- setup_extra_declares *= "const jlkpp_stepmethod = 'remake'\n"
+ setup_extra_declares *= "const jlkpp_stepmethod = \"remake\"\n"
 
  # Model parameters write into main.jl as necessary
  setup_timesteps *= string("# These timesteps are set here by KPP.jl for convenience. If you are building a model, you should be able to read these from a configuration file.\n")
@@ -157,7 +167,7 @@ end
  kpp_compile = """
  function jlkpp_Compile(rs::AbstractReactionNetwork, u_scratch::AbstractArray{Float64, 1})::ODEProblem
  # println("KPP.jl: Running ODE solver for type specialization")
- ODEProblem(rs, u_scratch, (0.0, 600.0), t -> (0, 300.0, 2.4476e13))
+ ODEProblem(rs, u_scratch, (0.0, 600.0), (0, 300.0, 2.4476e13))
  # solve(oprob, alg_hints=[:stiff], dense=false, save_on=false, calck=false)
  # println("KPP.jl: Finished")
  # oprob # Return the ODEProblem
@@ -210,7 +220,11 @@ end
  end
 
  timestep_param *= " # SUN, TEMP, CFACTOR\n"
- timestep_param *= string(" p = t -> (jlkpp_SUN(t), T_chm, ", modelparams["CFACTOR"], timestep_param_extra_code, ")\n")
+ if driver == "Native"
+ timestep_param *= string(" p = t -> (jlkpp_SUN(t), T_chm, ", modelparams["CFACTOR"], timestep_param_extra_code, ")\n")
+ else
+ timestep_param *= string(" p = (jlkpp_SUN(t), T_chm, ", modelparams["CFACTOR"], timestep_param_extra_code, ")\n")
+ end
 
  # ----------------------- END PREP ------------------------- #
 

optimizers/santillana2010.jl

@@ -18,7 +18,7 @@ function jlkpp_adapt_diagn(u::Array{Float64, 1}, t::Float64)
 
  # Outs
  _optim_slow_count = 0
- for spc in 1:length(u)
+ for spc in 1:jlkpp_nspecies
  jlkpp_adapt_slow[spc] = (jlkpp_scratch_p[spc] < δ && jlkpp_scratch_l[spc] < δ)
  jlkpp_adapt_slow[spc] && (_optim_slow_count += 1)
  end
@@ -35,9 +35,21 @@ function jlkpp_adapt_slow!(u::Array{Float64, 1}, dt::Float64)
  # so we just assume ki* = Li/Ci, which is fair enough as most equations are like this.
 
  # We use the global scratchpad
- for spc in 1:length(u)
- if jlkpp_adapt_slow[spc] == true # is slow spc
- u[spc] = (jlkpp_scratch_p[spc]/jlkpp_scratch_l[spc] + (1.0-jlkpp_scratch_p[spc])*exp(-jlkpp_scratch_l[spc]*dt/u[spc]))*u[spc]
+ for spc in 1:jlkpp_nspecies
+ if jlkpp_adapt_slow[spc] # is slow spc
+ if jlkpp_scratch_l[spc] != 0
+ # print("* jlkpp_adapt_slow! updating slow spc ", spc, " from ", u[spc])
+ # print(" p/l=", jlkpp_scratch_p[spc]/jlkpp_scratch_l[spc], " exp=", exp(-jlkpp_scratch_l[spc]*dt/u[spc]))
+ u[spc] = (jlkpp_scratch_p[spc]/jlkpp_scratch_l[spc] + (1.0-jlkpp_scratch_p[spc]/jlkpp_scratch_l[spc])*exp(-jlkpp_scratch_l[spc]*dt/u[spc]))*u[spc]
+ # println("-> to ", u[spc])
+ else
+ # print("* jlkpp_adapt_slow! updating slow spc ", spc, " from ", u[spc])
+ u[spc] = jlkpp_scratch_p[spc]*dt + u[spc]
+ # println("-> using linear method to ", u[spc])
+ end
+
+ # Never allow negative species conc
+ u[spc] < 0 && (u[spc] = 0)
  end
  end
 end

parser.jl

@@ -93,13 +93,11 @@ function parsekcoeff(str::AbstractString)
  # Params: TEMP
  rgx_ARR = r"ARR\(\s*(?<A0>[\d\.e\-\+]+)\s*,\s*(?<B0>[\d\.e\-\+]+)\s*,\s*(?<C0>[\d\.e\-\+]+)\s*\)"
  tpl_ARR = s"jlkpp_ARR(TEMP, \g<A0>, \g<B0>, \g<C0>)"
- # tpl_ARR = s"((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>))"
 
  # Simplified Arrhenius, with two arguments
  # Params: TEMP
  rgx_ARR2 = r"ARR2\(\s*(?<A0>[\d\.e\-\+]+)\s*,\s*(?<B0>[\d\.e\-\+]+)\)"
  tpl_ARR2 = s"jlkpp_ARR2(TEMP, \g<A0>, \g<B0>)"
- # tpl_ARR2 = s"((\g<A0>) * exp((\g<B0>)/TEMP))"
 
  # EP2
  # K0 = (:A0) * exp(-(:C0)/TEMP)
@@ -108,22 +106,28 @@ function parsekcoeff(str::AbstractString)
  # Params: TEMP, CFACTOR
  rgx_EP2 = r"EP2\(\s*(?<A0>[\d\.e\-\+]+)\s*,\s*(?<C0>[\d\.e\-\+]+)\s*,\s*(?<A2>[\d\.e\-\+]+)\s*,\s*(?<C2>[\d\.e\-\+]+)\s*,\s*(?<A3>[\d\.e\-\+]+)\s*,\s*(?<C3>[\d\.e\-\+]+)\s*\)"
  tpl_EP2 = s"jlkpp_EP2(TEMP, CFACTOR, \g<A0>, \g<C0>, \g<A2>, \g<C2>, \g<A3>, \g<C3>)"
- # tpl_EP2 = s"((\g<A0>) * exp(-(\g<C0>)/TEMP) + ((\g<A3>) * exp(-(\g<C3>)/TEMP) * 1.0e6 * CFACTOR)/(1.0 + ((\g<A3>) * exp(-(\g<C3>)/TEMP) * 1.0e6 * CFACTOR)/((\g<A2>) * exp(-(\g<C2>)/TEMP))))"
 
  # EP3
  # K1 = (:A1) * exp(-(:C1)/TEMP)
  # K2 = (:A2) * exp(-(:C2)/TEMP)
  rgx_EP3 = r"EP3\(\s*(?<A1>[\d\.e\-\+]+)\s*,\s*(?<C1>[\d\.e\-\+]+)\s*,\s*(?<A2>[\d\.e\-\+]+)\s*,\s*(?<C2>[\d\.e\-\+]+)\s*\)"
  tpl_EP3 = s"jlkpp_EP3(TEMP, CFACTOR, \g<A1>, \g<C1>, \g<A2>, \g<C2>)"
- # tpl_EP3 = s"((\g<A1>) * exp(-(\g<C1>)/TEMP) + ((\g<A2>) * exp(-(\g<C2>)/TEMP)) * 1.0e6 * CFACTOR)"
 
  # FALL
  # K0 = (:A0) * exp(-(:B0)/TEMP) * (TEMP/300.0)^(:C0) * CFACTOR * 1.0e6
  # K1 = (:A1) * exp(-(:B1)/TEMP) * (TEMP/300.0)^(:C1)
  # K1' = ((:A0) * exp(-(:B0)/TEMP) * (TEMP/300.0)^(:C0) * CFACTOR * 1.0e6)/((:A1) * exp(-(:B1)/TEMP) * (TEMP/300.0)^(:C1))
  rgx_FALL = r"FALL\(\s*(?<A0>[\d\.e\-\+]+)\s*,\s*(?<B0>[\d\.e\-\+]+)\s*,\s*(?<C0>[\d\.e\-\+]+)\s*,\s*(?<A1>[\d\.e\-\+]+)\s*,\s*(?<B1>[\d\.e\-\+]+)\s*,\s*(?<C1>[\d\.e\-\+]+)\s*,\s*(?<CF>[\d\.e\-\+]+)\s*\)"
  tpl_FALL = s"jlkpp_FALL(TEMP, CFACTOR, \g<A0>, \g<B0>, \g<C0>, \g<A1>, \g<B1>, \g<C1>, \g<CF>)"
- # tpl_FALL = s"((((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/(1.0 + ((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/((\g<A1>) * exp(-(\g<B1>)/TEMP) * (TEMP/300.0)^(\g<C1>)))) * (\g<CF>)^(1.0/(1.0 + log10((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/((\g<A1>) * exp(-(\g<B1>)/TEMP) * (TEMP/300.0)^(\g<C1>)))^2))"
+
+ if driver == "DiffBioEq"
+ tpl_ARR = s"((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>))"
+ tpl_ARR2 = s"((\g<A0>) * exp((\g<B0>)/TEMP))"
+ tpl_EP2 = s"((\g<A0>) * exp(-(\g<C0>)/TEMP) + ((\g<A3>) * exp(-(\g<C3>)/TEMP) * 1.0e6 * CFACTOR)/(1.0 + ((\g<A3>) * exp(-(\g<C3>)/TEMP) * 1.0e6 * CFACTOR)/((\g<A2>) * exp(-(\g<C2>)/TEMP))))"
+ tpl_EP3 = s"((\g<A1>) * exp(-(\g<C1>)/TEMP) + ((\g<A2>) * exp(-(\g<C2>)/TEMP)) * 1.0e6 * CFACTOR)"
+ tpl_FALL = s"((((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/(1.0 + ((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/((\g<A1>) * exp(-(\g<B1>)/TEMP) * (TEMP/300.0)^(\g<C1>)))) * (\g<CF>)^(1.0/(1.0 + log10((\g<A0>) * exp(-(\g<B0>)/TEMP) * (TEMP/300.0)^(\g<C0>) * CFACTOR * 1.0e6)/((\g<A1>) * exp(-(\g<B1>)/TEMP) * (TEMP/300.0)^(\g<C1>)))^2))"
+ # FIXME HACK: DiffBioEq has problems interpolating functions so we have to use the old way...
+ end
 
  # replace away, then tokenize
  str = replace(str, rgx_ARR => tpl_ARR)
@@ -133,7 +137,7 @@ function parsekcoeff(str::AbstractString)
  str = replace(str, rgx_FALL => tpl_FALL)
 
  if driver == "DiffBioEq"
- string(":(", str, ")") # As meta-expression
+ string(":(", str, ")")
  elseif driver == "Native"
  string("(", str, ")")
  end
@@ -278,18 +282,21 @@ function readfile_eqn(spclist::Array{String, 1}, spcfixID::Array{Int, 1}, modelp
  netrates[ridx] *= " + $rcoef * $rate_expr"
  prod[ridx] *= " + $rcoef * $rate_expr"
  end
- end
- end
-
- codeeqn *= join(netrates, "\r\n") * "\r\n"
- prodloss_diag_kpp_eqn *= join(prod, "\r\n") * "\r\n\r\n"
- prodloss_diag_kpp_eqn *= join(loss, "\r\n") * "\r\n"
+ end # driver == "Native"
+ end # for each match
 
- # For cleaner code we can remove the initial "= 0 +" or "= 0 -" if they exist
- codeeqn = replace(codeeqn, "= 0 +" => "=")
- codeeqn = replace(codeeqn, "= 0 -" => "= -")
- prodloss_diag_kpp_eqn = replace(prodloss_diag_kpp_eqn, "= 0 +" => "=")
- prodloss_diag_kpp_eqn = replace(prodloss_diag_kpp_eqn, "= 0 -" => "= -")
+ if driver == "Native"
+ # Build the ODE and Prod/Loss Diagn for the Native driver...
+ codeeqn *= join(netrates, "\r\n") * "\r\n"
+ prodloss_diag_kpp_eqn *= join(prod, "\r\n") * "\r\n\r\n"
+ prodloss_diag_kpp_eqn *= join(loss, "\r\n") * "\r\n"
+
+ # For cleaner code we can remove the initial "= 0 +" or "= 0 -" if they exist
+ codeeqn = replace(codeeqn, "= 0 +" => "=")
+ codeeqn = replace(codeeqn, "= 0 -" => "= -")
+ prodloss_diag_kpp_eqn = replace(prodloss_diag_kpp_eqn, "= 0 +" => "=")
+ prodloss_diag_kpp_eqn = replace(prodloss_diag_kpp_eqn, "= 0 -" => "= -")
+ end
 
  modelparams["reactions"] = rxn_count
 

template/Core/kpp.jl

@@ -9,7 +9,7 @@ from the original KPP specification file and imported into the registry
 into the given array, by grid box
 """
 function jlkpp_Initialize_Defaults(u::Array{Float64, 1})
- @inbounds for spc in 1:length(jlkpp_bg_ic)
+ @inbounds for spc in 1:jlkpp_nspecies
  u[spc] = jlkpp_bg_ic[spc]
  end
 end
@@ -80,6 +80,13 @@ function jlkpp_Timestep!(
 
  # $TIMESTEP_AFTER_SOLVE$ #
 
+ # Update species concentrations as this is the mutating form.
+ # Allow hooks
+ for spc in 1:jlkpp_nspecies
+ # $TIMESTEP_SPC_LOOP_COND_LEFT$ #
+ u[spc] = _osol.u[end][spc]
+ # $TIMESTEP_SPC_LOOP_COND_RIGHT$ #
+ end
 end
 
 # $KPP_CODE_EXTRA$ #