Update README.md

README.md

@@ -9,50 +9,52 @@ Author: Brian Groenke (brian.groenke@awi.de)
 
 ### Quick start
 
-```julia
-using CryoGrid
-using CryoGrid.Models
-using Dates
-using Plots
+Single layer heat conduction model with free water freeze curve and air temperature upper boundary condition:
 
-forcings = loadforcings("input/FORCING_JSONfiles/FORCING_ULC_126_72.json", :Tair => u"°C");
-# use air temperature as upper boundary forcing
+```julia
+# load provided forcing data from Samoylov;
+# The forcing file will be automatically downloaded to the input/ folder if not already present.
+forcings = loadforcings(Models.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
+# use air temperature (converted to Kelvin) as upper boundary forcing
 tair = TimeSeriesForcing(ustrip.(u"K", forcings.data.Tair), forcings.timestamps, :Tair);
 # basic 1-layer heat conduction model (defaults to free water freezing scheme)
 model = Models.SoilHeat(TemperatureGradient(tair), SamoylovDefault)
-# define time span
-tspan = (DateTime(2010,1,1),DateTime(2010,12,31))
+# define time span (1 year)
+tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
 # CryoGrid front-end for ODEProblem
 prob = CryoGridProblem(model,tspan)
 # solve discretized system, saving every 6 hours;
-# ROS3P is a third order Rosenbrock solver that should work well without a freeze curve.
-out = @time solve(prob, ROS3P(), abstol=1e-2, saveat=6*3600.0) |> CryoGridOutput;
-zs = [1:10...,20:10:100...]
-cg = Plots.cgrad(:berlin,rev=true)
-plot(out.soil.T[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, dpi=150)
+# Trapezoid on a discretized PDE is analogous to the well known Crank-Nicolson method.
+out = @time solve(prob, Trapezoid(), abstol=1e-3, reltol=1e-4, saveat=6*3600.0, progress=true) |> CryoGridOutput;
+zs = [1.0,5,10,20,30,50,100,500,1000]u"cm"
+cg = Plots.cgrad(:copper,rev=true)
+plot(out.soil.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, dpi=150)
+```
+![Ts_output_freew](res/Ts_H_tair_freeW_2010-2011.png)
 
-# Alternatively, we can use a van Genuchten freeze curve
+Alternatively, we can use a van Genuchten freeze curve:
+
+```julia
 model = Models.SoilHeat(TemperatureGradient(tair), SamoylovDefault, freezecurve=SFCC(VanGenuchten()))
 # Set-up parameters
 p = copy(model.pproto)
 p.soil.α .= 4.0
 p.soil.n .= 2.0
 p.soil.Tₘ .= 273.15
-tspan = (DateTime(2010,1,1),DateTime(2010,12,31))
+tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
 prob = CryoGridProblem(model,tspan,p)
 # stiff solvers don't work well with van Genuchten due to the ill-conditioned Jacobian;
-# Thus, we use forward Euler instead
-out = @time solve(prob, Euler(), dt=2*60.0, saveat=6*3600.0) |> CryoGridOutput;
-zs = [1:10...,20:10:100...]
-cg = Plots.cgrad(:berlin,rev=true)
-plot(out.soil.T[Z(zs)], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, dpi=150)
+# We can just forward Euler instead.
+out = @time solve(prob, Euler(), dt=120.0, saveat=6*3600.0, progress=true) |> CryoGridOutput;
+plot(out.soil.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, dpi=150)
 ```
+![Ts_output_vgfc](res/Ts_H_tair_vg_2010-2011.png)
 
 Note that `SoilHeat` uses energy as the state variable by default. To use temperature as the state variable instead:
 
 ```julia
 # Note that this will work with any freeze curve, here we use Westermann (2011).
-# u"K" is the unit (Kelvin) for temperature, u"J" (Joules) represents energy.
+# :T is the variable name for temperature, :H represents enthalpy/energy.
 # This is used in the specification of the Heat process type.
-model = Models.SoilHeat(u"K", TemperatureGradient(tair), SamoylovDefault, freezecurve=SFCC(Westermann()))
+model = Models.SoilHeat(:T, TemperatureGradient(tair), SamoylovDefault, freezecurve=SFCC(Westermann()))
 ```