Source code for sfepy.terms.terms_diffusion

import numpy as nm

from sfepy.base.base import assert_
from sfepy.linalg import dot_sequences
from sfepy.terms.terms import Term, terms
from sfepy.terms.terms_dot import ScalarDotMGradScalarTerm

[docs] class DiffusionTerm(Term): r""" General diffusion term with permeability :math:`K_{ij}`. Can be evaluated. Can use derivatives. :Definition: .. math:: \int_{\Omega} K_{ij} \nabla_i q \nabla_j p :Arguments: - material: :math:`K_{ij}` - virtual/parameter_1: :math:`q` - state/parameter_2: :math:`p` """ name = 'dw_diffusion' arg_types = (('material', 'virtual', 'state'), ('material', 'parameter_1', 'parameter_2')) arg_shapes = {'material' : 'D, D', 'virtual' : (1, 'state'), 'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1} modes = ('weak', 'eval') symbolic = {'expression': 'div( K * grad( u ) )', 'map' : {'u' : 'state', 'K' : 'material'}}
[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: if self.name in ('dw_laplace', 'dw_st_pspg_p'): n_el, n_qp, _, _, _ = 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') 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') return grad1, grad2, mat, vg 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.dw_diffusion else: self.function = terms.d_diffusion
[docs] class SDDiffusionTerm(Term): r""" Diffusion sensitivity analysis term. :Definition: .. math:: \int_{\Omega} \hat{K}_{ij} \nabla_i q\, \nabla_j p .. math:: \hat{K}_{ij} = K_{ij}\left( \delta_{ik}\delta_{jl} \nabla \cdot \ul{\Vcal} - \delta_{ik}{\partial \Vcal_j \over \partial x_l} - \delta_{jl}{\partial \Vcal_i \over \partial x_k}\right) :Arguments: - material: :math:`K_{ij}` - parameter_q: :math:`q` - parameter_p: :math:`p` - parameter_mv: :math:`\ul{\Vcal}` """ name = 'ev_sd_diffusion' arg_types = ('material', 'parameter_q', 'parameter_p', 'parameter_mv') arg_shapes = {'material' : 'D, D', 'parameter_q' : 1, 'parameter_p' : 1, 'parameter_mv' : 'D'} function = staticmethod(terms.d_sd_diffusion)
[docs] def get_fargs(self, mat, parameter_q, parameter_p, parameter_mv, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(parameter_p) grad_q = self.get(parameter_q, 'grad') grad_p = self.get(parameter_p, 'grad') grad_mv = self.get(parameter_mv, 'grad') div_mv = self.get(parameter_mv, 'div') return grad_q, grad_p, grad_mv, div_mv, mat, vg
[docs] def get_eval_shape(self, mat, parameter_q, parameter_p, parameter_mv, mode=None, term_mode=None, diff_var=None, **kwargs): n_el, n_qp, dim, n_en, n_c = self.get_data_shape(parameter_q) return (n_el, 1, 1, 1), parameter_q.dtype
[docs] class LaplaceTerm(DiffusionTerm): r""" Laplace term with :math:`c` coefficient. Can be evaluated. Can use derivatives. :Definition: .. math:: \int_{\Omega} c \nabla q \cdot \nabla p :Arguments 1: - material: :math:`c` - virtual/parameter_1: :math:`q` - state/parameter_2: :math:`p` """ name = 'dw_laplace' arg_types = (('opt_material', 'virtual', 'state'), ('opt_material', 'parameter_1', 'parameter_2')) arg_shapes = [{'opt_material' : '1, 1', 'virtual' : (1, 'state'), 'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1}, {'opt_material' : None}] modes = ('weak', 'eval') symbolic = {'expression': 'c * div( grad( u ) )', 'map' : {'u' : 'state', 'c' : 'opt_material'}}
[docs] def set_arg_types(self): if self.mode == 'weak': self.function = terms.dw_laplace else: self.function = terms.d_laplace
[docs] class DiffusionRTerm(Term): r""" Diffusion-like term with material parameter :math:`K_{j}` (to use on the right-hand side). :Definition: .. math:: \int_{\Omega} K_{j} \nabla_j q :Arguments: - material : :math:`K_j` - virtual : :math:`q` """ name = 'dw_diffusion_r' arg_types = ('material', 'virtual') arg_shapes = {'material' : 'D, 1', 'virtual' : (1, None)} function = staticmethod(terms.dw_diffusion_r)
[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 DiffusionCoupling(Term): r""" Diffusion copupling term with material parameter :math:`K_{j}`. :Definition: .. math:: \int_{\Omega} p K_{j} \nabla_j q \mbox{ , } \int_{\Omega} q K_{j} \nabla_j p :Arguments: - material : :math:`K_{j}` - virtual : :math:`q` - state : :math:`p` """ name = 'dw_diffusion_coupling' arg_types = (('material', 'virtual', 'state'), ('material', 'state', 'virtual'), ('material', 'parameter_1', 'parameter_2')) arg_shapes = {'material' : 'D, 1', 'virtual' : (1, 'state'), 'state' : 1, 'parameter_1' : 1, 'parameter_2' : 1} modes = ('weak0', 'weak1', 'eval')
[docs] @staticmethod def d_fun(out, mat, val, grad, vg): out_qp = val * dot_sequences(mat, grad, 'ATB') status = vg.integrate(out, out_qp) return status
[docs] @staticmethod def dw_fun(out, val, mat, bf, vg, fmode): if fmode == 0: status = terms.mulAB_integrate(out, vg.bfg, mat * val, vg.cmap, mode='ATB') elif fmode == 1: status = terms.mulAB_integrate(out, bf * mat, val, vg.cmap, mode='ATB') elif fmode == 2: status = terms.mulAB_integrate(out, vg.bfg, mat * bf, vg.cmap, mode='ATB') elif fmode == 3: status = terms.mulAB_integrate(out, mat * bf, vg.bfg, vg.cmap, mode='ATB') return status
[docs] def get_fargs( self, mat, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): vg, _ = self.get_mapping(virtual) if mode == 'weak': vgs, _ = self.get_mapping(state) if diff_var is None: if self.mode == 'weak0': val = self.get(state, 'val') fmode = 0 else: val = self.get(virtual, 'grad') fmode = 1 else: val = nm.array([0], ndmin=4, dtype=nm.float64) if self.mode == 'weak0': fmode = 2 else: fmode = 3 return val, mat, vgs.bf, vg, fmode elif mode == 'eval': grad = self.get(virtual, 'grad') val = self.get(state, 'val') return mat, val, grad, vg 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[:-1] == 'weak': self.function = self.dw_fun else: self.function = self.d_fun
[docs] class DiffusionVelocityTerm( Term ): r""" Evaluate diffusion velocity. Supports 'eval', 'el_avg' and 'qp' evaluation modes. :Definition: .. math:: - \int_{\cal{D}} K_{ij} \nabla_j p :Arguments: - material : :math:`K_{ij}` - parameter : :math:`p` """ name = 'ev_diffusion_velocity' arg_types = ('material', 'parameter') arg_shapes = {'material' : 'D, D', 'parameter' : 1} integration = ('cell', 'facet_extra')
[docs] @staticmethod def function(out, grad, mat, vg, fmode): dvel = dot_sequences(mat, grad) if fmode == 2: out[:] = dvel status = 0 else: status = vg.integrate(out, dvel, fmode) out *= -1.0 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') fmode = {'eval' : 0, 'el_avg' : 1, 'qp' : 2}.get(mode, 1) return grad, mat, 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, 1), parameter.dtype
[docs] class SurfaceFluxTerm(Term): r""" Surface flux term. Supports 'eval', 'el_eval' and 'el_avg' evaluation modes. :Definition: .. math:: \int_{\Gamma} \ul{n} \cdot K_{ij} \nabla_j p :Arguments: - material: :math:`\ul{K}` - parameter: :math:`p`, """ name = 'ev_surface_flux' arg_types = ('material', 'parameter') arg_shapes = {'material' : 'D, D', 'parameter' : 1} integration = 'facet_extra' function = staticmethod(terms.d_surface_flux)
[docs] def get_fargs(self, mat, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(parameter) grad = self.get(parameter, 'grad') fmode = {'eval' : 0, 'el_avg' : 1}.get(mode, 0) return grad, mat, sg, fmode
[docs] def get_eval_shape(self, mat, parameter, mode=None, term_mode=None, diff_var=None, **kwargs): n_fa, n_qp, dim, n_en, n_c = self.get_data_shape(parameter) return (n_fa, 1, 1, 1), parameter.dtype
[docs] class SurfaceFluxOperatorTerm(Term): r""" Surface flux operator term. :Definition: .. math:: \int_{\Gamma} q \ul{n} \cdot \ull{K} \cdot \nabla p :Arguments: - material : :math:`\ull{K}` - virtual : :math:`q` - state : :math:`p` """ name = 'dw_surface_flux' arg_types = ('opt_material', 'virtual', 'state') arg_shapes = [{'opt_material' : 'D, D', 'virtual' : (1, 'state'), 'state' : 1}, {'opt_material' : None}] integration = 'facet_extra' function = staticmethod(terms.dw_surface_flux)
[docs] def get_fargs(self, mat, virtual, state, mode=None, term_mode=None, diff_var=None, **kwargs): sg, _ = self.get_mapping(state) sd = state.field.extra_data[f'sd_{self.region.name}'] self.get_mapping(virtual) # Creates BQP for the basis. bf = virtual.field.eval_basis(sd.bkey, 0, self.integral) if mat is None: _, n_qp, dim, _, _ = self.get_data_shape(state) mat = nm.empty((1, n_qp, dim, dim), dtype=nm.float64) mat[..., :, :] = nm.eye(dim, dtype=nm.float64) if diff_var is None: grad = self.get(state, 'grad', integration='facet_extra') fmode = 0 else: grad = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 1 return grad, mat, bf, sg, sd.fis, fmode
[docs] class ConvectVGradSTerm(Term): r""" Scalar gradient term with convective velocity. :Definition: .. math:: \int_{\Omega} q (\ul{u} \cdot \nabla p) :Arguments: - virtual : :math:`q` - state_v : :math:`\ul{u}` - state_s : :math:`p` """ name = 'dw_convect_v_grad_s' arg_types = ('virtual', 'state_v', 'state_s') arg_shapes = [{'virtual' : (1, 'state_s'), 'state_v' : 'D', 'state_s' : 1}] function = staticmethod(terms.dw_convect_v_grad_s)
[docs] def get_fargs(self, virtual, state_v, state_s, mode=None, term_mode=None, diff_var=None, **kwargs): vvg, _ = self.get_mapping(state_v) svg, _ = self.get_mapping(state_s) if diff_var is None: grad_s = self.get(state_s, 'grad') val_v = self.get(state_v, 'val') fmode = 0 elif diff_var == state_s.name: grad_s = nm.array([0], ndmin=4, dtype=nm.float64) val_v = self.get(state_v, 'val') fmode = 1 else: assert_(diff_var == state_v.name) grad_s = self.get(state_s, 'grad') val_v = nm.array([0], ndmin=4, dtype=nm.float64) fmode = 2 return val_v, grad_s, vvg, svg, fmode
[docs] class AdvectDivFreeTerm(ScalarDotMGradScalarTerm): r""" Advection of a scalar quantity :math:`p` with the advection velocity :math:`\ul{y}` given as a material parameter (a known function of space and time). The advection velocity has to be divergence-free! :Definition: .. math:: \int_{\Omega} \nabla \cdot (\ul{y} p) q = \int_{\Omega} (\underbrace{(\nabla \cdot \ul{y})}_{\equiv 0} + \ul{y} \cdot \nabla) p) q :Arguments: - material : :math:`\ul{y}` - virtual : :math:`q` - state : :math:`p` """ name = 'dw_advect_div_free' arg_types = ('material', 'virtual', 'state') arg_shapes = {'material' : 'D, 1', 'virtual' : ('1', 'state'), 'state' : '1'} mode = 'grad_state'
[docs] class NonlinearDiffusionTerm(Term): r""" The diffusion term with a scalar coefficient given by a user supplied function of the state variable. :Definition: .. math:: \int_{\Omega} \nabla q \cdot \nabla p f(p) :Arguments: - fun : :math:`f(p)` - dfun : :math:`\partial f(p) / \partial p` - virtual : :math:`q` - state : :math:`p` """ name = 'dw_nl_diffusion' arg_types = ('fun', 'dfun', 'virtual', 'state') arg_shapes = {'fun' : lambda x: x, 'dfun' : lambda x: x, 'virtual' : (1, 'state'), 'state' : 1}
[docs] @staticmethod def function(out, out_qp, geo): status = geo.integrate(out, out_qp) return status
[docs] def get_fargs(self, fun, dfun, var1, var2, mode=None, term_mode=None, diff_var=None, **kwargs): vg1, _ = self.get_mapping(var1) vg2, _ = self.get_mapping(var2) if diff_var is None: geo = vg1 val_grad_qp = self.get(var2, 'grad') val_qp = fun(self.get(var2, 'val')) out_qp = dot_sequences(vg1.bfg, val_grad_qp*val_qp,'ATB') else: geo = vg1 val_grad_qp = self.get(var2, 'grad') val_d_qp = dfun(self.get(var2, 'val')) val_qp = fun(self.get(var2, 'val')) out_qp = (dot_sequences(vg1.bfg, vg2.bfg*val_qp,'ATB') + dot_sequences(vg1.bfg, val_grad_qp*val_d_qp,'ATB')*vg2.bf) return out_qp, geo