Skip to content

Commit 5a46616

Browse files
committed
Update comments of Jacobian derivation for kLAkH interface
1 parent 1746dae commit 5a46616

File tree

1 file changed

+121
-86
lines changed

1 file changed

+121
-86
lines changed

src/Domain.jl

Lines changed: 121 additions & 86 deletions
Original file line numberDiff line numberDiff line change
@@ -2076,22 +2076,21 @@ end
20762076
kHs = map.(inter.kHs,inter.T)
20772077

20782078
# evaporation
2079-
# inlet
2080-
# d/dV(dni/dt) = dflow_i/dV
2081-
# flow = sum(inter.kLAs.*inter.cs)/V
2082-
# flow_i = inter.kLAs[i]*inter.cs[i]/V
2083-
# d/dV(dni/dt) = dflow_i/dV = -inter.kLAs[i]*inter.cs[i]/(V*V)
2084-
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .+= -kLAs.*inter.cs/(V*V)
2079+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2080+
# dV/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/P
2081+
# d/dni(dn/dt) += 0
2082+
# d/dV(dn/dt) += 0
2083+
# d/dni(dV/dt) += 0
2084+
# d/dV(dV/dt) += 0
20852085

20862086
# condensation
2087-
# outlet
2088-
# d/dni(dV/dt) = dflow/dni*R*T/P = kLAs[i]/kHs[i]*R*T/P
2089-
# d/dV(dV/dt) = dflow/dV *R*T/P = 0
2090-
# d/dV(dni/dt) = dflow_i/dV = 0
2091-
# flow = sum(kLAs.*ns./kHs)
2092-
# dflow/dni = kLAs[i]/kHs[i]
2093-
# dflow/dV = 0
2094-
# dflow_i/dV = 0
2087+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2088+
# dV/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/P
2089+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2090+
# d/dV(dn/dt) -= kLAs.*inter.V.*cs*(-1/V)*R*T./kHs
2091+
# d/dni(dV/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/P
2092+
# -= d/dni(dni/dt)*R*T/P
2093+
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
20952094
@simd for i in domain.indexes[1]:domain.indexes[2]
20962095
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
20972096
end
@@ -2191,11 +2190,17 @@ end
21912190
evap = kLAs.*inter.V.*inter.cs
21922191
cond = kLAs.*inter.V.*cs*R*T./kHs
21932192

2194-
#evaporation
2195-
# inlet
2196-
# dTdt = flow*(inter.H - dot(Us,ns)/N)/(N*Cvave)
2197-
# ddnidTdt = flow*(-Us[i]/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2198-
# d/dni (dP/dt) = P/T * d/dni(dT/dt)
2193+
# evaporation
2194+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2195+
# dT/dt += sum(kLAs.*inter.V.*inter.cs)*(inter.H - dot(Us,ns)/N)/(N*Cvave)
2196+
# dP/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/V + P/T*dT/dt
2197+
# d/dni(dni/dt) += 0
2198+
# d/dni(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(-Us[i]/N)/(N*Cvave) - dT/dt * (-1/(N*Cvave)) * d/dni(N*Cvave)
2199+
# += sum(kLAs.*inter.V.*inter.cs)*(-Us[i]/N)/(N*Cvave) - dT/dt * (-1/(Cvave)) * d/dni(Cvave)
2200+
# Note: Cvave = dot(cpdivR,ns)*R/N-R
2201+
# Note: d/dni(Cvave) = cpdivR[i]*R/N
2202+
# += sum(kLAs.*inter.V.*inter.cs)*(-Us[i]/N)/(N*Cvave) - dT/dt * (dCvavedni/Cvave)
2203+
# d/dni(dP/dt) += P/T*d/dni(dT/dt)
21992204
flow = sum(evap)
22002205
@fastmath dTdt = flow*(inter.H - dot(Us,ns)/N)/(N*Cvave)
22012206
@simd for i in domain.indexes[1]:domain.indexes[2]
@@ -2206,12 +2211,18 @@ end
22062211
end
22072212

