Skip to content
Snippets Groups Projects
Commit 246e1a5c authored by Marco Matthies's avatar Marco Matthies
Browse files

Add mock testing setup for direct translation of ALMaSS crop model

parent 1336bd49
No related branches found
No related tags found
No related merge requests found
abstract type SimulationModel end
struct MockModel <: SimulationModel
end
mintemp(model::MockModel) = 10.0
maxtemp(model::MockModel) = 15.0
function get_plant_growth_data(model)
dds = [ # [5, n_inflections]
[0.0, 10.0, 20.0], # janfirst
[0.0, 10.0, 20.0], # sow
[0.0, 10.0, 20.0], # marchfirst
[0.0, 10.0, 20.0], # harvest1
[0.0, 10.0, 20.0], # harvest2
]
slopes = [ # [5, 3, n_inflections]
[ # janfirst
[1.0, 2.0, 3.0], # LAgreen
[1.0, 2.0, 3.0], # LAtotal
[1.0, 2.0, 3.0], # height
],
[ # sow
[1.0, 2.0, 3.0], # LAgreen
[1.0, 2.0, 3.0], # LAtotal
[1.0, 2.0, 3.0], # height
],
[ # marchfirst
[1.0, 2.0, 3.0], # LAgreen
[1.0, 2.0, 3.0], # LAtotal
[1.0, 2.0, 3.0], # height
],
[ # harvest1
[1.0, 2.0, 3.0], # LAgreen
[1.0, 2.0, 3.0], # LAtotal
[1.0, 2.0, 3.0], # height
],
[ # harvest2
[1.0, 2.0, 3.0], # LAgreen
[1.0, 2.0, 3.0], # LAtotal
[1.0, 2.0, 3.0], # height
],
]
start = [ # [5, 3]
[ # sow
0.0, # LAgreen
0.0, # LAtotal
0.0, # height
],
[ # janfirst
0.0, # LAgreen
0.0, # LAtotal
0.0, # height
],
[ # marchfirst
0.0, # LAgreen
0.0, # LAtotal
0.0, # height
],
[ # harvest1
0.0, # LAgreen
0.0, # LAtotal
0.0, # height
],
[ # harvest2
0.0, # LAgreen
0.0, # LAtotal
0.0, # height
],
]
start_valid = [
true, # janfirst
true, # sow
true, # marchfirst
true, # harvest1
true, # harvest2
]
cg = CropGrowth(; dds, slopes, start, start_valid)
pgd = PlantGrowthData([cg], [1])
return pgd
end
almass_SupplyGlobalRadiation(model::MockModel) = 1.0
almass_SupplyTemp(model::MockModel) = 12.0
almass_SolarConversion(model::MockModel, is_c4_plant, index) = 1.0
almass_SupplyFarmIntensity(model::MockModel, poly::Int) = 1.0
function diff(a::T, b::T) where {T}
if ! isstructtype(T)
error("Only works for struct types, type is $T")
end
println("Diff of $T structs:")
for field in fieldnames(T)
fa = getfield(a, field)
fb = getfield(b, field)
if fa != fb
println(" - $field: ", fa)
println(" + $field: ", fb)
end
end
end
const tov_Maize = 10
const tov_OMaizeSilage = 11
const tov_MaizeSilage = 12
const tov_MaizeStrigling = 13
include("new-almass.jl")
ve = VegElement()
model = MockModel()
ve_old = deepcopy(ve)
println(ve)
println()
almass_DoDevelopment(ve, model)
println(ve)
println()
diff(ve_old, ve)
......@@ -30,32 +30,52 @@
# `class VegElement : public LE` from `Landscape/Elements.h`, line 601
struct VegElement
@kwdef mutable struct VegElement
# Note: this class has a lot of fields, only the fields actually
# used are here
# new_weed_growth::Float64
force_growth::Bool
yddegs::Float64
ddegs::Float64
vegddegs::Float64
curve_num::Int
veg_phase::Int # TODO: enum type?
growth_scaler::Float64
LAgreen::Float64
LAtotal::Float64
veg_height::Float64
force_LAgreen::Float64
force_LAtotal::Float64
force_veg_height::Float64
newgrowth::Float64
veg_cover::Float64
force_growth::Bool = false
yddegs::Float64 = 0
ddegs::Float64 = 0
vegddegs::Float64 = 0
curve_num::Int = 1
veg_phase::Int = 1 # TODO: enum type?
growth_scaler::Float64 = 1
LAgreen::Float64 = 0
LAtotal::Float64 = 0
veg_height::Float64 = 0
force_LAgreen::Float64 = 0
force_LAtotal::Float64 = 0
force_veg_height::Float64 = 0
newgrowth::Float64 = 0
veg_cover::Float64 = 0
vege_type::Int = 1
oldLAtotal::Float64 = 0
biomass_scale::Vector{Float64} = [1.0]
owner_index::Int = 1
poly::Int = 1
veg_biomass::Float64 = 0
area::Float64 = 1
total_biomass::Float64 = 0
total_biomass_old::Float64 = 0
veg_density::Float64 = 0
green_biomass::Float64 = 0
dead_biomass::Float64 = 0
end
function Base.show(io::IO, ve::VegElement)
println(io, "VegElement:")
for field in fieldnames(typeof(ve))
println(io, " $field: ", getfield(ve, field))
end
end
# `typedef enum { ... } Growth_Phases` from `Landscape/Plants.h`, line 51
@enum Growth_Phases begin
janfirst = 0
janfirst = 1
sow
marchfirst
harvest1
......@@ -63,11 +83,11 @@ end
end
# `class CropGrowth` from `Landscape/Plants.h`, line 67
struct CropGrowth
@kwdef struct CropGrowth
# Note:
# - 5 growth phases (enum Growth_Phases)
# - 3 curves (leaf_area_green = 0, leaf_area_total, height)
lownut::Bool # = false
lownut::Bool = false
dds::Vector{Vector{Float64}} # ddegs at inflection, double m_dds[5][MaxNoInflections] = { 0 }
slopes::Vector{Vector{Vector{Float64}}} # double m_slopes[5][3][MaxNoInflections] = { 0 }
start::Vector{Vector{Float64}} # double m_start[5][3] = { 0 }
......@@ -77,16 +97,16 @@ end
# `class PlantGrowthData` from `Landscape/Plants.h`, line 79
struct PlantGrowthData
growth::Vector{CropGrowth}
final_ddeg::Vector{Vector{Vector{Float64}}}
#final_ddeg::Vector{Vector{Vector{Float64}}}
numbers::Vector{Int}
num_crops::Int
# num_crops::Int
# Note: some more fields ommitted related to weeds and bugs
end
# `void VegElement::DoDevelopment(void)` from `Landscape/Elements.cpp`, line 2254
function almass_DoDevelopment(ve::VegElement, model::SimulationModel)
# ve.new_weed_growth = 0.0
if ! ve.force.growth
if ! ve.force_growth
# first do the day degree calculations
ve.yddegs = ve.ddegs
pos_temp_today = almass_GetDDDegs(model)
......@@ -97,9 +117,10 @@ function almass_DoDevelopment(ve::VegElement, model::SimulationModel)
# sum up the phase day degrees
ve.ddegs = pos_temp_today + ve.yddegs
dLAG = almass_GetLAgreenDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
dLAT = almass_GetLAtotalDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
dHgt = almass_GetHeightDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
pgd = get_plant_growth_data(model)
dLAG = almass_GetLAgreenDiffScaled(pgd, ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
dLAT = almass_GetLAtotalDiffScaled(pgd, ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
dHgt = almass_GetHeightDiffScaled(pgd, ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
ve.LAgreen += dLAG
if ve.LAgreen < 0.0
......@@ -142,7 +163,8 @@ end
# `double Weather::GetDDDegs(long a_date)` from `Landscape/Weather.cpp`, line 205
function almass_GetDDDegs(model::SimulationModel)
basetemp = 0.0 # TODO: doesn't seem to use a min growth temp ???
return bounds((maxtemp(model) + mintemp(model)) / 2 - basetemp)
#return bounds((maxtemp(model) + mintemp(model)) / 2 - basetemp)
return max((maxtemp(model) + mintemp(model)) / 2 - basetemp, 0.0)
end
# `double PlantGrowthData::GetLAgreenDiffScaled(...)` from `Landscape/Plants.h`, line 136
......@@ -166,22 +188,22 @@ end
# `double PlantGrowthData::GetLAgreenDiff(...)` from `Landscape/Plants.h`, line 106
function almass_GetLAgreenDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
curve_num::Int, veg_phase::Int)
# TODO: hardcoded 0
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 0)
# TODO: hardcoded 1, in original code it is 0
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 1)
end
# `double PlantGrowthData::GetLAtotalDiff(...)` from `Landscape/Plants.h`, line 116
function almass_GetLAtotalDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
curve_num::Int, veg_phase::Int)
# TODO: hardcoded 1
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 1)
# TODO: hardcoded 2, in original code it is 1
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 2)
end
# `double PlantGrowthData::GetHeightDiff(...)` from `Landscape/Plants.h`, line 125
function almass_GetHeightDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
curve_num::Int, veg_phase::Int)
# TODO: hardcoded 2
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 2)
# TODO: hardcoded 3, in original code it is 2
return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 3)
end
# `double PlantGrowthData::FindDiff(...)` from `Landscape/Plants.cpp`, line 45
......@@ -301,7 +323,7 @@ function almass_RecalculateBugsNStuff(ve::VegElement, model::SimulationModel)
index = Int(floor(0.5 + almass_SupplyTemp(model)) + 30 + 1) # // There are 30 negative temps
# TODO: model needs this field???
# Note: +1 because Julia arrays are 1-based
radconv = model.SolarConversion[is_c4_plant + 1][index + 1]
radconv = almass_SolarConversion(model, is_c4_plant, index)
if ve.LAtotal >= ve.oldLAtotal
# we are in positive growth so grow depending on our equation
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment