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} \nabla \ul{v} : \nabla \ul{u} :Arguments: - material: :math:`\nu` (viscosity, optional) - virtualparameter_1: :math:`\ul{v}` - state/parameter_2: :math:`\ul{u}` """ 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])) sh = grad2.shape 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): vgs, _ = self.get_mapping(state) vgv, _ = self.get_mapping(virtual) 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, vgv, vgs, 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, vgv, 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{w} \cdot \nabla) \ul{u}) \cdot \ul{v} .. math:: ((\ul{w} \cdot \nabla) \ul{u})|_{qp} :Arguments: - virtual : :math:`\ul{v}` - parameter : :math:`\ul{w}` - 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{c} \cdot \nabla) \ul{u}) \cdot \ul{v} .. math:: ((\ul{c} \cdot \nabla) \ul{u})|_{qp} :Arguments: - material : :math:`\ul{c}` - 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/parameter_v: :math:`\ul{v}` - state/parameter_s: :math:`p` :Arguments 2: - material : :math:`c` (optional) - state : :math:`\ul{u}` - virtual : :math:`q` """ 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_{\cal{D}} \nabla p \mbox{ or } \int_{\cal{D}} \nabla \ul{u}\\ \int_{\cal{D}} c \nabla p \mbox{ or } \int_{\cal{D}} c \nabla \ul{u} :Arguments: - parameter : :math:`p` or :math:`\ul{u}` """ name = 'ev_grad' arg_types = ('opt_material', 'parameter') arg_shapes = [{'opt_material': '1, 1', 'parameter': 'N'}, {'opt_material': None}] integration = ('cell', 'facet_extra')
[docs] @staticmethod def function(out, mat, grad, vg, fmode): if fmode == 2: out[:] = grad if mat is None else mat * grad status = 0 else: status = vg.integrate(out, grad, fmode) if mat is None else\ vg.integrate(out, mat * grad, fmode) return status
[docs] def get_fargs(self, mat, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(parameter) grad = self.get(parameter, 'grad', integration=self.act_integration) fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return mat, grad, vg, fmode
[docs] def get_eval_shape(self, mat, 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 DivTerm(Term): r""" Evaluate divergence of a vector field. Supports 'eval', 'el_avg' and 'qp' evaluation modes. :Definition: .. math:: \int_{\cal{D}} \nabla \cdot \ul{u} \mbox { , } \int_{\cal{D}} c \nabla \cdot \ul{u} :Arguments: - parameter : :math:`\ul{u}` """ name = 'ev_div' arg_types = ('opt_material', 'parameter') arg_shapes = [{'opt_material': '1, 1', 'parameter': 'D'}, {'opt_material': None}] integration = ('cell', 'facet_extra')
[docs] @staticmethod def function(out, mat, div, vg, fmode): if fmode == 2: out[:] = div if mat is None else mat * div status = 0 else: status = vg.integrate(out, div, fmode) if mat is None else\ vg.integrate(out, mat * div, fmode) return status
[docs] def get_fargs(self, mat, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(parameter) div = self.get(parameter, 'div', integration=self.act_integration) fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return mat, div, vg, fmode
[docs] def get_eval_shape(self, mat, 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 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('cell', 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('cell', 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_econn(self.act_integration, self.region) 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_econn(self.act_integration, self.region) if diff_var is None: fmode = 0 else: fmode = 1 return val_qp, state(), delta, vg, conn, fmode