Source code for opentidalfarm.solvers.ipcs_sw_solver

import os.path

from dolfin import *
from dolfin_adjoint import *

from .. import finite_elements
from ..problems import SWProblem
from ..problems import SteadySWProblem
from ..helpers import FrozenClass
from .solver import Solver
from .les import LES

[docs]class IPCSSWSolverParameters(FrozenClass): """ A set of parameters for a :class:`IPCSSWSolver`. Performance parameters: :ivar dolfin_solver: The dictionary with parameters for the dolfin Newton solver. A list of valid entries can be printed with: .. code-block:: python info(NonlinearVariationalSolver.default_parameters(), True) By default, the MUMPS direct solver is used for the linear system. If not availabe, the default solver and preconditioner of FEniCS is used. :ivar quadrature_degree: The quadrature degree for the matrix assembly. Default: -1 (auto) :ivar cpp_flags: A list of cpp compiler options for the code generation. Default: ["-O3", "-ffast-math", "-march=native"] Large eddy simulation parameters: :ivar les_model: De-/Activates the LES model. Default: True :ivar les_model_parameters: A dictionary with parameters for the LES model. Default: {'smagorinsky_coefficient': 1e-2} """ dolfin_solver = {"newton_solver": {}} # Large eddy simulation les_model = True les_parameters = {'smagorinsky_coefficient': 1e-2} # Performance settings quadrature_degree = -1 cpp_flags = ["-O3", "-ffast-math", "-march=native"] def __init__(self): linear_solver = 'mumps' if ('mumps' in [x[0] for x in linear_solver_methods()]) else 'default' preconditioner = 'default' self.dolfin_solver["newton_solver"]["linear_solver"] = linear_solver self.dolfin_solver["newton_solver"]["preconditioner"] = preconditioner self.dolfin_solver["newton_solver"]["maximum_iterations"] = 20
[docs]class IPCSSWSolver(Solver): r""" This incremental pressure correction scheme (IPCS) is an operator splitting scheme that follows the idea of Goda [1]_ and Simo [2]_. This scheme preserves the exact same stability properties as Navier-Stokes and hence does not introduce additional dissipation in the flow (FIXME: needs verification). The idea is to replace the unknown free-surface with an approximation. This is chosen as the free-surface solution from the previous solution. The time discretization is done using a :math:`\theta`-scheme, the convection, friction and divergence are handled semi-implicitly. Thus, we have a discretized version of the shallow water equations as .. math:: \frac{1}{\Delta t}\left( u^{n+1}-u^{n} \right)-\nabla\cdot\nu\nabla u^{n+\theta}+u^*\cdot\nabla u^{n+\theta}+g\nabla \eta^{n+\theta} + \frac{c_b + c_t}{H^n} \| u^n\| u^{n+\theta} &= f_u^{n+\theta}, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \nabla \cdot \left( H^{n} u^{n+\theta} \right) &= 0, where :math:`\square^{n+\theta} = \theta{\square}^{n+1}+(1-\theta)\square^n, \theta \in [0, 1]` and :math:`u^* = \frac{3}{2}u^n - \frac{1}{2}u^{n-1}`. This convection term is unconditionally stable, and with :math:`\theta=0.5`, this equation is second order in time and space [2]_ (FIXME: Needs verification). For the operator splitting, we use the free-surface solution from the previous timestep as an estimation, giving an equation for a tentative velocity, :math:`\tilde{u}^{n+1}`: .. math:: \frac{1}{\Delta t}\left( \tilde{u}^{n+1}-u^{n} \right) - \nabla\cdot\nu\nabla \tilde{u}^{n+\theta} + u^*\cdot\nabla \tilde u^{n+\theta}+g\nabla \eta^{n} + \frac{c_b + c_t}{H^n} \| u^n\| \tilde u^{n+\theta} = f_u^{n+\theta}. This tenative velocity does not satisfy the divergence equation, and thus we define a velocity correction :math:`u^c=u^{n+1}-\tilde{u}^{n+1}`. Substracting the second equation from the first, we see that .. math:: \frac{1}{\Delta t}u^c - \theta \nabla\cdot\nu\nabla u^c + \theta u^*\cdot\nabla u^{c} + g\theta \nabla\left( \eta^{n+1} - \eta^n\right) + \theta \frac{c_b + c_t}{H^n} \| u^n\| u^{c} &= 0, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \theta \nabla \cdot \left( H^{n} u^c \right) &= -\nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right). The operator splitting is a first order approximation, :math:`O(\Delta t)`, so we can, without reducing the order of the approximation simplify the above to .. math:: \frac{1}{\Delta t}u^c + g\theta \nabla\left( \eta^{n+1} - \eta^n\right) &= 0, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \theta \nabla \cdot \left( H^{n} u^c \right) &= -\nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right), which is reducible to the problem: .. math:: \eta^{n+1}-\eta^{n} - g \Delta t^2 \theta^2 \nabla \cdot \left( H^{n+1} \nabla\left( \eta^{n+1} - \eta^n\right) \right) = -\Delta t \nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right). The corrected velocity is then easily calculated from .. math:: u^{n+1} = \tilde{u}^{n+1} - \Delta tg\theta\nabla\left(\eta^{n+1}-\eta^n\right) The scheme can be summarized in the following steps: #. Replace the pressure with a known approximation and solve for a tenative velocity :math:`\tilde u^{n+1}`. #. Solve a free-surface correction equation for the free-surface, :math:`\eta^{n+1}` #. Use the corrected pressure to find the velocity correction and calculate :math:`u^{n+1}` #. Update t, and repeat. Remarks: - This solver only works with transient problems, that is with :class:`opentidalfarm.problems.sw.SWProblem`. - This solver supports large eddy simulation (LES). The LES model is implemented via the :class:`opentidalfarm.solvers.les.LES` class. .. [1] Goda, Katuhiko. *A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows.* Journal of Computational Physics 30.1 (1979): 76-95. .. [2] Simo, J. C., and F. Armero. *Unconditional stability and long-term behavior of transient algorithms for the incompressible Navier-Stokes and Euler equations.* Computer Methods in Applied Mechanics and Engineering 111.1 (1994): 111-154. """ def __init__(self, problem, parameters): if not isinstance(problem, (SWProblem, SteadySWProblem)): raise TypeError("problem must be of type Problem") if not isinstance(parameters, IPCSSWSolverParameters): raise TypeError("parameters must be of type \ IPCSSWSolverParameters.") super(IPCSSWSolver, self).__init__() self.problem = problem self.parameters = parameters # If cache_for_nonlinear_initial_guess is true, then we store all # intermediate state variables in this dictionary to be used for the # next solve self.state_cache = {} self.mesh = problem.parameters.domain.mesh ele_u, ele_eta = self.problem.parameters.finite_element() self.V = FunctionSpace(self.mesh, ele_u) self.Q = FunctionSpace(self.mesh, ele_eta) @staticmethod def default_parameters(): return IPCSSWSolverParameters() def _generate_strong_bcs(self, dgu): if dgu: bcu_method = "geometric" else: bcu_method = "topological" bcs = self.problem.parameters.bcs facet_ids = self.problem.parameters.domain.facet_ids # Generate velocity boundary conditions bcs_u = [] for _, expr, facet_id, _ in bcs.filter("u", "strong_dirichlet"): bc = DirichletBC(self.V, expr, facet_ids, facet_id, method=bcu_method) bcs_u.append(bc) # Generate free-surface boundary conditions bcs_eta = [] for _, expr, facet_id, _ in bcs.filter("eta", "strong_dirichlet"): bc = DirichletBC(self.Q, expr, facet_ids, facet_id) bcs_eta.append(bc) return bcs_u, bcs_eta def _finished(self, current_time, finish_time): if (hasattr(self.problem.parameters, 'finished')): return self.problem.parameters.finished(current_time) else: return float(current_time - finish_time) >= - 1e3*DOLFIN_EPS
[docs] def solve(self, annotate=True): ''' Solve the shallow water equations ''' ############################### Setting up the equations ########################### # Initialise solver settings if not type(self.problem) == SWProblem: raise TypeError("Do not know how to solve problem of type %s." % type(self.problem)) # Get parameters problem_params = self.problem.parameters solver_params = self.parameters farm = problem_params.tidal_farm if farm: turbine_friction = farm.friction_function # Performance settings parameters['form_compiler']['quadrature_degree'] = \ solver_params.quadrature_degree parameters['form_compiler']['cpp_optimize_flags'] = \ " ".join(solver_params.cpp_flags) parameters['form_compiler']['cpp_optimize'] = True parameters['form_compiler']['optimize'] = True # Get domain measures ds = problem_params.domain.ds dx = problem_params.domain.dx # Get the boundary normal direction n = FacetNormal(self.mesh) # Get temporal settings theta = Constant(problem_params.theta) dt = Constant(problem_params.dt) finish_time = Constant(problem_params.finish_time) t = Constant(problem_params.start_time) # Get equation settings g = problem_params.g h = problem_params.depth nu = problem_params.viscosity include_advection = problem_params.include_advection include_viscosity = problem_params.include_viscosity linear_divergence = problem_params.linear_divergence f_u = problem_params.f_u include_les = solver_params.les_model # Get boundary conditions bcs = problem_params.bcs # Get function spaces V, Q = self.V, self.Q dgu = "Discontinuous" in str(V) # Test and trial functions v = TestFunction(V) u = TrialFunction(V) q = TestFunction(Q) eta = TrialFunction(Q) # Functions u00 = Function(V) u0 = Function(V, name="u0") ut = Function(V) # Tentative velocity u1 = Function(V, name="u") eta0 = Function(Q, name="eta0") eta1 = Function(Q, name="eta") # Large eddy model if include_les: les_V = FunctionSpace(problem_params.domain.mesh, "CG", 1) les = LES(les_V, u0, solver_params.les_parameters['smagorinsky_coefficient']) eddy_viscosity = les.eddy_viscosity nu += eddy_viscosity else: eddy_viscosity = None # Define the water depth if linear_divergence: H = h else: H = eta0 + h if f_u is None: f_u = Constant((0, 0)) # Bottom friction friction = problem_params.friction if farm: friction += Function(turbine_friction, name="turbine_friction", annotate=annotate) # Load initial conditions # Projection is necessary to obtain 2nd order convergence # FIXME: The problem should specify the ic for each component separately. u_ic = project(problem_params.initial_condition_u, self.V) u0.assign(u_ic, annotate=False) u00.assign(u_ic, annotate=False) eta_ic = project(problem_params.initial_condition_eta, self.Q) eta0.assign(eta_ic, annotate=False) eta1.assign(eta_ic, annotate=False) # Tentative velocity step u_mean = theta * u + (1. - theta) * u0 u_bash = 3./2 * u0 - 1./2 * u00 u_diff = u - u0 norm_u0 = inner(u0, u0)**0.5 F_u_tent = ((1/dt) * inner(v, u_diff) * dx() + inner(v, grad(u_bash)*u_mean) * dx() + g * inner(v, grad(eta0)) * dx() + friction / H * norm_u0 * inner(u_mean, v) * dx - inner(v, f_u) * dx()) # Viscosity term if dgu: # Taken from sigma = 1. # Penalty parameter. # Set tau=-1 for SIPG, tau=0 for IIPG, and tau=1 for NIPG tau = 0. edgelen = FacetArea(self.mesh)('+') # Facetarea is continuous, so # we can select either side alpha = sigma/edgelen F_u_tent += nu * inner(grad(v), grad(u_mean)) * dx() for d in range(2): F_u_tent += - nu * inner(avg(grad(u_mean[d])), jump(v[d], n))*dS F_u_tent += - nu * tau * inner(avg(grad(v[d])), jump(u[d], n))*dS F_u_tent += alpha * nu * inner(jump(u[d], n), jump(v[d], n))*dS else: F_u_tent += nu * inner(grad(v), grad(u_mean)) * dx() a_u_tent = lhs(F_u_tent) L_u_tent = rhs(F_u_tent) # Pressure correction eta_diff = eta - eta0 ut_mean = theta * ut + (1. - theta) * u0 F_p_corr = (q*eta_diff + g * dt**2 * theta**2 * H * inner(grad(q), grad(eta_diff)))*dx() + dt*q*div(H*ut_mean)*dx() a_p_corr = lhs(F_p_corr) L_p_corr = rhs(F_p_corr) # Velocity correction eta_diff = eta1 - eta0 a_u_corr = inner(v, u)*dx() L_u_corr = inner(v, ut)*dx() - dt*g*theta*inner(v, grad(eta_diff))*dx() bcu, bceta = self._generate_strong_bcs(dgu) # Assemble matrices A_u_corr = assemble(a_u_corr) for bc in bcu: bc.apply(A_u_corr) a_u_corr_solver = LUSolver(A_u_corr) a_u_corr_solver.parameters["reuse_factorization"] = True if linear_divergence: A_p_corr = assemble(a_p_corr) for bc in bceta: bc.apply(A_p_corr) a_p_corr_solver = LUSolver(A_p_corr) a_p_corr_solver.parameters["reuse_factorization"] = True yield({"time": float(t), "u": u0, "eta": eta0, "eddy_viscosity": eddy_viscosity, "is_final": self._finished(t, finish_time)}) log(INFO, "Start of time loop") adjointer.time.start(t) timestep = 0 # De/activate annotation annotate_orig = parameters["adjoint"]["stop_annotating"] parameters["adjoint"]["stop_annotating"] = not annotate while not self._finished(t, finish_time): # Update timestep timestep += 1 t = Constant(t + dt) # Update bc's t_theta = Constant(t - (1.0 - theta) * dt) bcs.update_time(t, only_type=["strong_dirichlet"]) bcs.update_time(t_theta, exclude_type=["strong_dirichlet"]) # Update source term if f_u is not None: f_u.t = Constant(t_theta) if include_les: log(PROGRESS, "Compute eddy viscosity.") les.solve() # Compute tentative velocity step log(PROGRESS, "Solve for tentative velocity.") A_u_tent = assemble(a_u_tent) b = assemble(L_u_tent) for bc in bcu: bc.apply(A_u_tent, b) solve(A_u_tent, ut.vector(), b) # Pressure correction log(PROGRESS, "Solve for pressure correction.") b = assemble(L_p_corr) for bc in bceta: bc.apply(b) if linear_divergence: a_p_corr_solver.solve(eta1.vector(), b) else: A_p_corr = assemble(a_p_corr) for bc in bceta: bc.apply(A_p_corr) solve(A_p_corr, eta1.vector(), b) # Velocity correction log(PROGRESS, "Solve for velocity update.") b = assemble(L_u_corr) for bc in bcu: bc.apply(b) a_u_corr_solver.solve(u1.vector(), b) # Rotate functions for next timestep u00.assign(u0) u0.assign(u1) eta0.assign(eta1) # Increase the adjoint timestep adj_inc_timestep(time=float(t), finished=self._finished(t, finish_time)) yield({"time": float(t), "u": u0, "eta": eta0, "eddy_viscosity": eddy_viscosity, "is_final": self._finished(t, finish_time)}) # Reset annotation flag parameters["adjoint"]["stop_annotating"] = annotate_orig log(INFO, "End of time loop.")