22082213
# condensation
2209-
# outlet
2210-
# flow = sum(kLAs.*ns./kHs)
2211-
# dflowdni = kLAs[i]/kHs[i]
2212-
# dTdt = flow*(P*V/N)/(N*Cvave)
2213-
# ddnidTdt = ( dflowdni *P*V/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave) = (kLAs[i]/kHs[i]*P*V/N)/(N*Cvave) -dTdt*(dCvavedni/Cvave)
2214-
# d/dni (dP/dt) = dflowdni *R*T/V + P/T * d/dni(dT/dt) = kLAs[i]/kHs[i]*R*T/V + P/T * d/dni(dT/dt)
2214+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2215+
# dT/dt -= (P*V/N*sum(kLAs.*inter.V.*cs*R*T./kHs))/(N*Cvave)
2216+
# dP/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/V + P/T*dT/dt
2217+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2218+
# d/dni(dT/dt) -= (P*V/N*kLAs[i]*inter.V/V*R*T/kHs[i])/(N*Cvave) + dT/dt * (-1/(N*Cvave)) * d/dni(N*Cvave)
2219+
# -= (P*V/N*kLAs[i]*inter.V/V*R*T/kHs[i])/(N*Cvave) + dT/dt * (-1/(Cvave)) * d/dni(Cvave)
2220+
# Note: Cvave = dot(cpdivR,ns)*R/N-R
2221+
# Note: d/dni(Cvave) = cpdivR[i]*R/N
2222+
# -= (P*V/N*kLAs[i]*inter.V/V*R*T/kHs[i])/(N*Cvave) - dT/dt * (dCvavedni/Cvave)
2223+
# -= (P*V/N*d/dni(dni/dt))/(N*Cvave) - dT/dt * (dCvavedni/Cvave)
2224+
# d/dni(dP/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/V + P/T*d/dni(dT/dt)
2225+
# -= d/dni(dni/dt)*R*T/V + P/T*d/dni(dT/dt)
22152226
flow = sum(cond)
22162227
@fastmath dTdt = (P*V/N*flow)/(N*Cvave)
22172228
@simd for i in domain.indexes[1]:domain.indexes[2]
@@ -2357,16 +2368,25 @@ end
23572368
evap = kLAs.*inter.V.*inter.cs
23582369

23592370
# evaporation
2360-
# inlet
2361-
# dTdt = flow*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2362-
# ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave)-dTdt*(dCpavedni/Cpave)
2363-
# d/dni (dV/dt) = V/T * d/dni(dT/dt)
2364-
# d/dV (dT/dt) = flow*(dot(Hs, ns)/N)/V/(N*Cpave)
2365-
# d/dV (dV/dt) = dflow/dV*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = flow/V*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = flow/N + dT/dt/T + V/T * d/dV(dT/dt)
2366-
# d/dV(dni/dt) = dflow_i/dV = kLAs[i]*inter.cs[i]
2367-
# dflowdV = sum(kLAs.*inter.cs) = flow/V
2368-
# flow_i = kLAs[i]*inter.cs[i]*V
2369-
# dflow_i/dV = kLAs[i]*inter.cs[i]
2371+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2372+
# dT/dt += sum(kLAs.*inter.V.*inter.cs)*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2373+
# dV/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/P + dT/dt*V/T
2374+
# d/dni(dn/dt) += 0
2375+
# d/dV(dn/dt) += 0
2376+
# d/dni(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(-Hs[i]/N)/(N*Cpave) - dTdt/(N*Cpave) * d/dni(N*Cpave)
2377+
# += sum(kLAs.*inter.V.*inter.cs)*(-Hs[i]/N)/(N*Cpave) - dTdt/(Cpave) * d/dni(Cpave)
2378+
# Note: Cpave = dot(cpdivR,ns)*R/N-R
2379+
# Note: d/dni(Cpave) = cpdivR[i]*R/N
2380+
# += sum(kLAs.*inter.V.*inter.cs)*(-Hs[i]/N)/(N*Cpave) - dTdt * (dCpavedni/Cpave)
2381+
# d/dV(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(dot(Hs,ns)/N^2*dN/dV)/(N*Cpave) + dT/dt * (-1/(N*Cpave)) * d/dV(N*Cpave)
2382+
# Note: dN/dV = d(PV/RT)/dV = P/RT
2383+
# Note: Cpave = dot(cpdivR,ns)*R/N = dot(cpdivR,ns)*R*(RT/PV)
2384+
# Note: d/dV(Cpave) = dot(cpdivR,ns)*R*(RT/PV)*(-1/V) = -Cpave/V
2385+
# Note: d/dV(N*Cpave) = dN/dV*Cpave + N*dCpave/dV = P/RT*Cpave - (PV/RT)*Cpave/V = 0
2386+
# += sum(kLAs.*inter.V.*inter.cs)*(dot(Hs,ns)/N^2*P/RT)/(N*Cpave) + 0
2387+
# += sum(kLAs.*inter.V.*inter.cs)*(dot(Hs,ns)/N/V)/(N*Cpave)
2388+
# d/dni(dV/dt) += d/dni(dT/dt)*V/T
2389+
# d/dV(dV/dt) += d/dV(dT/dt)*V/T + dT/dt/T
23702390
flow = sum(evap)
23712391
@fastmath H = dot(Hs,ns)/N
23722392
@fastmath dTdt = flow*(inter.H - H)/(N*Cpave)
@@ -2382,15 +2402,16 @@ end
23822402
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += flow/N + dTdt/T + V/T*ddVdTdt
23832403

