Skip to content

Improve Temperature Stabilization #100

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
eyoung55 wants to merge 4 commits into
base: dev_wrap
Choose a base branch
from
Open

Improve Temperature Stabilization #100

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
382b56a
add new options for temperature stabilization which include dt depend…
eyoung55 Jul 7, 2025
2941438
fix bug of hard coding tag 100
eyoung55 Jul 7, 2025
9f89b15
Merge branch 'dev' into temp-stab
eyoung55 Sep 10, 2025
70dc0a1
Merge branch 'dev_wrap' into temp-stab
eyoung55 Jan 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
240 changes: 224 additions & 16 deletions pvade/fluid/FlowManager.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -230,6 +230,7 @@ def build_forms(self, domain, params):
self.s = ufl.TestFunction(self.S)
self.theta_k = dolfinx.fem.Function(self.S, name="temperature ")
self.theta_k1 = dolfinx.fem.Function(self.S)
self.theta_k2 = dolfinx.fem.Function(self.S)
self.theta_r = dolfinx.fem.Function(self.S) # reference temperature

self.mesh_vel = dolfinx.fem.Function(self.V, name="mesh_displacement")
Expand Down Expand Up @@ -267,6 +268,7 @@ def build_forms(self, domain, params):
if self.thermal_analysis:
# initialize air temperature everywhere to be the ambient temperature
self.theta_k1.x.array[:] = PETSc.ScalarType(params.fluid.T_ambient)
self.theta_k2.x.array[:] = PETSc.ScalarType(params.fluid.T_ambient)
self.theta_k.x.array[:] = PETSc.ScalarType(params.fluid.T_ambient)
self.theta_r.x.array[:] = PETSc.ScalarType(params.fluid.T_ambient)

Expand Down Expand Up @@ -401,14 +403,6 @@ def sigma(u, p, nu, rho):
* ufl.dx
)
if self.thermal_analysis:
if self.stabilizing:
# Residual: the "strong" form of the governing equation
self.r = (
(1.0 / self.dt_c) * (self.theta - self.theta_k1)
+ ufl.dot(self.u_k, ufl.nabla_grad(self.theta))
- self.alpha_c * ufl.div(ufl.grad(self.theta))
)

# Define variational problem for step 4: temperature
self.F4 = (
(1.0 / self.dt_c)
Expand All @@ -422,12 +416,224 @@ def sigma(u, p, nu, rho):
)

if self.stabilizing:
# Donea and Huerta 2003 (Eq 2.64)
h = ufl.CellDiameter(domain.fluid.msh)
u_mag = ufl.sqrt(ufl.dot(self.u_k, self.u_k)) # magnitude of vector
Pe = u_mag * h / (2.0 * self.alpha_c) # Peclet number
tau = (h / (2.0 * u_mag)) * (1.0 + 1.0 / Pe) ** (-1)
stab = tau * ufl.dot(self.u_k, ufl.grad(self.s)) * self.r * ufl.dx
# Residual: the "strong" form of the governing equation
self.r = (
(1.0 / self.dt_c) * (self.theta - self.theta_k1)
+ ufl.dot(self.u_k, ufl.nabla_grad(self.theta))
- self.alpha_c * ufl.div(ufl.grad(self.theta))
)

def ufl_mag(input_var, min_val=1e-6):

# sqrt(x_1^2 + x_2^2 + ...)
input_var_mag = ufl.sqrt(ufl.dot(input_var, input_var))

# clip to a minimum value if specified
if min_val > 0.0:
input_var_mag = ufl.max_value(input_var_mag, min_val)

return input_var_mag

STABILIZING_METHOD = 2

if STABILIZING_METHOD == 1:
# Donea and Huerta "Finite Element Methods for Flow Problems" 2003 (starting on Eq 2.64)
h = ufl.CellDiameter(domain.fluid.msh)

u_mag = ufl_mag(self.u_k)

Pe = u_mag * h / (2.0 * self.alpha_c) # Peclet number

# h/(2a) * (coth(Pe) - 1/Pe)
tau_1 = h / (2.0 * u_mag) * (1.0 / ufl.tanh(Pe) - 1.0 / Pe)

