Skip to content

Commit cf4e1d8

Browse files
hwpanghwpang
committed
Calculate derivative and jacobian for volumetric flow rate for all domains
add missing - add missing - rename to Hpervolume
1 parent 48128ca commit cf4e1d8

File tree

1 file changed

+114
-0
lines changed

1 file changed

+114
-0
lines changed

‎src/Domain.jl

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1682,6 +1682,8 @@ export calcthermo
16821682
cond = kLAs.*inter.molefractions*inter.P/R/inter.T./kHs*V
16831683

16841684
dydt[d.indexes[1]:d.indexes[2]] .-= (evap .- cond)
1685+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1686+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
16851687
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
16861688
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
16871689
end
@@ -1707,6 +1709,9 @@ end
17071709
cond = kLAs.*ns./kHs
17081710

17091711
dydt[d.indexes[1]:d.indexes[2]] .+= (evap .- cond)
1712+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1713+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1714+
dydt[d.indexes[3]] += inter.Vin(t)
17101715
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
17111716
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
17121717
dydt[d.indexes[3]] -= inter.Vout(t)
@@ -1757,6 +1762,11 @@ end
17571762
dTdt = (P*V/N*flow)/(N*Cvave)
17581763
dydt[d.indexes[3]] -= dTdt
17591764
dydt[d.indexes[4]] -= flow*R*T/V + P/T*dTdt
1765+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1766+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1767+
dTdt = inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)
1768+
dydt[d.indexes[3]] += dTdt
1769+
dydt[d.indexes[4]] += inter.Vin(t)*P/V + P/T*dTdt
17601770
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
17611771
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
17621772
dTdt = (P*inter.Vout(t))/(N*Cvave)
@@ -1798,6 +1808,11 @@ end
17981808

17991809
flow = sum(cond)
18001810
dydt[d.indexes[4]] -= flow*R*T/P
1811+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1812+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1813+
dTdt = inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
1814+
dydt[d.indexes[3]] += dTdt
1815+
dydt[d.indexes[4]] += inter.Vin(t)
18011816
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
18021817
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
18031818
dydt[d.indexes[4]] -= inter.Vout(t)
@@ -1839,6 +1854,9 @@ end
18391854