23842404
# condensation
2385-
# outlet
2386-
# dTdt = 0
2387-
# d/dni (dV/dt) = dflow/dni *R*T/P = kLAs[i]/kHs[i]*R*T/P
2388-
# d/dV(dV/dt) = dflowdV *R*T/P = 0
2389-
# d/dV(dni/dt) = dflow_i/dV = 0
2390-
# flow_i = kLAs[i]*ns[i]/kHs[i]
2391-
# dflow/dni = kLAs[i]/kHs[i]
2392-
# dflow/dV = 0
2393-
# dflow_i/dV = 0
2405+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2406+
# dT/dt -= 0
2407+
# dV/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/P
2408+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2409+
# d/dV(dn/dt) .-= kLAs.*inter.V.*cs*(-1/V)*R*T./kHs
2410+
# d/dni(dT/dt) -= 0
2411+
# d/dV(dT/dt) -= 0
2412+
# d/dni(dV/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/P
2413+
# -= ddnidnidt*R*T/P
2414+
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
23942415
@simd for i in domain.indexes[1]:domain.indexes[2]
23952416
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
23962417
end
@@ -2464,7 +2485,7 @@ end
24642485
# dn/dt .-= inter.F(t)*ns./N
24652486
# dV/dt -= inter.F(t)*R*T/P
24662487
# d/dni(dni/dt) -= inter.F(t)/N
2467-
# d/dV(dni/dt) -= -inter.F(t)*ns/N^2 * dN/dV
2488+
# d/dV(dn/dt) -= -inter.F(t)*ns/N^2 * dN/dV
24682489
# Note: dN/dV = d(PV/RT)/dV = P/RT
24692490
# -= -inter.F(t)*ns/N^2 * P/RT
24702491
# -= -inter.F(t)*ns/N/V
@@ -2488,16 +2509,13 @@ end
24882509
@views @inbounds @fastmath jac[domain.indexes[1]:domain.indexes[2],domain.indexes[3]] .+= -kLAs.*inter.cs/(V*V)
24892510

24902511
# condensation
2491-
# outlet
2492-
# d/dni(dV/dt) = dflow/dni*R*T/P = kLAs[i]/kHs[i]*R*T/P
2493-
# d/dV(dV/dt) = dflow/dV *R*T/P = 0
2494-
# d/dV(dni/dt) = dflow_i/dV = 0
2495-
# dflow/dni = 0
2496-
# dflow/dV = 0
2497-
# flow = sum(kLAs.*ns./kHs)
2498-
# dflowdni = kLAs[i]/kHs[i]
2499-
# dflowdV = 0
2500-
# dflow_i/dV = 0
2512+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2513+
# dV/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/P
2514+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2515+
# d/dV(dn/dt) -= kLAs.*inter.V.*cs.*(-1/V)*R*T./kHs
2516+
# d/dni(dV/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/P
2517+
# -= d/dni(dni/dt)*R*T/P
2518+
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
25012519
@simd for i in domain.indexes[1]:domain.indexes[2]
25022520
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
25032521
end
@@ -2593,10 +2611,13 @@ end
25932611
cond = kLAs.*inter.V.*cs*R*T./kHs
25942612