# h/(2a) * (1 + 1/Pe) ^ (-1)
tau_2 = (
h / (2.0 * u_mag) * (1.0 + 1.0 / Pe) ** (-1)
) # from codina's algebraic analysis

# h/(2a) * (1 + 9/Pe^2) ^ (-1/2)
tau_3 = (
h / (2.0 * u_mag) * (1.0 + 9.0 / Pe**2) ** (-0.5)
) # shakib et al. fourth order accurate

stab = ufl.dot(tau_1 * self.u_k, ufl.grad(self.s)) * self.r * ufl.dx

elif STABILIZING_METHOD == 2:
# "Finite element stabilization parameters computed from element matrices and vectors"
# Tayfun E. Tezduyar, Yasuo Osawa, 1999, starting Eq 55

h = ufl.CellDiameter(domain.fluid.msh)

u_mag = ufl_mag(self.u_k)

tau_s1 = h / (2.0 * u_mag)
tau_s2 = self.dt_c / 2.0
tau_s3 = h**2 / (4.0 * self.alpha_c)

tau = (
1.0 / tau_s1**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

fac = ufl.sqrt(1.0 / self.dt_c)
# fac = 1.0

stab = (
ufl.dot(fac * tau * self.u_k, ufl.grad(self.s))
* self.r
* ufl.dx
)

elif STABILIZING_METHOD == 3:
# "Finite element stabilization parameters computed from element matrices and vectors"
# Tayfun E. Tezduyar, Yasuo Osawa, 1999, starting Eq 55
# Plus the method of Mizukami-Hughes with the calculation of u_k_hat

h = ufl.CellDiameter(domain.fluid.msh)

u_mag = ufl_mag(self.u_k)

tau_s1 = h / (2.0 * u_mag)
tau_s2 = self.dt_c / 2.0
tau_s3 = h**2 / (4.0 * self.alpha_c)

tau = (
1.0 / tau_s1**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

fac = ufl.sqrt(1.0 / self.dt_c)
fac = 1.0

stab = (
ufl.dot(fac * tau * self.u_k, ufl.grad(self.s))
* self.r
* ufl.dx
)

# Form u_k_hat, the projection of the fluid velocity onto the direction of temperature gradient
grad_T = ufl.grad(self.theta_k1)
grad_T_mag = ufl_mag(grad_T)
grad_T_unit_vec = grad_T / grad_T_mag

u_k_hat = ufl.dot(self.u_k, grad_T_unit_vec) * grad_T_unit_vec
u_k_hat_mag = ufl_mag(u_k_hat)

tau_s1_alt = h / (2.0 * u_k_hat_mag)

tau_hat = (
1.0 / tau_s1_alt**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

tau_hat = tau # too diffusive
# tau_hat = ufl.max_value(0.0, tau_hat - tau)

stab += (
ufl.dot(tau_hat * u_k_hat, ufl.grad(self.s) * self.r) * ufl.dx
)

elif STABILIZING_METHOD == 4:
# "Finite element stabilization parameters computed from element matrices and vectors"
# Tayfun E. Tezduyar, Yasuo Osawa, 1999, starting Eq 55

# m_mat = 1.0/self.dt_c * ufl.inner(self.theta - self.theta_k1, self.s)
c_mat = ufl.dot(self.u_k, ufl.grad(self.theta_k1)) * self.s
# k_mat = self.alpha_c * ufl.inner(ufl.grad(self.theta), ufl.grad(self.s))
k_tilde_mat = ufl.dot(self.u_k, ufl.grad(self.theta_k1)) * ufl.dot(
self.u_k, ufl.grad(self.s)
)
c_tilde_mat = (
1.0
/ self.dt_c
* ufl.dot(
self.u_k, ufl.grad(self.s) * (self.theta_k1 - self.theta_k2)
)
)

Re_element = (
ufl_mag(self.u_k) ** 2
/ self.alpha_c
* ufl_mag(c_mat)
/ ufl_mag(k_tilde_mat)
)

tau_s1 = ufl_mag(c_mat) / ufl_mag(k_tilde_mat)
tau_s2 = self.dt_c / 2.0 * ufl_mag(c_mat) / ufl_mag(c_tilde_mat)
tau_s3 = tau_s1 * Re_element

# NVM, the above isn't working, try u_hgn
# h_ugn = ufl.CellDiameter(domain.fluid.msh)
h_ugn = (
2.0
* ufl_mag(self.u_k)
* ufl.max_value(
1e-6, ufl.dot(self.u_k, ufl.grad(self.theta_k1))
)
** (-1)
)

tau_s1 = h_ugn / (2.0 * ufl_mag(self.u_k))
tau_s2 = self.dt_c / 2.0
tau_s3 = h_ugn**2 / (4.0 * self.alpha_c)

tau = (
1.0 / tau_s1**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

fac = 1.0 # /ufl.sqrt(1.0/self.dt_c)

stab = (
ufl.dot(fac * tau * self.u_k, ufl.grad(self.s))
* self.r
* ufl.dx
)

elif STABILIZING_METHOD == 5:
# "Finite element stabilization parameters computed from element matrices and vectors"
# Tayfun E. Tezduyar, Yasuo Osawa, 1999, starting Eq 55
# plus Do Carmo and Galeao definitions of \hat{a}

h = ufl.CellDiameter(domain.fluid.msh)

u_mag = ufl_mag(self.u_k)

tau_s1 = h / (2.0 * u_mag)
tau_s2 = self.dt_c / 2.0
tau_s3 = h**2 / (4.0 * self.alpha_c)

tau = (
1.0 / tau_s1**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

fac = ufl.sqrt(1.0 / self.dt_c)
fac = 1.0

stab = (
ufl.dot(fac * tau * self.u_k, ufl.grad(self.s))
* self.r
* ufl.dx
)

# Form u_k_hat, the projection of the fluid velocity onto the direction of temperature gradient
grad_T = ufl.grad(self.theta_k1)
grad_T_mag = ufl_mag(grad_T)
grad_T_unit_vec = grad_T / grad_T_mag

u_k_hat = self.r * grad_T_unit_vec
u_k_hat_mag = ufl_mag(u_k_hat)

tau_s1_alt = h / (2.0 * u_k_hat_mag)

tau_hat = (
1.0 / tau_s1_alt**2 + 1.0 / tau_s2**2 + 1.0 / tau_s3**2
) ** (-0.5)

tau_hat = tau # too diffusive?
# tau_hat = ufl.max_value(0.0, tau_hat - tau) # original galeao and do Carmo proposal, sigma = max(0, tau(zh) - tau(b))
# tau_hat = tau * ufl.max_value(0.0, u_mag/u_k_hat_mag - 1.0) # improved/simplified, sigma = tau(b) * max(0, |b|/|z| - 1)

stab += (
ufl.dot(tau_hat * u_k_hat, ufl.grad(self.s) * self.r) * ufl.dx
)

self.F4 += stab

Expand Down Expand Up @@ -680,8 +886,9 @@ def _assemble_system(self, params):
if self.thermal_analysis:
self.solver_4 = PETSc.KSP().create(self.comm)
self.solver_4.setOperators(self.A4)
self.solver_4.setType(params.solver.solver4_ksp)
self.solver_4.getPC().setType(params.solver.solver4_pc)
self.solver_4.setType("bcgs")
self.solver_4.getPC().setType("hypre")
self.solver_4.getPC().setHYPREType("boomeramg")
self.solver_4.setFromOptions()

self.solver_5 = PETSc.KSP().create(self.comm)
Expand Down Expand Up @@ -779,6 +986,7 @@ def solve(self, domain, params, current_time):
self.u_k1.x.array[:] = self.u_k.x.array
self.p_k1.x.array[:] = self.p_k.x.array
if self.thermal_analysis:
self.theta_k2.x.array[:] = self.theta_k1.x.array
self.theta_k1.x.array[:] = self.theta_k.x.array
self.mesh_vel_old.x.array[:] = self.mesh_vel.x.array

Expand Down
2 changes: 1 addition & 1 deletion pvade/geometry/panels2d/DomainCreation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ def build_FSI(self, params):

elif np.isclose(com[1], y_min_panel):
self._add_to_domain_markers(
f"bottom_{panel_id:.0f}", [100], "facet"
f"bottom_{panel_id:.0f}", [surf_id], "facet"
)

elif np.isclose(com[1], y_max_panel):
Expand Down
Loading