# Source code for opentidalfarm.solvers.les

# Adapted from the Firedrake-Fluids project.
from dolfin import *
from dolfin_adjoint import *

[docs]
class LES(object):
r""" A solver for computing the eddy viscosity by solving:

.. math:: e = (s w)^2I

where :math:e is the eddy viscosity, :math:s is the smagorinsky
coefficient, :math:w = \sqrt{\text{cell volume}} is the filter width, and :math:I is the second
invariant defined as:

.. math:: I = \sum_{1 \le i, j \le 2} 2S_{i, j}^2, \quad S = \frac{1}{2} \left(\nabla u + \nabla u^T\right)

Parameters:

:param V: The function space for the the eddy viscosity.
:param u: The velocity function.
:param smagorinsky_coefficient: The smagorinsky coefficient.

:ivar eddy_viscosity: The smagorinsky coefficient.
"""

def __init__(self, V, u, smagorinsky_coefficient):
self._V = V
self.eddy_viscosity = Function(V)

# Create a eddy viscosity solver
les_lhs, les_rhs = self._eddy_viscosity_eqn(u, smagorinsky_coefficient)

eddy_viscosity_problem = LinearVariationalProblem(les_lhs, les_rhs,
self.eddy_viscosity, bcs=[])
self._solver = LinearVariationalSolver(eddy_viscosity_problem)
self._solver.parameters["linear_solver"] = "lu"
self._solver.parameters["symmetric"] = True
self._solver.parameters["lu_solver"]["reuse_factorization"] = True

def _strain_rate_tensor(self, u):
return S

def _eddy_viscosity_eqn(self, u, smagorinsky_coefficient):

dim = len(u)
w = TestFunction(self._V)
eddy_viscosity = TrialFunction(self._V)

cell_vol = CellVolume(self._V.mesh())
filter_width = cell_vol**(1.0/dim)

S = self._strain_rate_tensor(u)
second_invariant = 0.0
for i in range(0, dim):
for j in range(0, dim):
second_invariant += 2.0*(S[i,j]**2)

second_invariant = sqrt(second_invariant)
rhs = (smagorinsky_coefficient*filter_width)**2*second_invariant

lhs = inner(w, eddy_viscosity)*dx
rhs = inner(w, rhs)*dx

return lhs, rhs

[docs]
def solve(self):
""" Update the eddy viscosity solution for the current velocity.

:returns: The eddy viscosity.
"""
self._solver.solve()
return self.eddy_viscosity