-
Notifications
You must be signed in to change notification settings - Fork 0
/
twp_navier_stokes_p.py
112 lines (97 loc) · 4.15 KB
/
twp_navier_stokes_p.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
from proteus.default_p import *
from proteus.mprans import RANS2P
import numpy as np
from proteus import Context
ct = Context.get()
domain = ct.domain
nd = domain.nd
mesh = domain.MeshOptions
genMesh = mesh.genMesh
movingDomain = ct.movingDomain
T = ct.T # might not be necessary
LevelModelType = RANS2P.LevelModel
if ct.useOnlyVF:
LS_model = None
else:
LS_model = 2
if ct.useRANS >= 1:
Closure_0_model = 5
Closure_1_model = 6
if ct.useOnlyVF:
Closure_0_model = 2
Closure_1_model = 3
if ct.movingDomain:
Closure_0_model += 1
Closure_1_model += 1
else:
Closure_0_model = None
Closure_1_model = None
# for absorption zones (defined as regions)
if hasattr(domain, 'porosityTypes'):
porosityTypes = domain.porosityTypes
dragAlphaTypes = domain.dragAlphaTypes
dragBetaTypes = domain.dragBetaTypes
epsFact_porous = domain.epsFact_porous
else:
porosityTypes = None
dragAlphaTypes = None
dragBetaTypes = None
epsFact_porous = None
coefficients = RANS2P.Coefficients(epsFact=ct.epsFact_viscosity,
sigma=0.0,
rho_0=ct.rho_0,
nu_0=ct.nu_0,
rho_1=ct.rho_1,
nu_1=ct.nu_1,
g=ct.g,
nd=nd,
ME_model=int(ct.movingDomain)+0,
VF_model=int(ct.movingDomain)+1,
LS_model=int(ct.movingDomain)+LS_model,
Closure_0_model=Closure_0_model,
Closure_1_model=Closure_1_model,
epsFact_density=ct.epsFact_density,
stokes=False,
useVF=ct.useVF,
useRBLES=ct.useRBLES,
useMetrics=ct.useMetrics,
eb_adjoint_sigma=1.0,
eb_penalty_constant=ct.weak_bc_penalty_constant,
forceStrongDirichlet=ct.ns_forceStrongDirichlet,
turbulenceClosureModel=ct.ns_closure,
movingDomain=ct.movingDomain,
porosityTypes=porosityTypes,
dragAlphaTypes=dragAlphaTypes,
dragBetaTypes=dragBetaTypes,
epsFact_porous=epsFact_porous,
barycenters=ct.domain.barycenters)
dirichletConditions = {0: lambda x, flag: domain.bc[flag].p_dirichlet.init_cython(),
1: lambda x, flag: domain.bc[flag].u_dirichlet.init_cython(),
2: lambda x, flag: domain.bc[flag].v_dirichlet.init_cython()}
advectiveFluxBoundaryConditions = {0: lambda x, flag: domain.bc[flag].p_advective.init_cython(),
1: lambda x, flag: domain.bc[flag].u_advective.init_cython(),
2: lambda x, flag: domain.bc[flag].v_advective.init_cython()}
diffusiveFluxBoundaryConditions = {0: {},
1: {1: lambda x, flag: domain.bc[flag].u_diffusive.init_cython()},
2: {2: lambda x, flag: domain.bc[flag].v_diffusive.init_cython()}}
if nd == 3:
dirichletConditions[3] = lambda x, flag: domain.bc[flag].w_dirichlet.init_cython()
advectiveFluxBoundaryConditions[3] = lambda x, flag: domain.bc[flag].w_advective.init_cython()
diffusiveFluxBoundaryConditions[3] = {3: lambda x, flag: domain.bc[flag].w_diffusive.init_cython()}
class P_IC:
def uOfXT(self, x, t):
return ct.twpflowPressure_init(x, t)
class U_IC:
def uOfXT(self, x, t):
return 0.0
class V_IC:
def uOfXT(self, x, t):
return 0.0
class W_IC:
def uOfXT(self, x, t):
return 0.0
initialConditions = {0: P_IC(),
1: U_IC(),
2: V_IC()}
if nd == 3:
initialConditions[3] = W_IC()