Source code for opencmp.models.multi_component_ins

########################################################################################################################
# Copyright 2021 the authors (see AUTHORS file for full list).                                                         #
#                                                                                                                      #
# This file is part of OpenCMP.                                                                                        #
#                                                                                                                      #
# OpenCMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public  #
# License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any     #
# later version.                                                                                                       #
#                                                                                                                      #
# OpenCMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied        #
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more  #
# details.                                                                                                             #
#                                                                                                                      #
# You should have received a copy of the GNU Lesser General Public License along with OpenCMP. If not, see             #
# <https://www.gnu.org/licenses/>.                                                                                     #
########################################################################################################################

from ..helpers.ngsolve_ import get_special_functions
from . import INS
from ..config_functions import ConfigParser
from typing import Dict, List, Optional, Union
from ngsolve import BilinearForm, FESpace, Grad, InnerProduct, LinearForm, GridFunction, dx, \
    Parameter, CoefficientFunction
from ngsolve.comp import ProxyFunction

import ngsolve

# TODO: These might be a good reference
# https://github.com/NGSolve/modeltemplates/tree/master/templates
# https://github.com/NGSolve/modeltemplates/blob/master/introduction/introduction.ipynb


[docs]class MultiComponentINS(INS): """ A single phase multicomponent incompressible Navier-Stokes model. """
[docs] def _add_multiple_components(self, config: ConfigParser) -> None: """ Function to add multiple components to the model from the config file. Args: config: The ConfigParser """ # Name of new components self.extra_components = config.get_list(['OTHER', 'component_names'], str) # Whether each component should be added to the local error calculation new_components_in_error_calc: Dict[str, bool] = config.get_dict(['OTHER', 'component_in_error_calc'], "") # Whether each component has a time derivative associated with it new_components_in_time_deriv: Dict[str, bool] = config.get_dict(['OTHER', 'component_in_time_deriv'], "") # Ensure that the user specified the required information about each variable assert len(self.extra_components) == len(new_components_in_error_calc) assert len(self.extra_components) == len(new_components_in_time_deriv) # Number of default variables num_existing = len(self.model_components) for i in range(len(self.extra_components)): component = self.extra_components[i] self.model_components[component] = i + num_existing self.model_local_error_components[component] = new_components_in_error_calc[component] self.time_derivative_components[0][component] = new_components_in_time_deriv[component]
def _construct_fes(self) -> FESpace: # Create the FE spaces for velocity and pressure fes_total = self._contruct_fes_helper() # Iterate over each component and add a fes for each for component in self.extra_components: fes_total.append( getattr(ngsolve, self.element[component])(self.mesh, order=self.interp_ord, dirichlet=self.dirichlet_names.get(component, ''), dgjumps=self.DG) ) return FESpace(fes_total, dgjumps=self.DG) def _define_bc_types(self) -> List[str]: # TODO: For James. Neumann BC or diffusive flux BC? Which makes more sense? # TODO: If we go with a Neumann BC need to add in the diffusion coefficient in the bilinear form. return super()._define_bc_types() + ['neumann', 'total_flux', 'surface_rxn'] def _post_init(self) -> None: if self.DIM or self.DG: raise NotImplementedError('DIM and DG are not yet implemented.') super()._post_init() def _pre_init(self) -> None: # Read in and add the multiple components self._add_multiple_components(self.config) def _set_model_parameters(self) -> None: # Read in kinematic viscosity and source terms super()._set_model_parameters() # Read in diffusion coefficient self.Ds = self.model_functions.model_parameters_dict['diffusion_coefficients'] # Ensure that a value was loaded for each extra component assert len(self.Ds) == len(self.extra_components) assert len(self.f) == len(self.extra_components) + 1 # Extra +1 since velocity also has a source term.
[docs] def construct_bilinear_time_coefficient(self, U: List[ProxyFunction], V: List[ProxyFunction], dt: Parameter, time_step: int) -> List[BilinearForm]: # Calculate the hydrodynamic contribution to the bilinear terms a = super().construct_bilinear_time_coefficient(U, V, dt, time_step)[0] w = self._get_wind(U, time_step) # TODO: Not sure that this has actually been split up correctly. # Convection term is currently the only term that must always be solved implicitly. # Transport of mixture components. if self.linearize == 'Oseen': for comp in self.extra_components: # Trial and test functions c = U[self.model_components[comp]] r = V[self.model_components[comp]] # Additional convection term. a += -dt * c * InnerProduct(w, Grad(r)) * dx return [a]
[docs] def construct_bilinear_time_ODE(self, U: Union[List[ProxyFunction], List[GridFunction]], V: List[ProxyFunction], dt: Parameter = Parameter(1.0), time_step: int = 0) -> List[BilinearForm]: # Calculate the hydrodynamic contribution to the bilinear terms a = super().construct_bilinear_time_ODE(U, V, dt, time_step)[0] w = self._get_wind(U, time_step) # Get unit normal n = get_special_functions(self.mesh, self.nu)[0] # Transport of the mixture components. for comp in self.extra_components: # Trial and test functions c = U[self.model_components[comp]] r = V[self.model_components[comp]] # Domain integrals. a += dt * self.Ds[comp][time_step] * InnerProduct(Grad(c), Grad(r)) * dx # NOTE: Reaction terms can go in either the bilinear form or the linear form depending on if the reaction # coefficients are functions of the trial functions. Basic type-checking can't be used since the # reaction coefficients will always be coefficientfunctions or gridfunctions. Instead, check if the # reaction coefficients are gridfunctions (go in linear form) and if they are coefficientfunctions # expand out the coefficientfunction tree and check for the presence of "trial-function" (goes in # bilinear form). # Surface reaction # TODO: Does not work for marker in self.BC.get('surface_rxn', {}).get(comp, {}): val = self.BC['surface_rxn'][comp][marker][time_step] if isinstance(val, ProxyFunction) or (isinstance(val, CoefficientFunction) and ('trial-function' in val.__str__())): a += -dt * val * r * self._ds(marker) # Bulk reaction # TODO: Temp solution to test, implement better once parser has been updated val = self.f[comp][time_step] if isinstance(val, ProxyFunction) or (isinstance(val, CoefficientFunction) and ('trial-function' in val.__str__())): # NOTE: The negative sign is here since the source term was brought over to bilinear side, it usually # does not have one since it is on the linear term side. if self.BC.get('surface_rxn', {}).get(comp, None) is not None: marker = list(self.BC['surface_rxn'][comp].keys())[0] a += -dt * val * r * self._ds(marker) else: a += -dt * val * r * dx if self.linearize == 'Oseen': # Neumann BC for C for marker in self.BC.get('neumann', {}).get(comp, {}): h = self.BC['neumann'][comp][marker][time_step] a += -dt * r.Trace() * (h - c * InnerProduct(w, n)) * self._ds(marker) return [a]
[docs] def construct_linear(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunction]], dt: Parameter, time_step: int) -> List[LinearForm]: # Calculate the hydrodynamic contribution to the linear terms L = super().construct_linear(V, gfu_0, dt, time_step)[0] # Transport of mixture components. for comp in self.extra_components: # Test function r = V[self.model_components[comp]] # Total Flux BC for marker in self.BC.get('total_flux', {}).get(comp, {}): h = self.BC['total_flux'][comp][marker][time_step] L += dt * r.Trace() * h * self._ds(marker) # NOTE: Reaction terms can go in either the bilinear form or the linear form depending on if the reaction # coefficients are functions of the trial functions. Basic type-checking can't be used since the # reaction coefficients will always be coefficientfunctions or gridfunctions. Instead, check if the # reaction coefficients are gridfunctions (go in linear form) and if they are coefficientfunctions # expand out the coefficientfunction tree and check for the presence of "trial-function" (goes in # bilinear form). # Surface reaction for marker in self.BC.get('surface_rxn', {}).get(comp, {}): val = self.BC['surface_rxn'][comp][marker][time_step] if not isinstance(val, ProxyFunction) and (not isinstance(val, CoefficientFunction) or ('trial-function' not in val.__str__())): L += -dt * val * r * self._ds(marker) # Bulk reaction # TODO: Temp solution to test, implement better once parser has been updated val = self.f[comp][time_step] if not isinstance(val, ProxyFunction) and (not isinstance(val, CoefficientFunction) or ('trial-function' not in val.__str__())): if self.BC.get('surface_rxn', {}).get(comp, None) is not None: marker = list(self.BC['surface_rxn'][comp].keys())[0] L += dt * val * r * self._ds(marker) else: L += dt * val * r * dx return [L]
[docs] def construct_imex_explicit(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunction]], dt: Parameter, time_step: int) -> List[LinearForm]: # Calculate the hydrodynamic contribution to the linear terms L = super().construct_imex_explicit(V, gfu_0, dt, time_step)[0] # Velocity linearization term gfu_u = gfu_0[self.model_components['u']] for comp in self.extra_components: # Previous concentration result c_prev = gfu_0[self.model_components[comp]] # Test function r = V[self.model_components[comp]] # Convection term L += -dt * c_prev * InnerProduct(gfu_u, Grad(r)) * dx # Neumann BC for C # TODO: Is this implemented correctly? for marker in self.BC.get('neumann', {}).get(comp, {}): h = self.BC['neumann'][comp][marker][time_step] L += -dt * r.Trace() * h * self._ds(marker) # Check that a total flux BC has not been incorrectly specified # TODO: Is this implemented correctly? if self.BC.get('total_flux', {}).get(comp, None) is not None: raise ValueError('Cannot specify total flux for IMEX schemes since convection term is fully explicit.') return [L]