25952613
# evaporation
2596-
# inlet
2597-
# dTdt = flow*(inter.H - dot(Us,ns)/N)/(N*Cvave)
2598-
# ddnidTdt = flow*(-Us[i]/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave)
2599-
# d/dni (dP/dt) = P/T * d/dni(dT/dt)
2614+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2615+
# dT/dt += sum(kLAs.*inter.V.*inter.cs)*(inter.H - dot(Us,ns)/N)/(N*Cvave)
2616+
# dP/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/V + dT/dt*P/T
2617+
# d/dni(dni/dt) += 0
2618+
# d/dni(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(-Us[i]/N)/(N*Cvave) + dT/dt * (-1/N*Cvave) * d/dni(N*Cvave)
2619+
# += sum(kLAs.*inter.V.*inter.cs)*(-Us[i]/N)/(N*Cvave) + dT/dt * (-1/Cvave) * d/dni(Cvave)
2620+
# d/dni(dP/dt) += d/dni(dT/dt)*P/T
26002621
flow = sum(evap)
26012622
@fastmath dTdt = flow*(inter.H - dot(Us,ns)/N)/(N*Cvave)
26022623
@simd for i in domain.indexes[1]:domain.indexes[2]
@@ -2606,14 +2627,14 @@ end
26062627
@inbounds @fastmath jac[domain.indexes[4],i] += P/T*ddnidTdt
26072628
end
26082629

2609-
26102630
# condensation
2611-
# outlet
2612-
# flow = sum(kLAs.*ns./kHs)
2613-
# dflowdni = kLAs[i]/kHs[i]
2614-
# dTdt = flow*(P*V/N)/(N*Cvave)
2615-
# ddnidTdt = ( dflowdni *P*V/N)/(N*Cvave)-dTdt*(dCvavedni/Cvave) = (kLAs[i]/kHs[i]*P*V/N)/(N*Cvave) -dTdt*(dCvavedni/Cvave)
2616-
# d/dni (dP/dt) = dflowdni *R*T/V + P/T * d/dni(dT/dt) = kLAs[i]/kHs[i]*R*T/V + P/T * d/dni(dT/dt)
2631+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2632+
# dT/dt -= (P*V/N*sum(kLAs.*inter.V.*cs*R*T./kHs))/(N*Cvave)
2633+
# dP/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/V + dT/dt*P/T
2634+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2635+
# d/dni(dT/dt) -= (P*V/N*ddnidnidt)/(N*Cvave) - dTdt*(dCvavedni/Cvave)
2636+
# d/dni(dP/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/V + P/T*ddnidTdt
2637+
# -= d/dni(dni/dt)*R*T/V + P/T*ddnidTdt
26172638
flow = sum(cond)
26182639
@fastmath dTdt = (P*V/N*flow)/(N*Cvave)
26192640
@simd for i in domain.indexes[1]:domain.indexes[2]
@@ -2750,16 +2771,23 @@ end
27502771
evap = kLAs.*inter.V.*inter.cs
27512772

