Skip to content

Commit 0bfd28e

Browse files
hwpanghwpang
committed
Rewrite and debug Jacobian expression for domains and interfaces
1 parent 1d5551e commit 0bfd28e

File tree

1 file changed

+85
-89
lines changed

1 file changed

+85
-89
lines changed

‎src/Domain.jl

Lines changed: 85 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -2105,21 +2105,17 @@ end
21052105
nothing
21062106

21072107
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2108-
# outlet
2109-
# d/dni(dV/dt) = dflow/dni*R*T/P = inter.Vout(t)/V*R*T/P = inter.Vout(t)/N
2110-
# d/dV(dV/dt) = dflow/dV *R*T/P = -inter.Vout(t)*sum(ns)/V^2 *R*T/P = -inter.Vout(t)/V*sum(ns)/N
2111-
# d/dV(dni/dt) = dflow_i/dV
2112-
# flow = inter.Vout(t)/V*sum(ns)
2113-
# dflow/dni = inter.Vout(t)/V
2114-
# dflow/dV = -inter.Vout(t)*sum(ns)/V^2
2115-
# flow_i = inter.Vout(t)/V*ns[i]
2116-
# dflow_i/dV = -inter.Vout(t)*ns[i]/V^2
2108+
# dn/dt .-= inter.Vout(t)*ns/V
2109+
# dV/dt -= inter.Vout(t)
2110+
# d/dni(dni/dt) -= inter.Vout(t)/V
2111+
# d/dV(dn/dt) -= -inter.Vout(t)*ns/(V*V)
2112+
# d/dni(dV/dt) -= 0
2113+
# d/dV(dV/dt) -= 0
2114+
21172115
@simd for i in domain.indexes[1]:domain.indexes[2]
21182116
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
21192117
end
2120-
@inbounds @fastmath jac[domain.indexes[3],domain.indexes[1]:domain.indexes[2]] .-= inter.Vout(t)/N
21212118
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= -inter.Vout(t)*ns/(V*V)
2122-
@inbounds @fastmath jac[domain.indexes[3],domain.indexes[3]] -= -inter.Vout(t)/V*sum(ns)/N
21232119
end
21242120
end
21252121
end
@@ -2243,19 +2239,23 @@ end
22432239
end
22442240

