Skip to content

Commit

Permalink
2D vertical models. Simplification of two-phase notation
Browse files Browse the repository at this point in the history
  • Loading branch information
@tlroy
tlroy committed Mar 5, 2020
1 parent 2120561 commit f3dce72
Show file tree
Hide file tree
Showing 6 changed files with 89 additions and 45 deletions.
13 changes: 7 additions & 6 deletions tests/test_homo_wells.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
}
Expand Down
5 changes: 3 additions & 2 deletions tests_twophase/test_homo_wells.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
5 changes: 4 additions & 1 deletion thermalporous/mixingcase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
18 changes: 16 additions & 2 deletions thermalporous/preconditioners.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
15 changes: 11 additions & 4 deletions thermalporous/singlephase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down
78 changes: 48 additions & 30 deletions thermalporous/twophase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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:
Expand All @@ -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
Expand Down

0 comments on commit f3dce72

Please sign in to comment.