burgers_assembler_carreau.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <algorithm>
#include <kernel/assembly/asm_traits.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/global/vector.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename DataType_, typename IndexType_, int dim_>
    class BurgersAssemblerCarreau
    {
    public:
      typedef DataType_ DataType;
 
      bool deformation;
 
 
      bool diff;
 
      bool frechet_diff;
 
      //DataType_ nu;
 
      DataType_ theta;
 
      DataType_ beta;
 
      DataType_ frechet_beta;
 
      DataType_ sd_delta;
 
      DataType_ sd_nu;
 
      DataType_ sd_v_norm;
 
      DataType_ mu_inf, mu_0, lambda, n, a, reg_eps;
 
      BurgersAssemblerCarreau() :
        deformation(true),
        diff(true),
        frechet_diff(false),
        theta(DataType_(0)),
        beta(DataType_(0)),
        frechet_beta(DataType_(0)),
        sd_delta(DataType_(0)),
        sd_nu(DataType_(0)),
        sd_v_norm(DataType_(0)),
        mu_inf(DataType_(0.0)),
        mu_0(DataType_(0.01)),
        lambda(DataType_(1.0)),
        n(DataType_(0.31)),
        a(DataType_(2.0)),
        reg_eps(DataType_(1e-100))
      {
      }
 
      template<typename Space_, typename CubatureFactory_>
      void assemble_matrix(
        LAFEM::SparseMatrixBCSR<DataType_, IndexType_, dim_, dim_>& matrix,
        const LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_>& convect,
        const Space_& space,
        const CubatureFactory_& cubature_factory,
        const DataType_ scale = DataType_(1)
        ) const
      {
        // validate matrix and vector dimensions
        XASSERTM(matrix.rows() == space.get_num_dofs(), "invalid matrix dimensions");
        XASSERTM(matrix.columns() == space.get_num_dofs(), "invalid matrix dimensions");
        XASSERTM(convect.size() == space.get_num_dofs(), "invalid vector size");
 
        typedef LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_> VectorType;
        typedef LAFEM::SparseMatrixBCSR<DataType_, IndexType_, dim_, dim_> MatrixType;
 
        // tolerance: sqrt(eps)
        const DataType tol_eps = Math::sqrt(Math::eps<DataType>());
 
        // first of all, let's see what we have to assemble
 
        const bool need_conv = (Math::abs(beta) > DataType(0));
        const bool need_conv_frechet = (Math::abs(frechet_beta) > DataType(0));
        const bool need_reac = (Math::abs(theta) > DataType(0));
        const bool need_streamdiff = (Math::abs(sd_delta) > DataType(0)) && (sd_v_norm > tol_eps);
 
        // define our assembly traits
        typedef AsmTraits1<DataType_, Space_, TrafoTags::jac_det, SpaceTags::value|SpaceTags::grad> AsmTraits;
 
        // fetch our trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // create cubature rule
        typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
 
        // create matrix scatter-axpy
        typename MatrixType::ScatterAxpy scatter_matrix(matrix);
 
        // create convection gather-axpy
        typename VectorType::GatherAxpy gather_conv(convect);
 
        // get maximum number of local dofs
        static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
        // create local matrix data
        typedef Tiny::Matrix<DataType, dim_, dim_> MatrixValue;
        typedef Tiny::Matrix<MatrixValue, max_local_dofs, max_local_dofs> LocalMatrixType;
        LocalMatrixType local_matrix;
 
        // create local vector data
        typedef Tiny::Vector<DataType, dim_> VectorValue;
        typedef Tiny::Vector<VectorValue, max_local_dofs> LocalVectorType;
 
        // local convection field dofs
        LocalVectorType local_conv_dofs;
 
        // our local velocity value
        VectorValue loc_v, mean_v;
 
        // our local velocity gradient
        MatrixValue loc_grad_v;
        MatrixValue strain_rate_tensor_2;
 
        loc_v.format();
        mean_v.format();
        loc_grad_v.format();
 
        // compute reference element barycentre
        Tiny::Vector<DataType, dim_> barycentre;
        for(int i(0); i < dim_; ++i)
          barycentre[i] = Shape::ReferenceCell<typename Space_::ShapeType>::template centre<DataType>(i);
 
        // our local streamline diffusion coefficients
        Tiny::Vector<DataType, max_local_dofs> streamdiff_coeffs;
        streamdiff_coeffs.format();
 
        // our local delta for streamline diffusion
        DataType local_delta = DataType(0);
 
        // loop over all cells of the mesh
        for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
        {
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // gather our local convection dofs
          local_conv_dofs.format();
          gather_conv(local_conv_dofs, dof_mapping);
 
          // compute mesh width if necessary
          if(need_streamdiff)
          {
            // reset local delta
            local_delta = DataType(0);
 
            // evaluate trafo and space at barycentre
            trafo_eval(trafo_data, barycentre);
            space_eval(space_data, trafo_data);
 
            // compute velocity at barycentre
            mean_v.format();
            for(int i(0); i < num_loc_dofs; ++i)
              mean_v.axpy(space_data.phi[i].value, local_conv_dofs[i]);
 
            // compute norm of mean velocity
            const DataType local_norm_v = mean_v.norm_euclid();
 
            // do we have a non-zero velocity?
            if(local_norm_v > tol_eps)
            {
              // compute local mesh width w.r.t. mean velocity
              const DataType local_h = trafo_eval.width_directed(mean_v) * local_norm_v;
 
              // compute local Re_T
              const DataType local_re = (local_norm_v * local_h) / this->sd_nu;
 
              // compute local delta
              local_delta = this->sd_delta * (local_h / this->sd_v_norm) * (DataType(2)*local_re) / (DataType(1) + local_re);
            }
          }
 
          // format our local matrix and vector
          local_matrix.format();
 
          // loop over all quadrature points and integrate
          for(int point(0); point < cubature_rule.get_num_points(); ++point)
          {
            // compute trafo data
            trafo_eval(trafo_data, cubature_rule.get_point(point));
 
            // compute basis function data
            space_eval(space_data, trafo_data);
 
            // pre-compute cubature weight
            const DataType weight = trafo_data.jac_det * cubature_rule.get_weight(point);
 
            // evaluate convection function and its gradient (if required)
            if(need_conv || need_streamdiff)
            {
              loc_v.format();
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // update velocity value
                loc_v.axpy(space_data.phi[i].value, local_conv_dofs[i]);
              }
            }
            // always needed for material model
            loc_grad_v.format();
            for(int i(0); i < num_loc_dofs; ++i)
            {
              // update velocity gradient
              loc_grad_v.add_outer_product(local_conv_dofs[i], space_data.phi[i].grad);
            }
 
 
            // evaluate streamline diffusion coefficients
            if(need_streamdiff)
            {
              for(int i(0); i < num_loc_dofs; ++i)
              {
                streamdiff_coeffs[i] = Tiny::dot(loc_v, space_data.phi[i].grad);
              }
            }
 
            //calculate nu_loc
 
            strain_rate_tensor_2.set_transpose(loc_grad_v);
            strain_rate_tensor_2.axpy(DataType(1.0), loc_grad_v);
            DataType_ gamma_dot = Math::sqrt(0.5)*strain_rate_tensor_2.norm_frobenius();
            DataType_ nu_loc = mu_inf + (mu_0 - mu_inf) * Math::pow((DataType_(1.0) + Math::pow(lambda*gamma_dot, a)), (n-DataType(1.0)) / a);
 
 
            // assemble diffusion matrix?
            if(diff && !deformation)
            {
              // assemble gradient-tensor diffusion
 
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = nu_loc * weight * Tiny::dot(space_data.phi[i].grad, space_data.phi[j].grad);
 
                  // update local matrix
                  local_matrix[i][j].add_scalar_main_diag(value);
                }
              }
            }
            else if(diff && deformation)
            {
 
              // assemble deformation-tensor diffusion
 
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute inner product of grad(phi) and grad(psi)
                  const DataType value = nu_loc * weight * Tiny::dot(space_data.phi[j].grad, space_data.phi[i].grad);
 
                  // update local matrix
                  local_matrix[i][j].add_scalar_main_diag(value);
 
                  // add outer product of grad(phi) and grad(psi)
                  local_matrix[i][j].add_outer_product(space_data.phi[j].grad, space_data.phi[i].grad, nu_loc * weight);
                }
              }
            }
 
            if(frechet_diff)
            {
              XASSERT(deformation);
              const DataType_ mu_prime = (mu_0 - mu_inf) * ((n-DataType_(1.0)) / a) *
                Math::pow( DataType(1.0) + Math::pow(lambda*gamma_dot, a), ((n-DataType_(1.0) -a) / a)) *
                a * lambda * Math::pow(lambda*gamma_dot, a-DataType(1.0));
              const DataType fac =  mu_prime / (gamma_dot + reg_eps);
              //std::cout << fac ;
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                VectorValue du_grad_phi;
                du_grad_phi.set_mat_vec_mult(strain_rate_tensor_2, space_data.phi[i].grad);
 
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  VectorValue du_grad_psi;
                  du_grad_psi.set_mat_vec_mult(strain_rate_tensor_2, space_data.phi[j].grad);
                  // add outer product of grad(phi) and grad(psi)
                  local_matrix[i][j].add_outer_product(du_grad_phi, du_grad_psi, fac * weight);
                }
              }
            }
 
 
 
 
            // assemble convection?
            if(need_conv)
            {
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = beta * weight * space_data.phi[i].value * Tiny::dot(loc_v, space_data.phi[j].grad);
 
                  // update local matrix
                  local_matrix[i][j].add_scalar_main_diag(value);
                }
              }
            }
 
            // assemble convection Frechet?
            if(need_conv_frechet)
            {
 
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = frechet_beta * weight * space_data.phi[i].value * space_data.phi[j].value;
 
                  // update local matrix
                  local_matrix[i][j].axpy(value, loc_grad_v);
                }
              }
            }
 
            // assemble reaction?
            if(need_reac)
            {
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = theta * weight *  space_data.phi[i].value * space_data.phi[j].value;
 
                  // update local matrix
                  local_matrix[i][j].add_scalar_main_diag(value);
                }
              }
            }
 
            // assemble streamline diffusion?
            if(need_streamdiff && (local_delta > tol_eps))
            {
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = local_delta * weight * streamdiff_coeffs[i] * streamdiff_coeffs[j];
 
                  // update local matrix
                  local_matrix[i][j].add_scalar_main_diag(value);
                }
              }
            }
 
            // continue with next cubature point
          }
 
          // scatter into matrix
          scatter_matrix(local_matrix, dof_mapping, dof_mapping, scale);
 
          // finish dof mapping
          dof_mapping.finish();
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
        }
      }
 
      template<typename Space_, typename CubatureFactory_>
      void assemble_vector(
        LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_>& vector,
        const LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_>& convect,
        const LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_>& primal,
        const Space_& space,
        const CubatureFactory_& cubature_factory,
        const DataType_ scale = DataType_(1)
        ) const
      {
        // validate matrix and vector dimensions
        XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector size");
        XASSERTM(convect.size() == space.get_num_dofs(), "invalid vector size");
 
        typedef LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_> VectorType;
 
 
 
        // first of all, let's see what we have to assemble
        const bool need_conv = (Math::abs(beta) > DataType(0));
        //const bool need_conv_frechet = (frechet_beta > DataType(0));
        const bool need_reac = (Math::abs(theta) > DataType(0));
 
        // define our assembly traits
        typedef AsmTraits1<DataType_, Space_, TrafoTags::jac_det, SpaceTags::value|SpaceTags::grad> AsmTraits;
 
        // fetch our trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // create cubature rule
        typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
 
        // create vector-scatter-axpy (if needed)
        typename VectorType::ScatterAxpy scatter_vector(vector);
 
        // create convection gather-axpy
        typename VectorType::GatherAxpy gather_conv(convect);
 
        // create primal gather-axpy
        typename VectorType::GatherAxpy gather_prim(primal);
 
        // get maximum number of local dofs
        static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
        // create local vector data
        typedef Tiny::Vector<DataType, dim_> VectorValue;
        typedef Tiny::Vector<VectorValue, max_local_dofs> LocalVectorType;
        LocalVectorType local_vector;
 
        // local convection field dofs
        LocalVectorType local_conv_dofs;
 
        // local primal vector dofs
        LocalVectorType local_prim_dofs;
 
        // our local velocity value
        Tiny::Vector<DataType, dim_> loc_v;
 
        // our local velocity gradient
        Tiny::Matrix<DataType, dim_, dim_> loc_grad_v;
 
        Tiny::Matrix<DataType, dim_, dim_> strain_rate_tensor_2;
 
        loc_v.format();
        loc_grad_v.format();
 
        // loop over all cells of the mesh
        for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
        {
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // gather our local convection dofs
          local_conv_dofs.format();
          gather_conv(local_conv_dofs, dof_mapping);
 
          // gather our local primal dofs
          local_prim_dofs.format();
          gather_prim(local_prim_dofs, dof_mapping);
 
          // format our local vector
          local_vector.format();
 
          // loop over all quadrature points and integrate
          for(int point(0); point < cubature_rule.get_num_points(); ++point)
          {
            // compute trafo data
            trafo_eval(trafo_data, cubature_rule.get_point(point));
 
            // compute basis function data
            space_eval(space_data, trafo_data);
 
            // pre-compute cubature weight
            const DataType weight = trafo_data.jac_det * cubature_rule.get_weight(point);
 
            // evaluate convection function and its gradient (if required)
            if(need_conv)
            {
              loc_v.format();
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // update velocity value
                loc_v.axpy(space_data.phi[i].value, local_conv_dofs[i]);
              }
            }
            if(diff || frechet_diff)
            {
 
              loc_grad_v.format();
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // update velocity gradient
                loc_grad_v.add_outer_product(local_conv_dofs[i], space_data.phi[i].grad);
              }
 
            }
 
            loc_grad_v.format();
            for(int i(0); i < num_loc_dofs; ++i)
            {
              // update velocity gradient
              loc_grad_v.add_outer_product(local_conv_dofs[i], space_data.phi[i].grad);
            }
 
            //calculate nu_loc
            strain_rate_tensor_2.set_transpose(loc_grad_v);
            strain_rate_tensor_2.axpy(DataType(1.0), loc_grad_v);
            DataType_ gamma_dot = Math::sqrt(0.5)*strain_rate_tensor_2.norm_frobenius();
            DataType_ nu_loc = mu_inf + (mu_0 - mu_inf) * Math::pow((DataType_(1.0) + Math::pow(lambda*gamma_dot, a)), (n-DataType(1.0)) / a);
 
 
            // assemble diffusion matrix?
            if(diff && !deformation)
            {
              // assemble gradient-tensor diffusion
 
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = nu_loc * weight * Tiny::dot(space_data.phi[i].grad, space_data.phi[j].grad);
 
                  // update local vector
                  local_vector[i].axpy(value, local_prim_dofs[j]);
                }
              }
            }
            else if(diff && deformation)
            {
              // assemble deformation-tensor diffusion
 
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute inner product of grad(phi) and grad(psi)
                  const DataType value1 = nu_loc * weight * Tiny::dot(space_data.phi[i].grad, space_data.phi[j].grad);
 
                  // compute outer product of grad(phi) and grad(psi)
                  const DataType value2 = nu_loc * weight * Tiny::dot(local_prim_dofs[j], space_data.phi[i].grad);
 
                  // update local vector
                  local_vector[i].axpy(value1, local_prim_dofs[j]);
                  local_vector[i].axpy(value2, space_data.phi[j].grad);
                }
              }
            }
 
            // assemble convection?
            if(need_conv)
            {
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = beta * weight * space_data.phi[i].value * Tiny::dot(loc_v, space_data.phi[j].grad);
 
                  // update local vector
                  local_vector[i].axpy(value, local_prim_dofs[j]);
                }
              }
            }
 
            // assemble reaction?
            if(need_reac)
            {
              // test function loop
              for(int i(0); i < num_loc_dofs; ++i)
              {
                // trial function loop
                for(int j(0); j < num_loc_dofs; ++j)
                {
                  // compute scalar value
                  const DataType value = theta * weight *  space_data.phi[i].value * space_data.phi[j].value;
 
                  // update local vector
                  local_vector[i].axpy(value, local_prim_dofs[j]);
                }
              }
            }
 
            // continue with next cubature point
          }
 
          // scatter into vector
          scatter_vector(local_vector, dof_mapping, scale);
 
          // finish dof mapping
          dof_mapping.finish();
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
        }
      }
 
      void set_sd_v_norm(const LAFEM::DenseVectorBlocked<DataType_, IndexType_, dim_>& convect)
      {
        const auto* vals = convect.elements();
        DataType_ r = DataType(0);
        for(Index i(0); i < convect.size(); ++i)
          r = Math::max(r, vals[i].norm_euclid());
        this->sd_v_norm = r;
      }
 
      template<typename LocalVector_, typename Mirror_>
      void set_sd_v_norm(const Global::Vector<LocalVector_, Mirror_>& convect)
      {
        this->set_sd_v_norm(convect.local());
        const auto* gate = convect.get_gate();
        if(gate != nullptr)
          this->sd_v_norm = gate->max(this->sd_v_norm);
      }
    }; // class BurgersAssembler<...>
  } // namespace Assembly
} // namespace FEAT