Source code for sfepy.terms.terms_navier_stokes

import numpy as nm

from sfepy.linalg import dot_sequences, insert_strided_axis
from sfepy.terms.terms import Term, terms

[docs]class DivGradTerm(Term): r""" Diffusion term. :Definition: .. math:: \int_{\Omega} \nu\ \nabla \ul{v} : \nabla \ul{u} \mbox{ , } \int_{\Omega} \nu\ \nabla \ul{u} : \nabla \ul{w} \\ \int_{\Omega} \nabla \ul{v} : \nabla \ul{u} \mbox{ , } \int_{\Omega} \nabla \ul{u} : \nabla \ul{w} :Arguments 1: - material : :math:`\nu` (viscosity, optional) - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` :Arguments 2: - material : :math:`\nu` (viscosity, optional) - parameter_1 : :math:`\ul{u}` - parameter_2 : :math:`\ul{w}` """ name = 'dw_div_grad' arg_types = (('opt_material', 'virtual', 'state'), ('opt_material', 'parameter_1', 'parameter_2')) arg_shapes = [{'opt_material' : '1, 1', 'virtual' : ('D', 'state'), 'state' : 'D', 'parameter_1' : 'D', 'parameter_2' : 'D'}, {'opt_material' : None}] modes = ('weak', 'eval') function = staticmethod(terms.term_ns_asm_div_grad)
[docs] def d_div_grad(self, out, grad1, grad2, mat, vg, fmode): sh = grad1.shape g1 = grad1.reshape((sh[0], sh[1], sh[2] * sh[3])) g2 = grad2.reshape((sh[0], sh[1], sh[2] * sh[3])) aux = mat * dot_sequences(g1[..., None], g2, 'ATB')[..., None] if fmode == 2: out[:] = aux status = 0 else: status = vg.integrate(out, aux, fmode) return status
[docs] def get_fargs(self, mat, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(state) if mat is None: n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state) mat = nm.ones((1, n_qp, 1, 1), dtype=nm.float64) if mode == 'weak': if diff_var is None: grad = self.get(state, 'grad').transpose((0, 1, 3, 2)) sh = grad.shape grad = grad.reshape((sh[0], sh[1], sh[2] * sh[3], 1)) fmode = 0 else: grad = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return grad, mat, vg, fmode elif mode == 'eval': grad1 = self.get(virtual, 'grad') grad2 = self.get(state, 'grad') fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return grad1, grad2, mat, vg, fmode else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs] def get_eval_shape(self, mat, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(state) return (n_el, 1, 1, 1), state.dtype
[docs] def set_arg_types(self): if self.mode == 'weak': self.function = terms.term_ns_asm_div_grad else: self.function = self.d_div_grad
[docs]class ConvectTerm(Term): r""" Nonlinear convective term. :Definition: .. math:: \int_{\Omega} ((\ul{u} \cdot \nabla) \ul{u}) \cdot \ul{v} :Arguments: - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_convect' arg_types = ('virtual', 'state') arg_shapes = {'virtual' : ('D', 'state'), 'state' : 'D'} function = staticmethod(terms.term_ns_asm_convect)
[docs] def get_fargs(self, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(state) grad = self.get(state, 'grad').transpose((0, 1, 3, 2)).copy() val_qp = self.get(state, 'val') fmode = diff_var is not None return grad, val_qp, vg, fmode
[docs]class LinearConvectTerm(Term): r""" Linearized convective term. :Definition: .. math:: \int_{\Omega} ((\ul{b} \cdot \nabla) \ul{u}) \cdot \ul{v} .. math:: ((\ul{b} \cdot \nabla) \ul{u})|_{qp} :Arguments: - virtual : :math:`\ul{v}` - parameter : :math:`\ul{b}` - state : :math:`\ul{u}` """ name = 'dw_lin_convect' arg_types = ('virtual', 'parameter', 'state') arg_shapes = {'virtual' : ('D', 'state'), 'parameter' : 'D', 'state' : 'D'} function = staticmethod(terms.dw_lin_convect)
[docs] def get_fargs(self, virtual, parameter, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(state) val_qp = self.get(parameter, 'val') if mode == 'weak': if diff_var is None: grad = self.get(state, 'grad').transpose((0, 1, 3, 2)).copy() fmode = 0 else: grad = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return grad, val_qp, vg, fmode elif mode == 'qp': grad = self.get(state, 'grad').transpose((0, 1, 3, 2)).copy() fmode = 2 return grad, val_qp, vg, fmode else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs]class LinearConvect2Term(Term): r""" Linearized convective term with the convection velocity given as a material parameter. :Definition: .. math:: \int_{\Omega} ((\ul{b} \cdot \nabla) \ul{u}) \cdot \ul{v} .. math:: ((\ul{b} \cdot \nabla) \ul{u})|_{qp} :Arguments: - material : :math:`\ul{b}` - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_lin_convect2' arg_types = ('material', 'virtual', 'state') arg_shapes = {'material' : 'D, 1', 'virtual' : ('D', 'state'), 'state' : 'D'} function = staticmethod(terms.dw_lin_convect)
[docs] def get_fargs(self, material, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(state) if mode == 'weak': if diff_var is None: grad = self.get(state, 'grad').transpose((0, 1, 3, 2)).copy() fmode = 0 else: grad = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return grad, material, vg, fmode elif mode == 'qp': grad = self.get(state, 'grad').transpose((0, 1, 3, 2)).copy() fmode = 2 return grad, material, vg, fmode else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs]class StokesTerm(Term): r""" Stokes problem coupling term. Corresponds to weak forms of gradient and divergence terms. Can be evaluated. :Definition: .. math:: \int_{\Omega} p\ \nabla \cdot \ul{v} \mbox{ , } \int_{\Omega} q\ \nabla \cdot \ul{u} \mbox{ or } \int_{\Omega} c\ p\ \nabla \cdot \ul{v} \mbox{ , } \int_{\Omega} c\ q\ \nabla \cdot \ul{u} :Arguments 1: - material : :math:`c` (optional) - virtual : :math:`\ul{v}` - state : :math:`p` :Arguments 2: - material : :math:`c` (optional) - state : :math:`\ul{u}` - virtual : :math:`q` :Arguments 3: - material : :math:`c` (optional) - parameter_v : :math:`\ul{u}` - parameter_s : :math:`p` """ name = 'dw_stokes' arg_types = (('opt_material', 'virtual', 'state'), ('opt_material', 'state', 'virtual'), ('opt_material', 'parameter_v', 'parameter_s')) arg_shapes = [{'opt_material' : '1, 1', 'virtual/grad' : ('D', None), 'state/grad' : 1, 'virtual/div' : (1, None), 'state/div' : 'D', 'parameter_v' : 'D', 'parameter_s' : 1}, {'opt_material' : None}] modes = ('grad', 'div', 'eval')
[docs] @staticmethod def d_eval(out, coef, vec_qp, div, vvg): out_qp = coef * vec_qp * div status = vvg.integrate(out, out_qp) return status
[docs] def get_fargs(self, coef, vvar, svar, mode=None, term_mode=None, diff_var=None, **kwargs): if self.mode == 'grad': qp_var, qp_name = svar, 'val' else: qp_var, qp_name = vvar, 'div' n_el, n_qp, dim, n_en, n_c = self.get_data_shape(vvar) if coef is None: coef = nm.ones((1, n_qp, 1, 1), dtype=nm.float64) if mode == 'weak': vvg, _ = self.get_mapping(vvar) svg, _ = self.get_mapping(svar) if diff_var is None: val_qp = self.get(qp_var, qp_name) fmode = 0 else: val_qp = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return coef, val_qp, svg, vvg, fmode elif mode == 'eval': vvg, _ = self.get_mapping(vvar) div = self.get(vvar, 'div') vec_qp = self.get(svar, 'val') return coef, vec_qp, div, vvg else: raise ValueError('unsupported evaluation mode in %s! (%s)' % (self.name, mode))
[docs] def get_eval_shape(self, coef, vvar, svar, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(vvar) return (n_el, 1, 1, 1), vvar.dtype
[docs] def set_arg_types(self): self.function = { 'grad' : terms.dw_grad, 'div' : terms.dw_div, 'eval' : self.d_eval, }[self.mode]
[docs]class GradTerm(Term): r""" Evaluate gradient of a scalar or vector field. Supports 'eval', 'el_avg' and 'qp' evaluation modes. :Definition: .. math:: \int_{\Omega} \nabla p \mbox{ or } \int_{\Omega} \nabla \ul{w} .. math:: \mbox{vector for } K \from \Ical_h: \int_{T_K} \nabla p / \int_{T_K} 1 \mbox{ or } \int_{T_K} \nabla \ul{w} / \int_{T_K} 1 .. math:: (\nabla p)|_{qp} \mbox{ or } \nabla \ul{w}|_{qp} :Arguments: - parameter : :math:`p` or :math:`\ul{w}` """ name = 'ev_grad' arg_types = ('parameter',) arg_shapes = {'parameter' : 'N'}
[docs] @staticmethod def function(out, grad, vg, fmode): if fmode == 2: out[:] = grad status = 0 else: status = vg.integrate(out, grad, fmode) return status
[docs] def get_fargs(self, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(parameter) grad = self.get(parameter, 'grad', integration=self.integration) fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return grad, vg, fmode
[docs] def get_eval_shape(self, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(parameter) if mode != 'qp': n_qp = 1 return (n_el, n_qp, dim, n_c), parameter.dtype
[docs]class SurfaceGradTerm(GradTerm): __doc__ = GradTerm.__doc__.replace('Omega', 'Gamma') name = 'ev_surface_grad' integration = 'surface_extra'
[docs]class DivTerm(Term): r""" Evaluate divergence of a vector field. Supports 'eval', 'el_avg' and 'qp' evaluation modes. :Definition: .. math:: \int_{\Omega} \nabla \cdot \ul{u} .. math:: \mbox{vector for } K \from \Ical_h: \int_{T_K} \nabla \cdot \ul{u} / \int_{T_K} 1 .. math:: (\nabla \cdot \ul{u})|_{qp} :Arguments: - parameter : :math:`\ul{u}` """ name = 'ev_div' arg_types = ('parameter',) arg_shapes = {'parameter' : 'D'}
[docs] @staticmethod def function(out, div, vg, fmode): if fmode == 2: out[:] = div status = 0 else: status = vg.integrate(out, div, fmode) return status
[docs] def get_fargs(self, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(parameter) div = self.get(parameter, 'div', integration=self.integration) fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return div, vg, fmode
[docs] def get_eval_shape(self, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(parameter) if mode != 'qp': n_qp = 1 return (n_el, n_qp, 1, 1), parameter.dtype
[docs]class SurfaceDivTerm(DivTerm): __doc__ = DivTerm.__doc__.replace('Omega', 'Gamma') name = 'ev_surface_div' integration = 'surface_extra'
[docs]class DivOperatorTerm(Term): r""" Weighted divergence term of a test function. :Definition: .. math:: \int_{\Omega} \nabla \cdot \ul{v} \mbox { or } \int_{\Omega} c \nabla \cdot \ul{v} :Arguments: - material : :math:`c` (optional) - virtual : :math:`\ul{v}` """ name = 'dw_div' arg_types = ('opt_material', 'virtual') arg_shapes = [{'opt_material' : '1, 1', 'virtual' : ('D', None)}, {'opt_material' : None}]
[docs] @staticmethod def function(out, mat, vg): div_bf = vg.bfg n_el, n_qp, dim, n_ep = div_bf.shape div_bf = div_bf.reshape((n_el, n_qp, dim * n_ep, 1)) div_bf = nm.ascontiguousarray(div_bf) if mat is not None: status = vg.integrate(out, mat * div_bf) else: status = vg.integrate(out, div_bf) return status
[docs] def get_fargs(self, mat, virtual, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(virtual) return mat, vg
[docs]class StokesWaveTerm(Term): r""" Stokes dispersion term with the wave vector :math:`\ul{\kappa}`. :Definition: .. math:: \int_{\Omega} (\ul{\kappa} \cdot \ul{v}) (\ul{\kappa} \cdot \ul{u}) :Arguments: - material : :math:`\ul{\kappa}` - virtual : :math:`\ul{v}` - statee : :math:`\ul{u}` """ name = 'dw_stokes_wave' arg_types = ('material', 'virtual', 'state') arg_shapes = {'material' : '.: D', 'virtual' : ('D', 'state'), 'state' : 'D'} geometries = ['2_3', '2_4', '3_4', '3_8']
[docs] @staticmethod def function(out, out_qp, geo, fmode): status = geo.integrate(out, out_qp) return status
[docs] def get_fargs(self, kappa, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): from sfepy.discrete.variables import create_adof_conn, expand_basis geo, _ = self.get_mapping(state) n_el, n_qp, dim, n_en, n_c = self.get_data_shape(virtual) ebf = expand_basis(geo.bf, dim) aux = nm.einsum('i,...ij->...j', kappa, ebf)[0, :, None, :] kebf = insert_strided_axis(aux, 0, n_el) if diff_var is None: econn = state.field.get_econn('volume', self.region) adc = create_adof_conn(nm.arange(state.n_dof, dtype=nm.int32), econn, n_c, 0) vals = state()[adc] aux = dot_sequences(kebf, vals[:, None, :, None]) out_qp = dot_sequences(kebf, aux, 'ATB') fmode = 0 else: out_qp = dot_sequences(kebf, kebf, 'ATB') fmode = 1 return out_qp, geo, fmode
[docs]class StokesWaveDivTerm(Term): r""" Stokes dispersion term with the wave vector :math:`\ul{\kappa}` and the divergence operator. :Definition: .. math:: \int_{\Omega} (\ul{\kappa} \cdot \ul{v}) (\nabla \cdot \ul{u}) \;, \int_{\Omega} (\ul{\kappa} \cdot \ul{u}) (\nabla \cdot \ul{v}) :Arguments 1: - material : :math:`\ul{\kappa}` - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` :Arguments 2: - material : :math:`\ul{\kappa}` - state : :math:`\ul{u}` - virtual : :math:`\ul{v}` """ name = 'dw_stokes_wave_div' arg_types = (('material', 'virtual', 'state'), ('material', 'state', 'virtual')) arg_shapes = {'material' : '.: D', 'virtual' : ('D', 'state'), 'state' : 'D'} geometries = ['2_3', '2_4', '3_4', '3_8'] modes = ('kd', 'dk')
[docs] @staticmethod def function(out, out_qp, geo, fmode): status = geo.integrate(out, out_qp) return status
[docs] def get_fargs(self, kappa, kvar, dvar, mode=None, term_mode=None, diff_var=None, **kwargs): from sfepy.discrete.variables import create_adof_conn, expand_basis geo, _ = self.get_mapping(dvar) n_el, n_qp, dim, n_en, n_c = self.get_data_shape(kvar) ebf = expand_basis(geo.bf, dim) aux = nm.einsum('i,...ij->...j', kappa, ebf)[0, :, None, :] kebf = insert_strided_axis(aux, 0, n_el) div_bf = geo.bfg div_bf = div_bf.reshape((n_el, n_qp, 1, dim * n_en)) div_bf = nm.ascontiguousarray(div_bf) if diff_var is None: avar = dvar if self.mode == 'kd' else kvar econn = avar.field.get_econn('volume', self.region) adc = create_adof_conn(nm.arange(avar.n_dof, dtype=nm.int32), econn, n_c, 0) vals = avar()[adc] if self.mode == 'kd': aux = dot_sequences(div_bf, vals[:, None, :, None]) out_qp = dot_sequences(kebf, aux, 'ATB') else: aux = dot_sequences(kebf, vals[:, None, :, None]) out_qp = dot_sequences(div_bf, aux, 'ATB') fmode = 0 else: if self.mode == 'kd': out_qp = dot_sequences(kebf, div_bf, 'ATB') else: out_qp = dot_sequences(div_bf, kebf, 'ATB') fmode = 1 return out_qp, geo, fmode
[docs]class GradDivStabilizationTerm(Term): r""" Grad-div stabilization term ( :math:`\gamma` is a global stabilization parameter). :Definition: .. math:: \gamma \int_{\Omega} (\nabla\cdot\ul{u}) \cdot (\nabla\cdot\ul{v}) :Arguments: - material : :math:`\gamma` - virtual : :math:`\ul{v}` - state : :math:`\ul{u}` """ name = 'dw_st_grad_div' arg_types = ('material', 'virtual', 'state') arg_shapes = {'material' : '1, 1', 'virtual' : ('D', 'state'), 'state' : 'D'} function = staticmethod(terms.dw_st_grad_div)
[docs] def get_fargs(self, gamma, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(state) if diff_var is None: div = self.get(state, 'div') fmode = 0 else: div = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return div, gamma, vg, fmode
from sfepy.terms.terms_diffusion import LaplaceTerm
[docs]class PSPGPStabilizationTerm(LaplaceTerm): r""" PSPG stabilization term, pressure part ( :math:`\tau` is a local stabilization parameter), alias to Laplace term dw_laplace. :Definition: .. math:: \sum_{K \in \Ical_h}\int_{T_K} \tau_K\ \nabla p \cdot \nabla q :Arguments: - material : :math:`\tau_K` - virtual : :math:`q` - state : :math:`p` """ name = 'dw_st_pspg_p'
[docs]class PSPGCStabilizationTerm(Term): r""" PSPG stabilization term, convective part ( :math:`\tau` is a local stabilization parameter). :Definition: .. math:: \sum_{K \in \Ical_h}\int_{T_K} \tau_K\ ((\ul{b} \cdot \nabla) \ul{u}) \cdot \nabla q :Arguments: - material : :math:`\tau_K` - virtual : :math:`q` - parameter : :math:`\ul{b}` - state : :math:`\ul{u}` """ name = 'dw_st_pspg_c' arg_types = ('material', 'virtual', 'parameter', 'state') arg_shapes = {'material' : '1, 1', 'virtual' : (1, None), 'parameter' : 'D', 'state' : 'D'} function = staticmethod(terms.dw_st_pspg_c)
[docs] def get_fargs(self, tau, virtual, parameter, state, mode=None, term_mode=None, diff_var=None, **kwargs): svg, _ = self.get_mapping(virtual) vvg, _ = self.get_mapping(state) val_qp = self.get(parameter, 'val') conn = state.field.get_connectivity(self.region, self.integration) if diff_var is None: fmode = 0 else: fmode = 1 return val_qp, state(), tau, svg, vvg, conn, fmode
[docs]class SUPGPStabilizationTerm(Term): r""" SUPG stabilization term, pressure part ( :math:`\delta` is a local stabilization parameter). :Definition: .. math:: \sum_{K \in \Ical_h}\int_{T_K} \delta_K\ \nabla p\cdot ((\ul{b} \cdot \nabla) \ul{v}) :Arguments: - material : :math:`\delta_K` - virtual : :math:`\ul{v}` - parameter : :math:`\ul{b}` - state : :math:`p` """ name = 'dw_st_supg_p' arg_types = ('material', 'virtual', 'parameter', 'state') arg_shapes = {'material' : '1, 1', 'virtual' : ('D', None), 'parameter' : 'D', 'state' : 1} function = staticmethod(terms.dw_st_supg_p)
[docs] def get_fargs(self, delta, virtual, parameter, state, mode=None, term_mode=None, diff_var=None, **kwargs): vvg, _ = self.get_mapping(virtual) svg, _ = self.get_mapping(state) val_qp = self.get(parameter, 'val') if diff_var is None: grad = self.get(state, 'grad') fmode = 0 else: grad = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return val_qp, grad, delta, vvg, svg, fmode
[docs]class SUPGCStabilizationTerm(Term): r""" SUPG stabilization term, convective part ( :math:`\delta` is a local stabilization parameter). :Definition: .. math:: \sum_{K \in \Ical_h}\int_{T_K} \delta_K\ ((\ul{b} \cdot \nabla) \ul{u})\cdot ((\ul{b} \cdot \nabla) \ul{v}) :Arguments: - material : :math:`\delta_K` - virtual : :math:`\ul{v}` - parameter : :math:`\ul{b}` - state : :math:`\ul{u}` """ name = 'dw_st_supg_c' arg_types = ('material', 'virtual', 'parameter', 'state') arg_shapes = {'material' : '1, 1', 'virtual' : ('D', 'state'), 'parameter' : 'D', 'state' : 'D'} function = staticmethod(terms.dw_st_supg_c)
[docs] def get_fargs(self, delta, virtual, parameter, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(virtual) val_qp = self.get(parameter, 'val') conn = virtual.field.get_connectivity(self.region, self.integration) if diff_var is None: fmode = 0 else: fmode = 1 return val_qp, state(), delta, vg, conn, fmode