diff --git a/tests/test_homo_wells.py b/tests/test_homo_wells.py index 6120569..8775585 100644 --- a/tests/test_homo_wells.py +++ b/tests/test_homo_wells.py @@ -34,10 +34,10 @@ # GEO N = sys.argv[3] -N = float(N) +N = int(N) Nx = N Ny = N -L = 20 +L = 20. geo = GeoModel(Nx, Ny, params, L, L) Length = geo.Length @@ -49,16 +49,17 @@ prod_points = [[Length/4, Length_y/2]] inj_points = [[3*Length/4, Length_y/2]] case = TestCase(params, geo, well_case = "test0", constant_rate = True)#, prod_points = prod_points, inj_points = inj_points) +#case = TestCase(params, geo, prod_points = [], inj_points = [], constant_rate = True) parameters = { "snes_type": "newtonls", - "snes_monitor": True, - "snes_converged_reason": True, + "snes_monitor": None, + "snes_converged_reason": None, "snes_max_it": 15, "ksp_type": "fgmres", - "ksp_converged_reason": True, - "ksp_view": False, + "ksp_converged_reason": None, + #"ksp_view": None, "ksp_max_it": 200, #"ksp_view": True, } diff --git a/tests_twophase/test_homo_wells.py b/tests_twophase/test_homo_wells.py index f304a1e..d4b3f10 100644 --- a/tests_twophase/test_homo_wells.py +++ b/tests_twophase/test_homo_wells.py @@ -30,16 +30,17 @@ # GEO N = sys.argv[3] -N = float(N) +N = int(N) Nx = N Ny = N -L = 50 +L = 50. geo = GeoModel(Nx, Ny, params, L, L) # CASE case = TestCase(params, geo, well_case = "test0", constant_rate = True)#, prod_points = prod_points, inj_points = inj_points) +case = TestCase(params, geo, well_case = None, prod_points = [], inj_points = [], constant_rate = True) parameters = { "snes_type": "newtonls", diff --git a/thermalporous/mixingcase.py b/thermalporous/mixingcase.py index f4d4518..3b6c252 100644 --- a/thermalporous/mixingcase.py +++ b/thermalporous/mixingcase.py @@ -34,7 +34,10 @@ def init_IC(self, phases = "Single phase"): ic.sub(1).assign(ic_T) elif self.mixing_case == "heavyonlight": if phases == "Single phase": raise SystemExit("Error: heavyonlight case not defined for Single Phase") - ic_S = interpolate(conditional(gt(x[2],self.geo.Length_z/2.), 0.0, 1.0), self.V) + if self.geo.dim == 2: + ic_S = interpolate(conditional(gt(x[1],self.geo.Length_y/2.), 0.0, 1.0), self.V) + elif self.geo.dim == 3: + ic_S = interpolate(conditional(gt(x[2],self.geo.Length_z/2.), 0.0, 1.0), self.V) ic.sub(0).assign(Constant(p_ref)) ic.sub(1).assign(Constant(T_cold)) ic.sub(2).assign(ic_S) diff --git a/thermalporous/preconditioners.py b/thermalporous/preconditioners.py index 7a61972..d9b3598 100644 --- a/thermalporous/preconditioners.py +++ b/thermalporous/preconditioners.py @@ -62,9 +62,16 @@ def initialize(self, pc): if geo.dim == 2: Delta_h = sqrt(jump(x_func)**2 + jump(y_func)**2) + + if geo.gravity2D: + g = params.g + flow = jump(p0)/Delta_h - g*avg(oil_rho(p0,T0))*(abs(n[1]('+'))+abs(n[1]('-')))/2 + else: + flow = jump(p0)/Delta_h + # conservation of energy equation a_Eaccum = phi*c_v*(oil_rho(p0,T0)*T)/dt*r*dx + (1-phi)*rho_r*c_r*(T)/dt*r*dx - a_advec = K_facet*conditional(gt(jump(p0), 0.0), T('+')*oil_rho(p0('+'),T0('+'))/oil_mu(T0('+')), T('-')*oil_rho(p0('-'),T0('-'))/oil_mu(T0('-')))*c_v*jump(r)*jump(p0)/Delta_h*dS + a_advec = K_facet*conditional(gt(flow, 0.0), T('+')*oil_rho(p0('+'),T0('+'))/oil_mu(T0('+')), T('-')*oil_rho(p0('-'),T0('-'))/oil_mu(T0('-')))*c_v*jump(r)*flow*dS a_diff = kT_facet*jump(T)/Delta_h*jump(r)*dS a = a_Eaccum + a_diff + a_advec @@ -217,9 +224,16 @@ def initialize(self, pc): if geo.dim == 2: Delta_h = sqrt(jump(x_func)**2 + jump(y_func)**2) + if geo.gravity2D: + g = params.g + flow_o = jump(p0)/Delta_h - g*avg(oil_rho(p0,T0))*(abs(n[1]('+'))+abs(n[1]('-')))/2 + flow_w = jump(p0)/Delta_h - g*avg(water_rho(p0,T0))*(abs(n[1]('+'))+abs(n[1]('-')))/2 + else: + flow_o = jump(p0)/Delta_h + flow_w = jump(p0)/Delta_h # conservation of energy equation a_Eaccum = phi*c_v_o*S0*(oil_rho(p0,T0)*T)/dt*r*dx + phi*c_v_w*(1.0 - S0)*(water_rho(p0,T0)*T)/dt*r*dx + (1-phi)*rho_r*c_r*(T)/dt*r*dx - a_advec = K_facet*conditional(gt(jump(p0), 0.0), T('+')*rel_perm_w(S0('+'))*water_rho(p0('+'),T0('+'))/water_mu(T0('+')), T('-')*rel_perm_w(S0('-'))*water_rho(p0('-'),T0('-'))/water_mu(T0('-')))*c_v_w*jump(r)*jump(p0)/Delta_h*dS + K_facet*conditional(gt(jump(p0), 0.0), T('+')*rel_perm_o(S0('+'))*oil_rho(p0('+'),T0('+'))/oil_mu(T0('+')), T('-')*rel_perm_o(S0('-'))*oil_rho(p0('-'),T0('-'))/oil_mu(T0('-')))*c_v_o*jump(r)*jump(p0)/Delta_h*dS + a_advec = K_facet*conditional(gt(flow_w, 0.0), T('+')*rel_perm_w(S0('+'))*water_rho(p0('+'),T0('+'))/water_mu(T0('+')), T('-')*rel_perm_w(S0('-'))*water_rho(p0('-'),T0('-'))/water_mu(T0('-')))*c_v_w*jump(r)*flow_w*dS + K_facet*conditional(gt(flow_o, 0.0), T('+')*rel_perm_o(S0('+'))*oil_rho(p0('+'),T0('+'))/oil_mu(T0('+')), T('-')*rel_perm_o(S0('-'))*oil_rho(p0('-'),T0('-'))/oil_mu(T0('-')))*c_v_o*jump(r)*flow_o*dS a_diff = kT_facet*jump(T)/Delta_h*jump(r)*dS a = a_Eaccum + a_diff + a_advec diff --git a/thermalporous/singlephase.py b/thermalporous/singlephase.py index cb0c17d..3028e1c 100644 --- a/thermalporous/singlephase.py +++ b/thermalporous/singlephase.py @@ -4,7 +4,7 @@ from firedrake.utils import cached_property class SinglePhase(ThermalModel): - def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_save = 2, small_dt_start = True, checkpointing = {}, solver_parameters = None, filename = "results/results.txt", dt_init_fact = 2**(-10), vector = False): + def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_save = 2, small_dt_start = True, checkpointing = {}, solver_parameters = None, filename = "results/results.txt", dt_init_fact = 2**(-10), vector = False, gravity2D = False): self.name = "Single phase" self.geo = geo self.case = case @@ -24,6 +24,7 @@ def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_ self.solver_parameters = solver_parameters self.init_solver_parameters() self.scaled_eqns = False + self.geo.gravity2D = gravity2D if self.geo.dim == 2: self.init_variational_form = self.init_variational_form_2D elif self.geo.dim == 3: @@ -96,7 +97,13 @@ def init_variational_form_2D(self): K_facet = (K_x_facet*(abs(n[0]('+'))+abs(n[0]('-')))/2 + K_y_facet*(abs(n[1]('+'))+abs(n[1]('-')))/2) #need form to be symmetric wrt to '+' and '-' - kT_facet = conditional(gt(avg(kT), 0.0), kT('+')*kT('-') / avg(kT), 0.0) + kT_facet = conditional(gt(avg(kT), 0.0), kT('+')*kT('-') / avg(kT), 0.0) + + if self.geo.gravity2D: + g = self.params.g + flow = jump(p)/Delta_h - g*avg(rho_o)*(abs(n[1]('+'))+abs(n[1]('-')))/2 + else: + flow = jump(p)/Delta_h # Weight for mass equation if self.scaled_eqns: @@ -108,10 +115,10 @@ def init_variational_form_2D(self): ## Solve a coupled problem # conservation of mass equation a_accum = phi*(rho_o - oil_rho(p_,T_))/self.dt*q*dx - a_flow = K_facet*conditional(gt(jump(p), 0.0), rho_o('+')/mu_o('+'), rho_o('-')/mu_o('-'))*jump(q)*jump(p)/Delta_h*dS + a_flow = K_facet*conditional(gt(flow, 0.0), rho_o('+')/mu_o('+'), rho_o('-')/mu_o('-'))*jump(q)*flow*dS # conservation of energy equation a_Eaccum = phi*c_v*(rho_o*T - oil_rho(p_,T_)*T_)/self.dt*r*dx + (1-phi)*rho_r*c_r*(T - T_)/self.dt*r*dx - a_advec = K_facet*conditional(gt(jump(p), 0.0), T('+')*rho_o('+')/mu_o('+'), T('-')*rho_o('-')/mu_o('-'))*c_v*jump(r)*jump(p)/Delta_h*dS + a_advec = K_facet*conditional(gt(flow, 0.0), T('+')*rho_o('+')/mu_o('+'), T('-')*rho_o('-')/mu_o('-'))*c_v*jump(r)*flow*dS a_diff = kT_facet*jump(T)/Delta_h*jump(r)*dS a = m_w*a_accum + m_w*a_flow + a_Eaccum + a_diff + a_advec diff --git a/thermalporous/twophase.py b/thermalporous/twophase.py index a941985..2432710 100644 --- a/thermalporous/twophase.py +++ b/thermalporous/twophase.py @@ -5,7 +5,7 @@ from firedrake.utils import cached_property class TwoPhase(ThermalModel): - def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_save = 2, small_dt_start = True, checkpointing = {}, solver_parameters = None, filename = "results/results.txt", dt_init_fact = 2**(-10), vector = False): + def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_save = 2, small_dt_start = True, checkpointing = {}, solver_parameters = None, filename = "results/results.txt", dt_init_fact = 2**(-10), vector = False, gravity2D = False): self.name = "Two-phase" self.geo = geo self.case = case @@ -28,6 +28,7 @@ def __init__(self, geo, case, params, end = 1.0, maxdt = 0.005, save = False, n_ self.small_dt_start = small_dt_start self.scaled_eqns = True # Weights equations such that they are of similar scale self.pressure_eqn = True # Water equation -> pressure equation. Needed for Schur complement approach. + self.geo.gravity2D = gravity2D if self.geo.dim == 2: self.init_variational_form = self.init_variational_form_2D elif self.geo.dim == 3: @@ -97,7 +98,12 @@ def init_variational_form_2D(self): (p_, T_, S_o_) = split(self.u_) q, r, s = TestFunctions(W) - + rho_o = oil_rho(p, T) + rho_w = water_rho(p, T) + mu_o = oil_mu(T) + mu_w = water_mu(T) + kr_o = rel_perm_o(S_o) + kr_w = rel_perm_w(S_o) # Define facet quantities n = FacetNormal(mesh) @@ -129,29 +135,37 @@ def init_variational_form_2D(self): p_w = 1.0 o_w = 1.0 + if self.geo.gravity2D: + g = self.params.g + flow_o = jump(p)/Delta_h - g*avg(rho_o)*(abs(n[1]('+'))+abs(n[1]('-')))/2 + flow_w = jump(p)/Delta_h - g*avg(rho_w)*(abs(n[1]('+'))+abs(n[1]('-')))/2 + else: + flow_o = jump(p)/Delta_h + flow_w = jump(p)/Delta_h + ## Solve a coupled problem # conservation of mass equation WATER - "pressure equation" - a_accum_w =phi*(water_rho(p,T)*(1.0 - S_o) - water_rho(p_,T_)*(1.0 - S_o_))/self.dt*q*dx - a_flow_w = K_facet*conditional(gt(jump(p), 0.0), rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*jump(q)*jump(p)/Delta_h*dS + a_accum_w =phi*(rho_w*(1.0 - S_o) - water_rho(p_,T_)*(1.0 - S_o_))/self.dt*q*dx + a_flow_w = K_facet*conditional(gt(flow_w, 0.0), kr_w('+')*rho_w('+')/mu_w('+'), kr_w('-')*rho_w('-')/mu_w('-'))*jump(q)*flow_w*dS # conservation of mass equation OIL - "saturation equation" - a_accum_o = phi*(oil_rho(p,T)*S_o - oil_rho(p_,T_)*S_o_)/self.dt*s*dx - a_flow_o = K_facet*conditional(gt(jump(p), 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(s)*jump(p)/Delta_h*dS + a_accum_o = phi*(rho_o*S_o - oil_rho(p_,T_)*S_o_)/self.dt*s*dx + a_flow_o = K_facet*conditional(gt(flow_o, 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(s)*flow_o*dS ## WEIGHTED SUM - pressure equation if self.pressure_eqn: - a_accum_w = c_v_w*a_accum_w + c_v_o*(phi*(oil_rho(p,T)*S_o - oil_rho(p_,T_)*S_o_)/self.dt*q*dx) - a_flow_w = c_v_w*a_flow_w + c_v_o*(K_facet*conditional(gt(jump(p), 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(q)*jump(p)/Delta_h*dS) + a_accum_w = c_v_w*a_accum_w + c_v_o*(phi*(rho_o*S_o - oil_rho(p_,T_)*S_o_)/self.dt*q*dx) + a_flow_w = c_v_w*a_flow_w + c_v_o*(K_facet*conditional(gt(flow_o, 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(q)*flow_o*dS) # conservation of energy equation - a_Eaccum = phi*c_v_w*(water_rho(p,T)*(1.0 - S_o)*T - water_rho(p_,T_)*(1.0 - S_o_)*T_)/self.dt*r*dx + phi*c_v_o*(oil_rho(p,T)*S_o*T - oil_rho(p_,T_)*S_o_*T_)/self.dt*r*dx + (1-phi)*rho_r*c_r*(T - T_)/self.dt*r*dx - a_advec = K_facet*conditional(gt(jump(p), 0.0), T('+')*rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), T('-')*rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*c_v_w*jump(r)*jump(p)/Delta_h*dS + K_facet*conditional(gt(jump(p), 0.0), T('+')*rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), T('-')*rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*c_v_o*jump(r)*jump(p)/Delta_h*dS + a_Eaccum = phi*c_v_w*(rho_w*(1.0 - S_o)*T - water_rho(p_,T_)*(1.0 - S_o_)*T_)/self.dt*r*dx + phi*c_v_o*(rho_o*S_o*T - oil_rho(p_,T_)*S_o_*T_)/self.dt*r*dx + (1-phi)*rho_r*c_r*(T - T_)/self.dt*r*dx + a_advec = K_facet*conditional(gt(flow_w, 0.0), T('+')*kr_w('+')*rho_w('+')/mu_w('+'), T('-')*kr_w('-')*rho_w('-')/mu_w('-'))*c_v_w*jump(r)*flow_w*dS + K_facet*conditional(gt(flow_o, 0.0), T('+')*kr_o('+')*rho_o('+')/mu_o('+'), T('-')*kr_o('-')*rho_o('-')/mu_o('-'))*c_v_o*jump(r)*flow_o*dS a_diff = kT_facet*jump(T)/Delta_h*jump(r)*dS a = p_w*a_accum_w + p_w*a_flow_w + o_w*a_accum_o + o_w*a_flow_o + a_Eaccum + a_diff + a_advec self.F = a - rhow_o = oil_rho(p, T) - rhow_w = water_rho(p, T) + rhow_o = rho_o + rhow_w =rho_w rhow = water_rho(p, T_inj) ## Source terms using global deltas @@ -246,8 +260,12 @@ def init_variational_form_3D(self): (p_, T_, S_o_) = split(self.u_) q, r, s = TestFunctions(W) - - + rho_o = oil_rho(p, T) + rho_w = water_rho(p, T) + mu_o = oil_mu(T) + mu_w = water_mu(T) + kr_o = rel_perm_o(S_o) + kr_w = rel_perm_w(S_o) # Define facet quantities n = FacetNormal(mesh) @@ -270,8 +288,8 @@ def init_variational_form_3D(self): kT_facet = conditional(gt(avg(kT), 0.0), kT('+')*kT('-') / avg(kT), 0.0) - z_flow_w = jump(p)/Delta_h - g*avg(water_rho(p,T)) - z_flow_o = jump(p)/Delta_h - g*avg(oil_rho(p,T)) + z_flow_w = jump(p)/Delta_h - g*avg(rho_w) + z_flow_o = jump(p)/Delta_h - g*avg(rho_o) # weights for pressure and oil equations. Only uses initial saturation if self.scaled_eqns: @@ -286,32 +304,32 @@ def init_variational_form_3D(self): ## Solve a coupled problem # conservation of mass equation WATER - "pressure equation" - a_accum_w = phi*(water_rho(p,T)*(1.0 - S_o) - water_rho(p_,T_)*(1.0 - S_o_))/self.dt*q*dx - a_flow_w = K_facet*conditional(gt(jump(p), 0.0), rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*jump(q)*jump(p)/Delta_h*dS_v - a_flow_w_z = K_z_facet*conditional(gt(z_flow_w, 0.0), rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*jump(q)*z_flow_w*dS_h + a_accum_w = phi*(rho_w*(1.0 - S_o) - water_rho(p_,T_)*(1.0 - S_o_))/self.dt*q*dx + a_flow_w = K_facet*conditional(gt(jump(p), 0.0), kr_w('+')*rho_w('+')/mu_w('+'), kr_w('-')*rho_w('-')/mu_w('-'))*jump(q)*jump(p)/Delta_h*dS_v + a_flow_w_z = K_z_facet*conditional(gt(z_flow_w, 0.0), kr_w('+')*rho_w('+')/mu_w('+'), kr_w('-')*rho_w('-')/mu_w('-'))*jump(q)*z_flow_w*dS_h # conservation of mass equation OIL - "saturation equation" - a_accum_o = phi*(oil_rho(p,T)*S_o - oil_rho(p_,T_)*S_o_)/self.dt*s*dx - a_flow_o = K_facet*conditional(gt(jump(p), 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(s)*jump(p)/Delta_h*dS_v - a_flow_o_z = K_z_facet*conditional(gt(z_flow_o, 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(s)*z_flow_o*dS_h + a_accum_o = phi*(rho_o*S_o - oil_rho(p_,T_)*S_o_)/self.dt*s*dx + a_flow_o = K_facet*conditional(gt(jump(p), 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(s)*jump(p)/Delta_h*dS_v + a_flow_o_z = K_z_facet*conditional(gt(z_flow_o, 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(s)*z_flow_o*dS_h ## WEIGHTED SUM - pressure equation if self.pressure_eqn: - a_accum_w = c_v_w*a_accum_w + c_v_o*(phi*(oil_rho(p,T)*S_o - oil_rho(p_,T_)*S_o_)/self.dt*q*dx) - a_flow_w = c_v_w*a_flow_w + c_v_o*(K_facet*conditional(gt(jump(p), 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(q)*jump(p)/Delta_h*dS_v) - a_flow_w_z = c_v_w*a_flow_w_z + c_v_o*K_z_facet*conditional(gt(z_flow_o, 0.0), rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*jump(q)*z_flow_o*dS_h + a_accum_w = c_v_w*a_accum_w + c_v_o*(phi*(rho_o*S_o - oil_rho(p_,T_)*S_o_)/self.dt*q*dx) + a_flow_w = c_v_w*a_flow_w + c_v_o*(K_facet*conditional(gt(jump(p), 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(q)*jump(p)/Delta_h*dS_v) + a_flow_w_z = c_v_w*a_flow_w_z + c_v_o*K_z_facet*conditional(gt(z_flow_o, 0.0), kr_o('+')*rho_o('+')/mu_o('+'), kr_o('-')*rho_o('-')/mu_o('-'))*jump(q)*z_flow_o*dS_h # conservation of energy equation - a_Eaccum = phi*c_v_w*(water_rho(p,T)*(1.0 - S_o)*T - water_rho(p_,T_)*(1.0 - S_o_)*T_)/self.dt*r*dx + phi*c_v_o*(oil_rho(p,T)*S_o*T - oil_rho(p_,T_)*S_o_*T_)/self.dt*r*dx + (1-phi)*rho_r*c_r*(T - T_)/self.dt*r*dx - a_advec = K_facet*conditional(gt(jump(p), 0.0), T('+')*rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), T('-')*rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*c_v_w*jump(r)*jump(p)/Delta_h*dS_v + K_facet*conditional(gt(jump(p), 0.0), T('+')*rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), T('-')*rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*c_v_o*jump(r)*jump(p)/Delta_h*dS_v - a_advec_z = K_z_facet*conditional(gt(z_flow_w, 0.0), T('+')*rel_perm_w(S_o('+'))*water_rho(p('+'),T('+'))/water_mu(T('+')), T('-')*rel_perm_w(S_o('-'))*water_rho(p('-'),T('-'))/water_mu(T('-')))*c_v_w*jump(r)*z_flow_w*dS_h + K_z_facet*conditional(gt(z_flow_o, 0.0), T('+')*rel_perm_o(S_o('+'))*oil_rho(p('+'),T('+'))/oil_mu(T('+')), T('-')*rel_perm_o(S_o('-'))*oil_rho(p('-'),T('-'))/oil_mu(T('-')))*c_v_o*jump(r)*z_flow_o*dS_h + a_Eaccum = phi*c_v_w*(rho_w*(1.0 - S_o)*T - water_rho(p_,T_)*(1.0 - S_o_)*T_)/self.dt*r*dx + phi*c_v_o*(rho_o*S_o*T - oil_rho(p_,T_)*S_o_*T_)/self.dt*r*dx + (1-phi)*rho_r*c_r*(T - T_)/self.dt*r*dx + a_advec = K_facet*conditional(gt(jump(p), 0.0), T('+')*kr_w('+')*rho_w('+')/mu_w('+'), T('-')*kr_w('-')*rho_w('-')/mu_w('-'))*c_v_w*jump(r)*jump(p)/Delta_h*dS_v + K_facet*conditional(gt(jump(p), 0.0), T('+')*kr_o('+')*rho_o('+')/mu_o('+'), T('-')*kr_o('-')*rho_o('-')/mu_o('-'))*c_v_o*jump(r)*jump(p)/Delta_h*dS_v + a_advec_z = K_z_facet*conditional(gt(z_flow_w, 0.0), T('+')*kr_w('+')*rho_w('+')/mu_w('+'), T('-')*kr_w('-')*rho_w('-')/mu_w('-'))*c_v_w*jump(r)*z_flow_w*dS_h + K_z_facet*conditional(gt(z_flow_o, 0.0), T('+')*kr_o('+')*rho_o('+')/mu_o('+'), T('-')*kr_o('-')*rho_o('-')/mu_o('-'))*c_v_o*jump(r)*z_flow_o*dS_h a_diff = kT_facet*jump(T)/Delta_h*jump(r)*(dS_v + dS_h) a = p_w*a_accum_w + p_w*a_flow_w + p_w*a_flow_w_z + o_w*a_accum_o + o_w*a_flow_o +o_w* a_flow_o_z + a_Eaccum + a_advec + a_advec_z + a_diff self.F = a - rhow_o = oil_rho(p, T) - rhow_w = water_rho(p, T) + rhow_o = rho_o + rhow_w = rho_w rhow = water_rho(p, T_inj) ## Source terms using global deltas