Merge pull request #247 from brown-ccv/246-updates-on-the-odechamber-…

…function

246 updates on the odechamber function

.gitignore

@@ -1,4 +1,5 @@
 /.vscode
 202*
 sysimage_chamber*.so
-docs/build/
+docs/build/
+src/pplot.jl

src/runcode-func.jl

@@ -36,21 +36,12 @@ function odeChamber(
     L_e, L_m = param.L_e, param.L_m
     mm_h2o, mm_co2 = param.mm_h2o, param.mm_co2
     T_in = param.T_in
-    P_lit_0, dP_lit_dt, dP_lit_dt_0, P_lit_drop_max = param.P_lit_0,
-    param.dP_lit_dt, param.dP_lit_dt_0,
-    param.P_lit_drop_max
 
-    P0plusDP = u[1]
+    P = u[1]
     T = u[2]
     eps_g = u[3]
     X_co2 = u[7]
-    P_lit = P_lit_0 + dP_lit_dt_0 * t
-    if P_lit < P_lit_0 - P_lit_drop_max
-        P_lit = P_lit_0 - P_lit_drop_max
-        dP_lit_dt = 0.0
-        param.dP_lit_dt = 0.0
-    end
-    P = P_lit + P0plusDP - P_lit_0
+
     # effective gas molar mass
     m_g = mm_co2 * X_co2 + mm_h2o * (1 - X_co2)
 
@@ -94,8 +85,10 @@ function odeChamber(
 
     if phase == 3
         C_co2 = C_co2_t
+        C_h2o = m_eq
     elseif phase == 2
-        C_co2 = m_co2
+        C_co2 = M_co2 / (rho_m * eps_m * V)
+        C_h2o = M_h2o / (rho_m * eps_m * V)
     end
 
     # specific heat of gas
@@ -147,7 +140,7 @@ function odeChamber(
         dm_eq_dP,
         rho_m,
         eps_m,
-        m_eq,
+        C_h2o,
         drho_m_dP,
         drc_dT,
         drho_x_dT,
@@ -161,7 +154,6 @@ function odeChamber(
         m_h2o,
         m_co2,
         deps_x_dmco2_t,
-        dP_lit_dt,
         Hdot_in,
         Hdot_out,
         c_x,
@@ -189,14 +181,13 @@ function odeChamber(
 
     # coefficients in the system of unknowns Ax = B, here x= [dP/dt dT/dt deps_g/dt dX_co2/dt]
     if phase == 3
-        dDP_dt, dT_dt, deps_g_dt, dX_co2_dt = A \ b
+        dP_dt, dT_dt, deps_g_dt, dX_co2_dt = A \ b
     elseif phase == 2
-        dDP_dt, dT_dt = A \ b
+        dP_dt, dT_dt = A \ b
         deps_g_dt, dX_co2_dt = 0.0, 0.0
     end
-    dP_dt = dDP_dt + dP_lit_dt
 
-    du[1] = dDP_dt
+    du[1] = dP_dt
     du[2] = dT_dt
     du[3] = deps_g_dt
     du[4] = dV_dP * dP_dt + dV_dT * dT_dt + V * P_loss
@@ -285,7 +276,8 @@ function stopChamber_MT(
     out[5] = eps_x - 0.5
     out[6] = m_h2o_melt - m_eq_max
     out[7] = -(P0plusDP - param.P_lit_0 + DP_crit)
-    return out[8] = m_co2_melt - C_co2_sat
+    out[8] = m_co2_melt - C_co2_sat
+    return out
 end
 
 """
@@ -312,7 +304,6 @@ function affect!(
     param_IC_Finder::ParamICFinder{Float64},
     erupt_saved::EruptSaved{Float64},
 )
-    @info("*event idx: ", idx)
     composition = param.composition
     storeTime = param_saved_var.storeTime
     storeTemp = param_saved_var.storeTemp
@@ -340,17 +331,20 @@ function affect!(
         sw.eruption = 1
         rho_g0 = eos_g_rho_g(P_0, int.u[2])
         record_erupt_start(int.t, int.u[3], eps_x0, int.u[5], int.u[6], rho_g0, erupt_saved)
-        @info("time: $(int.t), Reached critical pressure and need to start an eruption")
+        @info(
+            "*event idx: $idx\n  time: $(int.t), Reached critical pressure and need to start an eruption"
+        )
     elseif idx == 4
         sw.eruption = 0
         record_erupt_end(int.t, erupt_saved, param)
-        @info("time: $(int.t), Finished an eruption...")
+        @info("*event idx: $idx\n  time: $(int.t), Finished an eruption...")
     elseif idx == 6 || idx == 8
         phase_here = param_saved_var.phase
         @info(
-            "time: $(int.t), starting ic finder for conversion of phase, phase_here: $phase_here",
+            "*event idx: $idx\n  time: $(int.t), starting ic finder for conversion of phase, phase_here: $phase_here",
         )
-        eps_g_temp, X_co2_temp, C_co2_temp, phase = IC_Finder(
+        eps_g_temp, X_co2_temp, C_co2_temp, phase = ic_phase_conversion(
+            phase_here,
             composition,
             int.u[9],
             int.u[10],
@@ -361,65 +355,10 @@ function affect!(
             int.u[5],
             param_IC_Finder,
         )
-
         param_saved_var.phase = phase
-        if phase_here != phase
-            @info("1st try in IC Finder successful")
-            int.u[3] = eps_g_temp
-            int.u[7] = X_co2_temp
-            C_co2 = C_co2_temp
-        else
-            @info("1st try in IC Finder unsuccessful, trying new IC parameters...")
-            param_IC_Finder.max_count = 150
-            eps_g_temp, X_co2_temp, C_co2_temp, phase = IC_Finder(
-                composition,
-                int.u[9],
-                int.u[10],
-                int.u[8],
-                P_0,
-                int.u[2],
-                int.u[4],
-                int.u[5],
-                param_IC_Finder,
-            )
-            param_saved_var.phase = phase
-            ## change back to initial max_count
-            param_IC_Finder.max_count = 100
-            if phase_here != phase
-                @info("2nd try in IC Finder successful")
-                int.u[3] = eps_g_temp
-                int.u[7] = X_co2_temp
-                C_co2 = C_co2_temp
-            else
-                @info("2nd try in IC Finder unsuccessful, trying new IC parameters...")
-                param_IC_Finder.max_count = 100
-                param_IC_Finder.Tol = param_IC_Finder.Tol * 0.1
-                eps_g_temp, X_co2_temp, C_co2_temp, phase = IC_Finder(
-                    composition,
-                    int.u[9],
-                    int.u[10],
-                    int.u[8],
-                    P_0,
-                    int.u[2],
-                    int.u[4],
-                    int.u[5],
-                    param_IC_Finder,
-                )
-                param_saved_var.phase = phase
-                ## change back to initial Tol
-                param_IC_Finder.Tol = param_IC_Finder.Tol * 10
-                if phase_here != phase
-                    @info("3rd try in IC Finder successful")
-                    int.u[3] = eps_g_temp
-                    int.u[7] = X_co2_temp
-                    C_co2 = C_co2_temp
-                else
-                    @warn("3rd try in IC Finder unsuccessful")
-                end
-            end
-        end
-        @info("phase_here: $phase_here, new_phase: $phase")
-
+        int.u[3] = eps_g_temp
+        int.u[7] = X_co2_temp
+        C_co2 = C_co2_temp
     elseif idx == 1 || idx == 2 || idx == 5 || idx == 7 || idx === nothing
         if idx == 1
             @info("eps_x became 0.")