Skip to content

Commit 438eac2

Browse files
hwpanghwpang
committed
Update Jacobian of kLAkH interface based on new derivation
1 parent 5a46616 commit 438eac2

File tree

1 file changed

+41
-25
lines changed

1 file changed

+41
-25
lines changed

‎src/Domain.jl

Lines changed: 41 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -2092,9 +2092,13 @@ end
20922092
# -= d/dni(dni/dt)*R*T/P
20932093
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
20942094
@simd for i in domain.indexes[1]:domain.indexes[2]
2095-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2095+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2096+
@inbounds @fastmath jac[i,i] -= ddnidnidt
2097+
@inbounds @fastmath jac[domain.indexes[3],i] -= ddnidnidt*R*T/P
20962098
end
2097-
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= kLAs./kHs*R*T/P
2099+
ddVdndt = kLAs.*inter.V.*cs*(-1/V)*R*T./kHs
2100+
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= ddVdndt
2101+
@views @inbounds @fastmath jac[domain.indexes[3],domain.indexes[3]] -= sum(ddVdndt)*R*T/P
20982102
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
20992103
# dn/dt += inter.Vin(t)*inter.cs
21002104
# dV/dt += inter.Vin(t)
@@ -2226,11 +2230,12 @@ end
22262230
flow = sum(cond)
22272231
@fastmath dTdt = (P*V/N*flow)/(N*Cvave)
22282232
@simd for i in domain.indexes[1]:domain.indexes[2]
2229-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2233+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2234+
@inbounds @fastmath jac[i,i] -= ddnidnidt
22302235
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2231-
@fastmath ddnidTdt = (P*V/N*kLAs[i]/kHs[i])/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2236+
@fastmath ddnidTdt = (P*V/N*ddnidnidt)/(N*Cvave) - dTdt*(dCvavedni/Cvave)
22322237
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
2233-
@inbounds @fastmath jac[domain.indexes[4],i] -= kLAs[i]/kHs[i]*R*T/V + P/T*ddnidTdt
2238+
@inbounds @fastmath jac[domain.indexes[4],i] -= ddnidnidt*R*T/V + P/T*ddnidTdt
22342239
end
22352240
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
22362241
# dn/dt += inter.Vin(t)*inter.cs
@@ -2398,8 +2403,7 @@ end
23982403
end
23992404
@fastmath ddVdTdt = flow*H/V/(N*Cpave)
24002405
@inbounds jac[domain.indexes[3],domain.indexes[4]] += ddVdTdt
2401-
@inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .+= kLAs.*inter.cs
2402-
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += flow/N + dTdt/T + V/T*ddVdTdt
2406+
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += dTdt/T + V/T*ddVdTdt
24032407

24042408
# condensation
24052409
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
@@ -2413,9 +2417,13 @@ end
24132417
# -= ddnidnidt*R*T/P
24142418
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
24152419
@simd for i in domain.indexes[1]:domain.indexes[2]
2416-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2420+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2421+
@inbounds @fastmath jac[i,i] -= ddnidnidt
2422+
@inbounds @fastmath jac[domain.indexes[4],i] -= ddnidnidt*R*T/P
24172423
end
2418-
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2424+
ddVdndt = kLAs.*inter.V.*cs*(-1/V)*R*T./kHs
2425+
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .-= ddVdndt
2426+
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] .-= sum(ddVdndt)*R*T/P
24192427
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
24202428
# dn/dt .+= inter.Vin(t)*inter.cs
24212429
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)
@@ -2501,12 +2509,12 @@ end
25012509
kHs = map.(inter.kHs,inter.T)
25022510

25032511
# evaporaiton
2504-
# inlet
2505-
# d/dV (dni/dt) = dflow_i/dV = -kLAs[i]*inter.cs[i]/V^2
2506-
# flow = sum(kLAs.*inter.cs)/V
2507-
# flow_i = kLAs[i]*inter.cs[i]/V
2508-
# dflow_idV = -kLAs[i]*inter.cs[i]/V^2
2509-
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .+= -kLAs.*inter.cs/(V*V)
2512+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2513+
# dV/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/P
2514+
# d/dni(dni/dt) += 0
2515+
# d/dV(dni/dt) += 0
2516+
# d/dni(dV/dt) += 0
2517+
# d/dV(dV/dt) += 0
25102518