22452241
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2246-
# outlet
2247-
# flow = inter.Vout(t)*sum(ns)/V
2248-
# dflowdni = inter.Vout(t)/V
2249-
# dTdt = flow*(P*V/N)/(N*Cvave)
2250-
# ddnidTdt = ( dflowdni *P*V/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave) = (inter.Vout(t)/V*P*V/N)/(N*Cvave) -dTdt*(dCvavedni/Cvave)
2251-
# d/dni (dP/dt) = dflowdni *R*T/V + P/T * d/dni(dT/dt) = inter.Vout(t)/V*R*T/V + P/T * d/dni(dT/dt) =
2242+
# dn/dt .-= inter.Vout(t)*ns/V
2243+
# dT/dt -= (P*inter.Vout(t))/(N*Cvave)
2244+
# dP/dt -= inter.Vout(t)*P/V + P/T*dTdt
2245+
# d/dni(dni/dt) -= inter.Vout(t)/V
2246+
# d/dni(dT/dt) -= dT/dt * (-1/(N*Cvave)) * d/dni(N*Cvave)
2247+
# -= dT/dt * (-1/(Cvave)) * d/dni(Cvave)
2248+
# Note: Cvave = dot(cpdivR,ns)*R/N-R
2249+
# Note: d/dni(Cvave) = cpdivR[i]*R/N
2250+
# -= -dT/dt * (dCvavedni/Cvave)
2251+
# d/dni(dP/dt) -= P/T * d/dni(dT/dt)
22522252
@fastmath dTdt = (P*inter.Vout(t))/(N*Cvave)
22532253
@simd for i in domain.indexes[1]:domain.indexes[2]
22542254
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
22552255
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2256-
@fastmath ddnidTdt = (inter.Vout*P/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2256+
@fastmath ddnidTdt = -dTdt*(dCvavedni/Cvave)
22572257
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
2258-
@inbounds @fastmath jac[domain.indexes[4],i] -= inter.Vout(t)/V*R*T/V + P/T*ddnidTdt
2258+
@inbounds @fastmath jac[domain.indexes[4],i] -= P/T*ddnidTdt
22592259
end
22602260
end
22612261
end
@@ -2305,18 +2305,22 @@ end
23052305

23062306
@simd for inter in interfaces
23072307
if isa(inter,Inlet) && domain == inter.domain
2308-
# inlet
2309-
# dTdt = flow*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2310-
# ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave)-dTdt*(dCpavedni/Cpave)
2311-
# d/dni (dV/dt) = V/T * d/dni(dT/dt)
2312-
# d/dV (dT/dt) = flow*(dot(Hs, ns)/N)/V/(N*Cpave)
2313-
# d/dV (dV/dt) = dflow/dV*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = dT/dt/T + V/T * d/dV(dT/dt)
2314-
# d/dV(dni/dt) = dflow_i/dV = -flow*ns[i]/N/V
2315-
# dflow/dV = 0
2316-
# flow_i = flow*ns[i]/N
2317-
# N = P*V/(R*T)
2318-
# dN/dV = P/(R*T)
2319-
# dflow_i/dV = -flow*ns[i]/N^2 P/(R*T) = -flow*ns[i]/N/V
2308+
# dn/dt .+= inter.y.*inter.F(t)
2309+
# dT/dt += inter.F(t)*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2310+
# dV/dt += inter.F(t)*R*T/P + dT/dt*V/T
2311+
# d/dni(dn/dt) += 0
2312+
# d/dV(dn/dt) += 0
2313+
# d/dni(dT/dt) += inter.F(t)*(-Hs[i]/N)/(N*Cpave) - dTdt/(N*Cpave) * d/dni(N*Cpave)
2314+
# += inter.F(t)*(-Hs[i]/N)/(N*Cpave) - dTdt/(Cpave) * d/dni(Cpave)
2315+
# Note: Cpave = dot(cpdivR,ns)*R/N-R
2316+
# Note: d/dni(Cpave) = cpdivR[i]*R/N
2317+
# += inter.F(t)*(-Hs[i]/N)/(N*Cpave) - dTdt * (dCpavedni/Cpave)
2318+
# d/dV(dT/dt) += inter.F(t)*(dot(Hs,ns)/N^2*dN/dV)/(N*Cpave) + dT/dt * (-1/(N*Cpave)) * d/dV(N*Cpave)
2319+
# += inter.F(t)*(dot(Hs,ns)/N^2*P/RT)/(N*Cpave) + 0
2320+
# += inter.F(t)*(dot(Hs,ns)/N^2*P/RT)/(N*Cpave)
2321+
# += inter.F(t)*(dot(Hs,ns)/N/V)/(N*Cpave)
2322+
# d/dni(dV/dt) += d/dni(dT/dt)*V/T
2323+
# d/dV(dV/dt) += d/dV(dT/dt)*V/T
23202324
flow = inter.F(t)
23212325
@fastmath H = dot(Hs,ns)/N
23222326
@fastmath dTdt = flow*(inter.H - H)/(N*Cpave)
@@ -2328,8 +2332,7 @@ end
23282332
end
23292333
@fastmath ddVdTdt = flow*H/V/(N*Cpave)
23302334
@inbounds jac[domain.indexes[3],domain.indexes[4]] += ddVdTdt
2331-
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += dTdt/T + V/T*ddVdTdt
2332-
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .+= flow*ns/N/V
2335+
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += V/T*ddVdTdt
23332336
elseif isa(inter,Outlet) && domain == inter.domain
23342337
# dn/dt .-= inter.F(t).*ns./N
23352338
# dT/dt -= 0
@@ -2408,22 +2411,19 @@ end
24082411
@inbounds jac[domain.indexes[3],i] += inter.Vin(t)*(-Hs[i]/N)/(N*Cpave) - dTdt*(dCpavedni/Cpave)
24092412
end
24102413
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2411-
# outlet
2412-
# dTdt = 0
2413-
# d/dni (dV/dt) = dflow/dni *R*T/P = inter.Vout/V*R*T/P = inter.Vout/N
2414-
# d/dV(dV/dt) = dflowdV *R*T/P = -inter.Vout*sum(ns)/V^2*R*T/P = -inter.Vout/V*sum(ns)/N
2415-
# d/dV(dni/dt) = dflow_i/dV = -inter.Vout*ns[i]/V^2
2416-
# dflow/dni = inter.Vout/V
2417-
# dflowdV = -inter.Vout*sum(ns)/V^2
2418-
# dflow/dV = 0
2419-
# flow_i = inter.Vout*ns[i]/V
2420-
# dflow_i/dV = -inter.Vout*ns[i]/V^2
2414+
# dn/dt .-= inter.Vout(t)*ns/V
2415+
# dT/dt -= 0
2416+
# dV/dt -= inter.Vout(t)
2417+
# d/dni(dni/dt) -= inter.Vout(t)/V
2418+
# d/dV(dn/dt) -= -inter.Vout(t)*ns/(V^2)
2419+
# d/dni(dT/dt) -= 0
2420+
# d/dV(dT/dt) -= 0
2421+
# d/dni(dV/dt) -= 0
2422+
# d/dV(dV/dt) -= 0
24212423
@simd for i in domain.indexes[1]:domain.indexes[2]
24222424
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
24232425
end
2424-
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= inter.Vout(t)/N
24252426
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .-= -inter.Vout(t)/(V*V)*ns
2426-
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] -= -inter.Vout(t)/V*sum(ns)/N
24272427
end
24282428
end
24292429