18401855
flow = sum(cond)
18411856
dydt[d.indexes[3]] -= flow*R*T/P
1857+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1858+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1859+
dydt[d.indexes[3]] += inter.Vin(t)
18421860
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
18431861
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
18441862
dydt[d.indexes[3]] -= inter.Vout(t)
@@ -1890,6 +1908,11 @@ end
18901908
dTdt = (P*V/N*flow)/(N*Cvave)
18911909
dydt[d.indexes[3]] -= dTdt
18921910
dydt[d.indexes[4]] -= flow*R*T/V + dTdt*P/T
1911+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1912+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1913+
dTdt = inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)
1914+
dydt[d.indexes[3]] += dTdt
1915+
dydt[d.indexes[4]] += inter.Vin(t)*P/V + dTdt*P/T
18931916
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
18941917
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
18951918
dTdt = (P*inter.Vout(t))/(N*Cvave)
@@ -1933,6 +1956,11 @@ end
19331956
flow = sum(cond)
19341957
dydt[d.indexes[1]:d.indexes[2]] .-= flow.*ns./N
19351958
dydt[d.indexes[4]] -= flow*R*T/P
1959+
elseif isa(inter,VolumetricFlowRateInlet) && d == inter.domain
1960+
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t)*inter.cs
1961+
dTdt = inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
1962+
dydt[d.indexes[3]] += dTdt
1963+
dydt[d.indexes[4]] += inter.Vin(t)
19361964
elseif isa(inter,VolumetricFlowRateOutlet) && d == inter.domain
19371965
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t)*ns/V
19381966
dydt[d.indexes[4]] -= inter.Vout(t)
@@ -2092,6 +2120,15 @@ end
20922120
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
20932121
end
20942122
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= kLAs./kHs*R*T/P
2123+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2124+
# dn/dt += inter.Vin(t)*inter.cs
2125+
# dV/dt += inter.Vin(t)
2126+
# d/dn(dn/dt) += 0
2127+
# d/dV(dn/dt) += 0
2128+
# d/dn(dV/dt) += 0
2129+
# d/dV(dV/dt) += 0
2130+
nothing
2131+
20952132
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
20962133
# outlet
20972134
# d/dni(dV/dt) = dflow/dni*R*T/P = inter.Vout(t)/V*R*T/P = inter.Vout(t)/N
@@ -2274,6 +2311,24 @@ end
22742311
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
22752312
@inbounds @fastmath jac[domain.indexes[4],i] -= kLAs[i]/kHs[i]*R*T/V + P/T*ddnidTdt
22762313
end
2314+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2315+
# dn/dt += inter.Vin(t)*inter.cs
2316+
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)
2317+
# dP/dt += inter.Vin(t)*P/V + P/T*dTdt
2318+
# d/dni(dT/dt) += inter.Vin(t)*(-Us[i]/N)/(N*Cvave) - inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)^2 * d/dni(N*Cvave)
2319+
# += inter.Vin(t)*(-Us[i]/N)/(N*Cvave) - dTdt * d/dni(N*Cvave)/(N*Cvave)
2320+
# += inter.Vin(t)*(-Us[i]/N)/(N*Cvave) - dTdt * d/dni(Cvave)/Cvave
2321+
# Note: Cvave = sum(ns.*cpdivR)*R/N - R
2322+
# Note: d/dni(Cvave) = cpdivR[i]*R/N
2323+
# d/dni(dP/dt) += P/T * d/dni(dT/dt)
2324+
@fastmath dTdt = (inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V))/(N*Cvave)
2325+
@simd for i in domain.indexes[1]:domain.indexes[2]
2326+
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2327+
@inbounds @fastmath ddnidTdt = inter.Vin(t)*(-Us[i]/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2328+
@inbounds jac[domain.indexes[3],i] += ddnidTdt
2329+
@inbounds @fastmath jac[domain.indexes[4],i] += P/T*ddnidTdt
2330+
end
2331+
22772332
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
22782333
# outlet
22792334
# flow = inter.Vout(t)*sum(ns)/V
@@ -2503,6 +2558,22 @@ end
25032558
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
25042559
end
25052560
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2561+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2562+
# dn/dt .+= inter.Vin(t)*inter.cs
2563+
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
2564+
# dV/dt += inter.Vin(t)
2565+
# d/dn(dn/dt) = 0
2566+
# d/dV(dn/dt) = 0
2567+
# d/dni(dT/dt) += inter.Vin(t)*(-Hs[i]/N)/(N*Cpave) - dT/dt*(dCpavedni/Cpave)
2568+
# d/dV(dT/dt) += inter.Vin(t)*(dot(Hs, ns)/V^2)/(N*Cpave) - (dT/dt)/(N*Cpave) * d/dV((PV/RT)*Cpave)
2569+
# Note: Cpave = dot(cpdivR,ns)*R/N = dot(cpdivR,ns)*R*(RT/PV)
2570+
# Note: d/dV(Cpave) = dot(cpdivR,ns)*R*(RT/PV)*(-1/V) = -Cpave/V
2571+
# = inter.Vin(t)*(dot(Hs, ns)/V^2)/(N*Cpave) - dT/dt/(N*Cpave) * ((P/RT)*Cpave - (PV/RT)*Cpave/V)
2572+
# = 0
2573+
@simd for i in domain.indexes[1]:domain.indexes[2]
2574+
@inbounds @fastmath dCpavedni = cpdivR[i]*R/N
2575+
@inbounds jac[domain.indexes[3],i] += inter.Vin(t)*(-Hs[i]/N)/(N*Cpave) - dTdt*(dCpavedni/Cpave)
2576+
end
25062577
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
25072578
# outlet
25082579
# dTdt = 0
@@ -2587,6 +2658,10 @@ end
25872658
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
25882659
end
25892660
@views @inbounds @fastmath jac[domain.indexes[3],domain.indexs[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2661+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2662+
# dn/dt .+= inter.Vin(t)*inter.cs
2663+
# dV/dt += inter.Vin(t)
2664+
nothing
25902665
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
25912666
# outlet
25922667
# d/dni(dV/dt) = dflow/dni*R*T/P = inter.Vout(t)/V*R*T/P = inter.Vout(t)/N
@@ -2694,6 +2769,24 @@ end
26942769
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
26952770
@inbounds @fastmath jac[domain.indexes[4],i] -= kLAs[i]/kHs[i]*R*T/V + P/T*ddnidTdt
26962771
end
2772+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2773+
# dn/dt .+= inter.Vin(t)*inter.cs
2774+
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)
2775+
# dP/dt += inter.Vin(t)*P/V + dT/dt*P/T
2776+
# d/dn(dn/dt) += 0
2777+
# d/dni(dT/dt) += inter.Vin(t)*(-Us[i]/V)/(N*Cvave) - inter.Vin(t)*(inter.Hpervolume - dot(Us,ns)/V)/(N*Cvave)^2 * d/dni(N*Cvave)
2778+
# = inter.Vin(t)*(-Us[i]/V)/(N*Cvave) - dT/dt * d/dni(N*Cvave) / (N*Cvave)
2779+
# = inter.Vin(t)*(-Us[i]/V)/(N*Cvave) - dT/dt * dCvave/dni / Cvave
2780+
# Note: Cvave = sum(ns.*cpdivR)*R/N - R
2781+
# Note: dCvave/dni = cpdivR[i]*R/N
2782+
# d/dni(dP/dt) += d/dni(dT/dt)*P/T
2783+
@simd for i in domain.indexes[1]:domain.indexes[2]
2784+
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2785+
@fastmath ddnidTdt = (inter.Vout*P/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2786+
ddnidTdt = inter.Vin(t)*(-Us[i]/V)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2787+
@inbounds jac[domain.indexes[3],i] += ddnidTdt
2788+
@inbounds @fastmath jac[domain.indexes[4],i] += ddnidTdt*P/T
2789+
end
26972790
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
26982791
# outlet
26992792
# flow = inter.Vout(t)*sum(ns)/V
@@ -2826,6 +2919,27 @@ end
28262919
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
28272920
end
28282921
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2922+
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
2923+
# dn/dt .+= inter.Vin(t)*inter.cs
2924+
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
2925+
# dV/dt += inter.Vin(t)
2926+
# d/dn(dn/dt) += 0
2927+
# d/dV(dn/dt) += 0
2928+
# d/dni(dT/dt) += inter.V(t)*(-Hs[i]/V)/(N*Cpave) - dT/dt/(N*Cpave) * d/dni(N*Cpave)
2929+
# Note: Cpave = dot(cpdivR,ns)*R/N
2930+
# Note: d/dni(Cpave) = cpdivR[i]*R/N
2931+
# = inter.V(t)*(-Hs[i]/V)/(N*Cpave) - dT/dt * (dCpavedni/Cpave)
2932+
# d/dV(dT/dt) += inter.Vin(t)*(dot(Hs,ns)/V^2)/(N*Cpave) - dT/dt/(N*Cpave) * d/dV(PV/RT*Cpave)
2933+
# = inter.Vin(t)*(dot(Hs,ns)/V^2)/(N*Cpave) - dT/dt/(N*Cpave) * (P/RT*Cpave + PV/RT*dCpavedV)
2934+
# Note: Cpave = dot(cpdivR,ns)*R/N = dot(cpdivR,ns)*R*(RT/PV)
2935+
# Note: d/dV(Cpave) = dot(cpdivR,ns)*R*(RT/PV)*(-1/V) = -Cpave/V
2936+
# = inter.Vin(t)*(dot(Hs,ns)/V^2)/(N*Cpave) - dT/dt/(N*Cpave) * (P/RT*Cpave + PV/RT(-Cpave/V))
2937+
# = 0
2938+
dTdt = inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
2939+
@simd for i in domain.indexes[1]:domain.indexes[2]
2940+
@inbounds @fastmath dCpavedni = cpdivR[i]*R/N
2941+
@inbounds @fastmath jac[domain.indexes[3],i] += inter.Vin(t)*(-Hs[i]/V)/(N*Cpave) - dTdt*(dCpavedni/Cpave)
2942+
end
28292943
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
28302944
# outlet
28312945
# dTdt = 0

0 commit comments

Comments
 (0)