25112519
# condensation
25122520
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
@@ -2517,9 +2525,13 @@ end
25172525
# -= d/dni(dni/dt)*R*T/P
25182526
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
25192527
@simd for i in domain.indexes[1]:domain.indexes[2]
2520-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2528+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2529+
@inbounds @fastmath jac[i,i] -= ddnidnidt
2530+
@inbounds @fastmath jac[domain.indexes[3],i] -= ddnidnidt*R*T/P
25212531
end
2522-
@views @inbounds @fastmath jac[domain.indexes[3],domain.indexs[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2532+
ddVdndt = kLAs.*inter.V.*cs.*(-1/V)*R*T./kHs
2533+
@views @inbounds @fastmath jac[domain.indexs[1]:domain.indexes[2],domain.indexes[3]] .-= ddVdndt
2534+
@views @inbounds @fastmath jac[domain.indexes[3],domain.indexes[3]] -= sum(ddVdndt)*R*T/P
25232535
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
25242536
# dn/dt .+= inter.Vin(t)*inter.cs
25252537
# dV/dt += inter.Vin(t)
@@ -2638,11 +2650,12 @@ end
26382650
flow = sum(cond)
26392651
@fastmath dTdt = (P*V/N*flow)/(N*Cvave)
26402652
@simd for i in domain.indexes[1]:domain.indexes[2]
2641-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2653+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2654+
@inbounds @fastmath jac[i,i] -= ddnidnidt
26422655
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2643-
@fastmath ddnidTdt = (P*V/N*kLAs[i]/kHs[i])/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2656+
@fastmath ddnidTdt = (P*V/N*ddnidnidt)/(N*Cvave) - dTdt*(dCvavedni/Cvave)
26442657
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
2645-
@inbounds @fastmath jac[domain.indexes[4],i] -= kLAs[i]/kHs[i]*R*T/V + P/T*ddnidTdt
2658+
@inbounds @fastmath jac[domain.indexes[4],i] -= ddnidnidt*R*T/V + P/T*ddnidTdt
26462659
end
26472660
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
26482661
# dn/dt .+= inter.Vin(t)*inter.cs
@@ -2793,14 +2806,13 @@ end
27932806
@fastmath dTdt = flow*(inter.H - H)/(N*Cpave)
27942807
@simd for i in domain.indexes[1]:domain.indexes[2]
27952808
@inbounds @fastmath dCpavedni = cpdivR[i]*R/N
2796-
@inbounds @fastmath ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave)-dTdt*(dCpavedni/Cpave)
2809+
@inbounds @fastmath ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave) - dTdt*(dCpavedni/Cpave)
27972810
@inbounds jac[domain.indexes[3],i] += ddnidTdt
27982811
@inbounds @fastmath jac[domain.indexes[4],i] += V/T*ddnidTdt
27992812
end
28002813
@fastmath ddVdTdt = flow*H/V/(N*Cpave)
28012814
@inbounds jac[domain.indexes[3],domain.indexes[4]] += ddVdTdt
2802-
@inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .+= kLAs.*inter.cs
2803-
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += flow/N + dTdt/T + V/T*ddVdTdt
2815+
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += dTdt/T + V/T*ddVdTdt
28042816

28052817
# condensation
28062818
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
@@ -2814,9 +2826,13 @@ end
28142826
# -= d/dni(dni/dt)*R*T/P
28152827
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
28162828
@simd for i in domain.indexes[1]:domain.indexes[2]
2817-
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
2829+
ddnidnidt = kLAs[i]*inter.V/V*R*T/kHs[i]
2830+
@inbounds @fastmath jac[i,i] -= ddnidnidt
2831+
@inbounds @fastmath jac[domain.indexes[4],i] -= ddnidnidt*R*T/P
28182832
end
2819-
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= kLAs./kHs*R*T/P
2833+
ddVdndt = kLAs.*inter.V.*cs.*(-1/V)*R*T./kHs
2834+
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .-= ddVdndt
2835+
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] -= sum(ddVdndt)*R*T/P
28202836
elseif isa(inter,VolumetricFlowRateInlet) && domain == inter.domain
28212837
# dn/dt .+= inter.Vin(t)*inter.cs
28222838
# dT/dt += inter.Vin(t)*(inter.Hpervolume - dot(Hs,ns)/V)/(N*Cpave)

0 commit comments

Comments
 (0)