@@ -2506,21 +2506,16 @@ end
25062506
# dV/dt += inter.Vin(t)
25072507
nothing
25082508
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2509-
# outlet
2510-
# d/dni(dV/dt) = dflow/dni*R*T/P = inter.Vout(t)/V*R*T/P = inter.Vout(t)/N
2511-
# d/dV(dV/dt) = dflow/dV *R*T/P = -inter.Vout*sum(ns)/V^2*R*T/P = -inter.Vout/V*sum(ns)/N
2512-
# d/dV(dni/dt) = dflow_i/dV = -inter.Vout(t)*ns[i]/V^2
2513-
# flow = inter.Vout(t)/V*sum(ns)
2514-
# dflowdni = inter.Vout(t)/V
2515-
# dflowdV = -inter.Vout*sum(ns)/V^2
2516-
# dflow_i/dV = -inter.Vout(t)*ns[i]/V^2
2509+
# dn/dt .-= inter.Vout(t)*ns/V
2510+
# dV/dt -= inter.Vout(t)
2511+
# d/dni(dni/dt) -= inter.Vout(t)/V
2512+
# d/dV(dni/dt) -= -inter.Vout(t)*ns/(V*V)
2513+
# d/dni(dV/dt) -= 0
2514+
# d/dV(dV/dt) -= 0
25172515
@simd for i in domain.indexes[1]:domain.indexes[2]
25182516
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
25192517
end
25202518
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .-= -inter.Vout(t)*ns/(V*V)
2521-
@views @inbounds @fastmath jac[domain.indexes[3],domain.indexes[1]:domain.indexes[2]] .-= inter.Vout(t)/N
2522-
@inbounds @fastmath jac[domain.indexes[3],domain.indexes[3]] -= -inter.Vout/V*sum(ns)/N
2523-
25242519
end
25252520
end
25262521

