Corrected domain issues

Added a `ProblemDomain` type and a domain check function the generalize the GaussianPlumeSolution. Previously it was ignoring the mixing layer when deciding on the domain of the problem.
Added a limiter to plume rise, forcing plume rise below the mixing layer.

Author: aefarrell

src/models/gaussian_mixing_layer.jl

@@ -107,6 +107,9 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet()
         error("Unknown mixing layer method: $method")
     end
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ)
+
     return GaussianPlumeSolution(
     scenario,                               # scenario::Scenario
     :simplemixinglayer,                     # model::Symbol
@@ -119,8 +122,8 @@ function plume(scenario::Scenario, ::Type{GaussianMixingLayer}, eqs=DefaultSet()
     ml,
     NoPlumeRise(),                          # plume rise model
     _stability(scenario),                   # stability class
-    eqs                                     # equation set 
-    )
+    eqs,                                    # equation set 
+    dom)                                    # problem domain
 end
 
 @doc doc"""
@@ -200,6 +203,10 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     if plumerise == true
         Tᵣ = _release_temperature(scenario)
         plume = plume_rise(Dⱼ,uⱼ,Tᵣ,u,Tₐ,Γ,stab)
+        if plume.final_rise > hₘ-hᵣ
+            @warn "Plume rise $(plume.final_rise) exceeds mixing height $hₘ, plume rise is limited to the mixing height."
+            plume = _new_plume_rise(plume, hₘ)
+        end
     else
         plume = NoPlumeRise()
     end
@@ -217,6 +224,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
         error("Unknown mixing layer method: $method")
     end
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, hₘ)
+
     return GaussianPlumeSolution(
     scenario, #scenario::Scenario
     :simplemixinglayer, #model::Symbol
@@ -229,6 +239,6 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     ml,
     plume, #plume rise model
     stab,  #stability class
-    eqs    #equation set 
-    )
+    eqs,   #equation set 
+    dom)   #problem domain
 end

src/models/gaussian_plume.jl

@@ -6,7 +6,7 @@ abstract type GaussianVerticalTerm end
 struct GaussianPlume <: PlumeModel end
 
 # Solution to the gaussian plume
-struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet} <: Plume
+struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalTerm,P<:PlumeRise,S<:StabilityClass,E<:EquationSet,D<:ProblemDomain} <: Plume
     scenario::Scenario
     model::Symbol
     rate::F
@@ -19,8 +19,9 @@ struct GaussianPlumeSolution{F<:Number,C<:GaussianCrossTerm,V<:GaussianVerticalT
     plumerise::P
     stability::Type{S}
     equationset::E
+    domain::D
 end
-GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es)
+GaussianPlumeSolution(s,m,Q,c,ρ,u,h_eff,cross,vert,pr,stab,es,dom) = GaussianPlumeSolution(s,m,promote(Q,c,ρ,u,h_eff)...,cross,vert,pr,stab,es,dom)
 
 struct SimpleCrossTerm <:  GaussianCrossTerm end
 cross_term(y, σy, ::SimpleCrossTerm) = exp(-0.5*(y/σy)^2)/(√(2π)*σy)
@@ -67,20 +68,23 @@ function plume(scenario::Scenario, ::Type{GaussianPlume}, eqs=DefaultSet(); h_mi
     c_max = m/Qᵒ
     c_max = c_max/ρₐ
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf)
+
     return GaussianPlumeSolution(
     scenario,                               # scenario::Scenario
     :gaussian,                              # model::Symbol
     m,                                      # mass emission rate
     c_max,                                  # max_concentration
-    ρₐ,                                      # mass concentration to vol concentration
-    windspeed(scenario,max(hᵣ,h_min),eqs), # windspeed
+    ρₐ,                                     # mass concentration to vol concentration
+    windspeed(scenario,max(hᵣ,h_min),eqs),  # windspeed
     hᵣ,                                     # effective_stack_height::Number
     SimpleCrossTerm(),
     SimpleVerticalTerm(),
     NoPlumeRise(),                          # plume rise model
     _stability(scenario),                   # stability class
-    eqs                                     # equation set 
-    )
+    eqs,                                    # equation set 
+    dom)                                    # problem domain
 end
 
 @doc doc"""
@@ -151,6 +155,9 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
         plume = NoPlumeRise()
     end
 
