Add mass to rop calc for film domain

ReactionMechanismGenerator · Aug 8, 2023 · 1adebc2 · 1adebc2
1 parent 12072da
commit 1adebc2
Showing 1 changed file with 10 additions and 4 deletions.
diff --git a/src/Simulation.jl b/src/Simulation.jl
@@ -286,7 +286,7 @@ function rops(ssys::SystemSimulation, t)
     for inter in ssys.interfaces
         if inter isa FragmentBasedReactiveFilmGrowthInterfaceConstantT
             kfs, krevs = getkfskrevs(inter)
-            rops!(ropmat, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], start)
+            rops!(ropmat, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], inter.domainfilm.indexes[3], start)
             start += length(kfs)
         elseif hasproperty(inter, :reactions)
             kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
@@ -423,21 +423,27 @@ function rops!(ropvec, rarray::Array{Int64,2}, cs, kfs, krevs, V, start, ind)
 end
 
 
-function rops!(ropmat, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int64,2}, cs, kfs, krevs, V, start)
+function rops!(ropmat, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int64,2}, cs, kfs, krevs, V, Mws, fragmentindexes, massindex, start)
     numfragmentbasedreacprod, numrxns = size(fragmentbasedrxnarray)
     half = Int(numfragmentbasedreacprod / 2)
     for i = 1:length(kfs)
         R = getrate(rarray, cs, kfs, krevs, V, i)
 
         for j = 1:half
             if fragmentbasedrxnarray[j, i] != 0
-                @fastmath ropmat[fragmentbasedrxnarray[j, i]] -= R
+                @fastmath ropmat[i+start, fragmentbasedrxnarray[j, i]] -= R
+                if !(fragmentbasedrxnarray[j, i] in fragmentindexes)
+                    ropmat[i+start, massindex] += R * Mws[fragmentbasedrxnarray[j, i]]
+                end
             end
         end
 
         for j = half+1:numfragmentbasedreacprod
             if fragmentbasedrxnarray[j, i] != 0
-                @fastmath ropmat[fragmentbasedrxnarray[j, i]] += R
+                @fastmath ropmat[i+start, fragmentbasedrxnarray[j, i]] += R
+                if !(fragmentbasedrxnarray[j, i] in fragmentindexes)
+                    ropmat[i+start, massindex] -= R * Mws[fragmentbasedrxnarray[j, i]]
+                end
             end
         end
     end