@@ -2646,19 +2641,21 @@ end
26462641
@inbounds @fastmath jac[domain.indexes[4],i] += ddnidTdt*P/T
26472642
end
26482643
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2649-
# outlet
2650-
# flow = inter.Vout(t)*sum(ns)/V
2651-
# dflowdni = inter.Vout(t)/V
2652-
# dTdt = flow*(P*V/N)/(N*Cvave)
2653-
# ddnidTdt = ( dflowdni *P*V/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave) = (inter.Vout(t)/V*P*V/N)/(N*Cvave) -dTdt*(dCvavedni/Cvave)
2654-
# d/dni (dP/dt) = dflowdni *R*T/V + P/T * d/dni(dT/dt) = inter.Vout(t)/V*R*T/V + P/T * d/dni(dT/dt) =
2644+
# dn/dt .-= inter.Vout(t)*ns/V
2645+
# dT/dt -= (P*inter.Vout(t))/(N*Cvave)
2646+
# dP/dt -= inter.Vout(t)*P/V + dT/dt*P/T
2647+
# d/dni(dni/dt) -= inter.Vout(t)/V
2648+
# d/dni(dT/dt) -= dT/dt *(-1/(N*Cvave)) d/dni(N*Cvave)
2649+
# -= dT/dt *(-1/(Cvave)) d/dni(Cvave)
2650+
# -= - dT/dt * (dCvavedni/Cvave)
2651+
# d/dni(dP/dt) -= P/T * d/dni(dT/dt)
26552652
@fastmath dTdt = (P*inter.Vout(t))/(N*Cvave)
26562653
@simd for i in domain.indexes[1]:domain.indexes[2]
26572654
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
26582655
@inbounds @fastmath dCvavedni = cpdivR[i]*R/N
2659-
@fastmath ddnidTdt = (inter.Vout*P/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2656+
@fastmath ddnidTdt = -dTdt*(dCvavedni/Cvave)
26602657
@inbounds jac[domain.indexes[3],i] -= ddnidTdt
2661-
@inbounds @fastmath jac[domain.indexes[4],i] -= inter.Vout(t)/V*R*T/V + P/T*ddnidTdt
2658+
@inbounds @fastmath jac[domain.indexes[4],i] -= P/T*ddnidTdt
26622659
end
26632660
end
26642661
end
@@ -2702,18 +2699,21 @@ end
27022699

27032700
@simd for inter in interfaces
27042701
if isa(inter,Inlet) && domain == inter.domain
2705-
# inlet
2706-
# dTdt = flow*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2707-
# ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave)-dTdt*(dCpavedni/Cpave)
2708-
# d/dni (dV/dt) = V/T * d/dni(dT/dt)
2709-
# d/dV (dT/dt) = flow*(dot(Hs, ns)/N)/V/(N*Cpave)
2710-
# d/dV (dV/dt) = dflow/dV*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = dT/dt/T + V/T * d/dV(dT/dt)
2711-
# d/dV(dni/dt) = dflow_i/dV = -flow*ns[i]/N/V
2712-
# dflow/dV = 0
2713-
# flow_i = flow*ns[i]/N
2714-
# N = P*V/(R*T)
2715-
# dN/dV = P/(R*T)
2716-
# dflow_i/dV = -flow*ns[i]/N^2 P/(R*T) = -flow*ns[i]/N/V
2702+
# dn/dt .+= inter.y.*inter.F(t)
2703+
# dT/dt += inter.F(t)*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2704+
# dV/dt += inter.F(t)*R*T/P + dT/dt*V/T
2705+
# d/dni(dni/dt) += 0
2706+
# d/dV(dni/dt) += 0
2707+
# d/dni(dT/dt) += inter.F(t)*(-Hs[i]/N)/(N*Cpave) + dTdt * (-1/(N*Cpave)) * d/dni(N*Cpave)
2708+
# += inter.F(t)*(-Hs[i]/N)/(N*Cpave) + dTdt * (-1/(Cpave)) * d/dni(Cpave)
2709+
# Note: Cpave = = dot(cpdivR,ns)*R/N
2710+
# Note: dCpavedni = cpdivR[i]*R/N
2711+
# d/dV(dT/dt) += inter.F(t)*(dot(Hs,ns)/N^2*dN/dV)/(N*Cpave) + dTdt * (-1/(N*Cpave)) * d/dV(N*Cpave)
2712+
# Note: dN/dV = d(PV/RT)/dV = P/RT
2713+
# Note: d(N*Cpave)/dV = (dN/dV*Cpave + N*dCpave/dV) = P/RT*Cpave - Cpave * P/RT = 0
2714+
# += inter.F(t)*(dot(Hs,ns)/N/V)/(N*Cpave)
2715+
# d/dni(dV/dt) += V/T*d/dni(dT/dt)
2716+
# d/dV(dV/dt) += d/dV(dT/dt)*V/T + dT/dt/T
27172717
flow = inter.F(t)
27182718
@fastmath H = dot(Hs,ns)/N
27192719
@fastmath dTdt = flow*(inter.H - H)/(N*Cpave)
@@ -2726,7 +2726,6 @@ end
27262726
@fastmath ddVdTdt = flow*H/V/(N*Cpave)
27272727
@inbounds jac[domain.indexes[3],domain.indexes[4]] += ddVdTdt
27282728
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += dTdt/T + V/T*ddVdTdt
2729-
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .+= flow*ns/N/V
27302729
elseif isa(inter,Outlet) && domain == inter.domain
27312730
# dn/dt .-= inter.F(t).*ns./N
27322731
# dT/dt -= 0
@@ -2810,22 +2809,19 @@ end
28102809
@inbounds @fastmath jac[domain.indexes[3],i] += inter.Vin(t)*(-Hs[i]/V)/(N*Cpave) - dTdt*(dCpavedni/Cpave)
28112810
end
28122811
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
2813-
# outlet
2814-
# dTdt = 0
2815-
# d/dni (dV/dt) = dflow/dni *R*T/P = inter.Vout/V*R*T/P = inter.Vout/N
2816-
# d/dV(dV/dt) = dflowdV *R*T/P = -inter.Vout*sum(ns)/V^2*R*T/P = -inter.Vout/V*sum(ns)/N
2817-
# d/dV(dni/dt) = dflow_i/dV = -inter.Vout*ns[i]/V^2
2818-
# dflow/dni = inter.Vout/V
2819-
# dflowdV = -inter.Vout*sum(ns)/V^2
2820-
# dflow/dV = 0
2821-
# flow_i = inter.Vout*ns[i]/V
2822-
# dflow_i/dV = -inter.Vout*ns[i]/V^2
2812+
# dn/dt .-= inter.Vout(t)*ns/V
2813+
# dT/dt -= 0
2814+
# dV/dt -= inter.Vout(t)
2815+
# d/dni(dni/dt) -= inter.Vout(t)/V
2816+
# d/dV(dni/dt) -= -inter.Vout(t)*ns/V^2
2817+
# d/dni(dT/dt) -= 0
2818+
# d/dV(dT/dt) -= 0
2819+
# d/dni(dV/dt) -= 0
2820+
# d/dV(dV/dt) -= 0
28232821
@simd for i in domain.indexes[1]:domain.indexes[2]
28242822
@inbounds @fastmath jac[i,i] -= inter.Vout(t)/V
28252823
end
2826-
@views @inbounds @fastmath jac[domain.indexes[4],domain.indexes[1]:domain.indexes[2]] .-= inter.Vout(t)/N
28272824
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[4]] .-= -inter.Vout(t)/(V*V)*ns
2828-
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] -= -inter.Vout(t)/V*sum(ns)/N
28292825
end
28302826
end
28312827

0 commit comments

Comments
 (0)