27522773
# evaporation
2753-
# inlet
2754-
# dTdt = flow*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2755-
# ddnidTdt = flow*(-Hs[i]/N)/(N*Cpave)-dTdt*(dCpavedni/Cpave)
2756-
# d/dni (dV/dt) = V/T * d/dni(dT/dt)
2757-
# d/dV (dT/dt) = flow*(dot(Hs, ns)/N)/V/(N*Cpave)
2758-
# d/dV (dV/dt) = dflow/dV*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = flow/V*R*T/P + dT/dt/T + V/T * d/dV(dT/dt) = flow/N + dT/dt/T + V/T * d/dV(dT/dt)
2759-
# d/dV(dni/dt) = dflow_i/dV = kLAs[i]*inter.cs[i]
2760-
# dflowdV = sum(kLAs.*inter.cs) = flow/V
2761-
# flow_i = kLAs[i]*inter.cs[i]*V
2762-
# dflow_i/dV = kLAs[i]*inter.cs[i]
2774+
# dn/dt .+= kLAs.*inter.V.*inter.cs
2775+
# dT/dt += sum(kLAs.*inter.V.*inter.cs)*(inter.H - dot(Hs,ns)/N)/(N*Cpave)
2776+
# dV/dt += sum(kLAs.*inter.V.*inter.cs)*R*T/P + dT/dt*V/T
2777+
# d/dni(dni/dt) += 0
2778+
# d/dV(dni/dt) += 0
2779+
# d/dni(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(-Hs[i]/N)/(N*Cpave) + dTdt * (-1/(N*Cpave)) * d/dni(N*Cpave)
2780+
# += sum(kLAs.*inter.V.*inter.cs)*(-Hs[i]/N)/(N*Cpave) + dTdt * (-1/(Cpave)) * d/dni(Cpave)
2781+
# Note: Cpave = dot(cpdivR,ns)*R/N
2782+
# Note: dCpavedni = cpdivR[i]*R/N
2783+
# d/dV(dT/dt) += sum(kLAs.*inter.V.*inter.cs)*(dot(Hs,ns)/N^2*dN/dV)/(N*Cpave) + dTdt * (-1/(N*Cpave)) * d/dV(N*Cpave)
2784+
# Note: dN/dV = d(PV/RT)/dV = P/RT
2785+
# Note: Cpave = dot(cpdivR,ns)*R/N
2786+
# Note: d/dV(Cpave) = dot(cpdivR,ns)*R*(RT/PV)*(-1/V) = -Cpave/V
2787+
# Note: d(N*Cpave)/dV = (dN/dV*Cpave + N*dCpave/dV) = (P/RT*Cpave + PV/RT(-Cpave/V)) = 0
2788+
# += sum(kLAs.*inter.V.*inter.cs)*(dot(Hs,ns)/N/V)/(N*Cpave)
2789+
# d/dni(dV/dt) += V/T*d/dni(dT/dt)
2790+
# d/dV(dV/dt) += d/dV(dT/dt)*V/T + dT/dt/T
27632791
flow = sum(evap)
27642792
@fastmath H = dot(Hs,ns)/N
27652793
@fastmath dTdt = flow*(inter.H - H)/(N*Cpave)
@@ -2775,15 +2803,16 @@ end
27752803
@inbounds @fastmath jac[domain.indexes[4],domain.indexes[4]] += flow/N + dTdt/T + V/T*ddVdTdt
27762804

27772805
# condensation
2778-
# outlet
2779-
# dTdt = 0
2780-
# d/dni (dV/dt) = dflow/dni *R*T/P = kLAs[i]/kHs[i]*R*T/P
2781-
# d/dV(dV/dt) = dflowdV *R*T/P = 0
2782-
# d/dV(dni/dt) = dflow_i/dV = 0
2783-
# flow_i = kLAs[i]*ns[i]/kHs[i]
2784-
# dflow/dni = kLAs[i]/kHs[i]
2785-
# dflow/dV = 0
2786-
# dflow_i/dV = 0
2806+
# dn/dt .-= kLAs.*inter.V.*cs*R*T./kHs
2807+
# dT/dt -= 0
2808+
# dV/dt -= sum(kLAs.*inter.V.*cs*R*T./kHs)*R*T/P
2809+
# d/dni(dni/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]
2810+
# d/dV(dn/dt) -= kLAs.*inter.V.*cs.*(-1/V)*R*T./kHs
2811+
# d/dni(dT/dt) -= 0
2812+
# d/dV(dT/dt) -= 0
2813+
# d/dni(dV/dt) -= kLAs[i]*inter.V/V*R*T/kHs[i]*R*T/P
2814+
# -= d/dni(dni/dt)*R*T/P
2815+
# d/dV(dV/dt) -= sum(d/dV(dn/dt))*R*T/P
27872816
@simd for i in domain.indexes[1]:domain.indexes[2]
27882817
@inbounds @fastmath jac[i,i] -= kLAs[i]/kHs[i]
27892818
end
@@ -2853,10 +2882,16 @@ end
28532882
elseif isa(inter,kLAkHCondensationEvaporationWithReservoir) && domain == inter.domain
28542883
kLAs = map.(inter.kLAs,T)
28552884

2856-
#evaporation
2885+
# evaporation
2886+
# dn/dt .-= kLAs.*ns
2887+
# d/dni(dni/dt) -= kLAs[i]
28572888
@simd for i in domain.indexes[1]:domain.indexes[2]
28582889
@inbounds @fastmath jac[i,i] -= kLAs[i]
28592890
end
2891+
2892+
# condensation
2893+
# dn/dt .+= kLAs.*inter.molefractions*inter.P/kHs*V
2894+
nothing
28602895
elseif isa(inter,VolumetricFlowRateOutlet) && domain == inter.domain
28612896
# dn/dt .-= inter.Vout(t)*ns/V
28622897
# d/dni(dni/dt) -= inter.Vout(t)/V

0 commit comments

Comments
 (0)