+    # domain
+    dom = ProblemDomain(0.0, Inf, -Inf, Inf, 0.0, Inf)
+
     return GaussianPlumeSolution(
     scenario, #scenario::Scenario
     :gaussian, #model::Symbol
@@ -163,16 +170,16 @@ function plume(scenario::Scenario{<:AbstractSubstance,<:VerticalJet,<:Atmosphere
     SimpleVerticalTerm(),
     plume, #plume rise model
     stab,  #stability class
-    eqs    #equation set 
-    )
+    eqs,   #equation set 
+    dom)   #problem domain
 end
 
 function (g::GaussianPlumeSolution{F,C,V,NoPlumeRise,S,E})(x, y, z, t=0) where {F,C,V,S,E}
     # domain check
     h = g.effective_stack_height
     if (x==0)&&(y==0)&&(z==h)
         return g.max_concentration
-    elseif (x≤0)||(z<0)
+    elseif _in_domain(x,y,z,g.domain) == false
         return zero(F)
     else
         G = g.rate
@@ -197,7 +204,7 @@ function (g::GaussianPlumeSolution{F,C,V,<:BriggsModel,S,E})(x, y, z, t=0) where
     h = g.effective_stack_height
     if (x==0)&&(y==0)&&(z==h)
         return g.max_concentration
-    elseif (x≤0)||(z<0)
+    elseif _in_domain(x,y,z,g.domain) == false
         return zero(F)
     else
         G = g.rate

src/utils/plume_rise.jl

@@ -31,6 +31,16 @@ Base.isapprox(a::MomentumPlume, b::MomentumPlume) = all([
     if typeof(getproperty(a,k))<:Number ])
 
 
+# some helper functions to limit the plume rise, e.g. when a mixing layer prevents the plume from rising too high
+"""
+    _new_plume_rise(plume::BriggsModel, max_rise::Number)
+Creates a new plume rise model with the same parameters as `plume`, but limits the final
+rise to `max_rise`. This is useful when the plume rise exceeds the mixing height or
+other constraints in the model.
+"""
+_new_plume_rise(plume::BuoyantPlume, max_rise::Number) = BuoyantPlume(plume.Fb, plume.xf, plume.u, max_rise)
+_new_plume_rise(plume::MomentumPlume, max_rise::Number) = MomentumPlume(plume.Fm, plume.xf, plume.β, plume.u, plume.s, max_rise, plume.stab)
+
 """
     plume_rise(Dⱼ,uⱼ,Tᵣ,a::Atmosphere)
 Implements the Briggs plume rise equations for buoyancy and momentum driven

src/utils/utils.jl

@@ -27,4 +27,20 @@ include("equation_sets/isc3.jl")
 
 # Default correlations
 DefaultSet = BasicEquationSet{DefaultWind,Nothing,Defaultσy,Defaultσz}
-DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz}
+DefaultPuffSet = BasicEquationSet{DefaultWind,CCPSPuffσx,CCPSPuffσy,CCPSPuffσz}
+
+# a lazy check for if x,y,z are in the domain
+struct ProblemDomain{F<:Number}
+    xmin::F
+    xmax::F
+    ymin::F
+    ymax::F
+    zmin::F
+    zmax::F
+end
+
+function _in_domain(x::Number, y::Number, z::Number, domain::ProblemDomain)
+    return x ≥ domain.xmin && x ≤ domain.xmax &&
+           y ≥ domain.ymin && y ≤ domain.ymax &&
+           z ≥ domain.zmin && z ≤ domain.zmax
+end

test/models/gaussian_plume_tests.jl

@@ -30,7 +30,7 @@
     pl = plume(scn)
 
     # test default behaviour and type inheritance
-    @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet()) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}}
+    @test GasDispersion.GaussianPlumeSolution(scn,:test,1.0,2,3,4,5,GasDispersion.SimpleCrossTerm(),GasDispersion.SimpleVerticalTerm(),GasDispersion.NoPlumeRise(),ClassA,DefaultSet(),GasDispersion.ProblemDomain(0.0,Inf,-Inf,Inf,0.0,Inf)) isa GasDispersion.GaussianPlumeSolution{Float64, GasDispersion.SimpleCrossTerm, GasDispersion.SimpleVerticalTerm, GasDispersion.NoPlumeRise, ClassA, GasDispersion.BasicEquationSet{GasDispersion.DefaultWind, Nothing, GasDispersion.Defaultσy, GasDispersion.Defaultσz}, GasDispersion.ProblemDomain{Float64}}
     @test pl isa GasDispersion.GaussianPlumeSolution
     @test pl isa Plume
     @test pl(-1,0,0) == 0.0