naviersystem.c

一个用来实现偏微分方程中网格的计算库

/* $Id: naviersystem.C 2790 2008-04-13 16:43:55Z roystgnr $ */

/* The Next Great Finite Element Library. */
/* Copyright (C) 2003  Benjamin S. Kirk */

/* This library 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. */

/* This library 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 this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */

#include "getpot.h"

#include "naviersystem.h"

#include "boundary_info.h"
#include "fe_base.h"
#include "fe_interface.h"
#include "mesh.h"
#include "quadrature.h"


void NavierSystem::init_data ()
{
  const unsigned int dim = this->get_mesh().mesh_dimension();

  // Add the pressure variable "p". This will
  // be approximated with a first-order basis,
  // providing an LBB-stable pressure-velocity pair.
  // Add the velocity components "u" & "v".  They
  // will be approximated using second-order approximation.
  this->add_variable ("u", SECOND);
  u_var = 0;
  this->add_variable ("v", SECOND);
  v_var = 1;

  if (dim == 3)
    {
      this->add_variable ("w", SECOND);
      w_var = 2;
      p_var = 3;
    }
  else
    {
      w_var = u_var;
      p_var = 2;
    }

  this->add_variable ("p", FIRST);

  // Do the parent's initialization after variables are defined
  FEMSystem::init_data();

  // Tell the system to march velocity forward in time, but 
  // leave p as a constraint only
  this->time_evolving(u_var);
  this->time_evolving(v_var);
  if (dim == 3)
    this->time_evolving(w_var);

  // Get references to the finite elements we need
  fe_velocity = element_fe[this->variable_type(u_var)];
  fe_pressure = element_fe[this->variable_type(p_var)];
  fe_side_vel = side_fe[this->variable_type(u_var)];

  // To enable FE optimizations, we should prerequest all the data
  // we will need to build the linear system.
  fe_velocity->get_JxW();
  fe_velocity->get_phi();
  fe_velocity->get_dphi();
  fe_velocity->get_xyz();
  
  fe_pressure->get_phi();

  fe_side_vel->get_JxW();
  fe_side_vel->get_phi();

  // Check the input file for Reynolds number
  GetPot infile("navier.in");
  Reynolds = infile("Reynolds", 1.);
  application = infile("application", 0);

  // Useful debugging options
  // Set verify_analytic_jacobians to 1e-6 to use
  this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0);
  this->print_jacobians = infile("print_jacobians", false);
  this->print_element_jacobians = infile("print_element_jacobians", false);
}


bool NavierSystem::element_time_derivative (bool request_jacobian)
{
  // First we get some references to cell-specific data that
  // will be used to assemble the linear system.

  // Element Jacobian * quadrature weights for interior integration
  const std::vector<Real> &JxW = fe_velocity->get_JxW();

  // The velocity shape functions at interior quadrature points.
  const std::vector<std::vector<Real> >& phi = fe_velocity->get_phi();

  // The velocity shape function gradients at interior
  // quadrature points.
  const std::vector<std::vector<RealGradient> >& dphi =
    fe_velocity->get_dphi();

  // The pressure shape functions at interior
  // quadrature points.
  const std::vector<std::vector<Real> >& psi = fe_pressure->get_phi();

  // Physical location of the quadrature points
  const std::vector<Point>& qpoint = fe_velocity->get_xyz();
 
  // The number of local degrees of freedom in each variable
  const unsigned int n_p_dofs = dof_indices_var[p_var].size();
  const unsigned int n_u_dofs = dof_indices_var[u_var].size(); 
  libmesh_assert (n_u_dofs == dof_indices_var[v_var].size()); 

  // The subvectors and submatrices we need to fill:
  const unsigned int dim = this->get_mesh().mesh_dimension();
  DenseSubMatrix<Number> &Kuu = *elem_subjacobians[u_var][u_var];
  DenseSubMatrix<Number> &Kvv = *elem_subjacobians[v_var][v_var];
  DenseSubMatrix<Number> &Kww = *elem_subjacobians[w_var][w_var];
  DenseSubMatrix<Number> &Kuv = *elem_subjacobians[u_var][v_var];
  DenseSubMatrix<Number> &Kuw = *elem_subjacobians[u_var][w_var];
  DenseSubMatrix<Number> &Kvu = *elem_subjacobians[v_var][u_var];
  DenseSubMatrix<Number> &Kvw = *elem_subjacobians[v_var][w_var];
  DenseSubMatrix<Number> &Kwu = *elem_subjacobians[w_var][u_var];
  DenseSubMatrix<Number> &Kwv = *elem_subjacobians[w_var][v_var];
  DenseSubMatrix<Number> &Kup = *elem_subjacobians[u_var][p_var];
  DenseSubMatrix<Number> &Kvp = *elem_subjacobians[v_var][p_var];
  DenseSubMatrix<Number> &Kwp = *elem_subjacobians[w_var][p_var];
  DenseSubVector<Number> &Fu = *elem_subresiduals[u_var];
  DenseSubVector<Number> &Fv = *elem_subresiduals[v_var];
  DenseSubVector<Number> &Fw = *elem_subresiduals[w_var];
      
  // Now we will build the element Jacobian and residual.
  // Constructing the residual requires the solution and its
  // gradient from the previous timestep.  This must be
  // calculated at each quadrature point by summing the
  // solution degree-of-freedom values by the appropriate
  // weight functions.
  unsigned int n_qpoints = element_qrule->n_points();

  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      // Compute the solution & its gradient at the old Newton iterate
      Number p = interior_value(p_var, qp),
             u = interior_value(u_var, qp),
             v = interior_value(v_var, qp),
             w = interior_value(w_var, qp);
      Gradient grad_u = interior_gradient(u_var, qp),
               grad_v = interior_gradient(v_var, qp),
               grad_w = interior_gradient(w_var, qp);

      // Definitions for convenience.  It is sometimes simpler to do a
      // dot product if you have the full vector at your disposal.
      NumberVectorValue U     (u,     v);
      if (dim == 3)
        U(2) = w;
      const Number  u_x = grad_u(0);
      const Number  u_y = grad_u(1);
      const Number  u_z = (dim == 3)?grad_u(2):0;
      const Number  v_x = grad_v(0);
      const Number  v_y = grad_v(1);
      const Number  v_z = (dim == 3)?grad_v(2):0;
      const Number  w_x = (dim == 3)?grad_w(0):0;
      const Number  w_y = (dim == 3)?grad_w(1):0;
      const Number  w_z = (dim == 3)?grad_w(2):0;

      // Value of the forcing function at this quadrature point
      Point f = this->forcing(qpoint[qp]);

      // First, an i-loop over the velocity degrees of freedom.
      // We know that n_u_dofs == n_v_dofs so we can compute contributions
      // for both at the same time.
      for (unsigned int i=0; i != n_u_dofs; i++)
        {
          Fu(i) += JxW[qp] *
                   (-Reynolds*(U*grad_u)*phi[i][qp] + // convection term
                    p*dphi[i][qp](0) -                // pressure term
		    (grad_u*dphi[i][qp]) +            // diffusion term
		    f(0)*phi[i][qp]                   // forcing function
		    );

            
          Fv(i) += JxW[qp] *
                   (-Reynolds*(U*grad_v)*phi[i][qp] + // convection term
                    p*dphi[i][qp](1) -                // pressure term
		    (grad_v*dphi[i][qp]) +            // diffusion term
		    f(1)*phi[i][qp]                   // forcing function
		    );


          if (dim == 3)
          Fw(i) += JxW[qp] *
                   (-Reynolds*(U*grad_w)*phi[i][qp] + // convection term
                    p*dphi[i][qp](2) -                // pressure term
		    (grad_w*dphi[i][qp]) +            // diffusion term
		    f(2)*phi[i][qp]                   // forcing function
		    );


          // Note that the Fp block is identically zero unless we are using
          // some kind of artificial compressibility scheme...

          // Matrix contributions for the uu and vv couplings.
          if (request_jacobian)
            for (unsigned int j=0; j != n_u_dofs; j++)
              {
                Kuu(i,j) += JxW[qp] *
 /* convection term */      (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
 /* diffusion term  */       (dphi[i][qp]*dphi[j][qp]) -
 /* Newton term     */       Reynolds*u_x*phi[i][qp]*phi[j][qp]);

                Kuv(i,j) += JxW[qp] *
 /* Newton term     */      -Reynolds*u_y*phi[i][qp]*phi[j][qp];

                Kvv(i,j) += JxW[qp] *
 /* convection term */      (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
 /* diffusion term  */       (dphi[i][qp]*dphi[j][qp]) -
 /* Newton term     */       Reynolds*v_y*phi[i][qp]*phi[j][qp]);

                Kvu(i,j) += JxW[qp] * 
 /* Newton term     */      -Reynolds*v_x*phi[i][qp]*phi[j][qp];
                if (dim == 3)
                  {
                    Kww(i,j) += JxW[qp] *
 /* convection term */          (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
 /* diffusion term  */           (dphi[i][qp]*dphi[j][qp]) -
 /* Newton term     */           Reynolds*w_z*phi[i][qp]*phi[j][qp]);
                    Kuw(i,j) += JxW[qp] *
 /* Newton term     */      -Reynolds*u_z*phi[i][qp]*phi[j][qp];
                    Kvw(i,j) += JxW[qp] *
 /* Newton term     */      -Reynolds*v_z*phi[i][qp]*phi[j][qp];
                    Kwu(i,j) += JxW[qp] *
 /* Newton term     */      -Reynolds*w_x*phi[i][qp]*phi[j][qp];
                    Kwv(i,j) += JxW[qp] *
 /* Newton term     */      -Reynolds*w_y*phi[i][qp]*phi[j][qp];
                  }
              }

          // Matrix contributions for the up and vp couplings.
          if (request_jacobian)
            for (unsigned int j=0; j != n_p_dofs; j++)
              {
                Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
                Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
                if (dim == 3)
                  Kwp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](2);
              }
        }
    } // end of the quadrature point qp-loop
  
  return request_jacobian;
}



bool NavierSystem::element_constraint (bool request_jacobian)
{
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.

  // Element Jacobian * quadrature weight for interior integration
  const std::vector<Real> &JxW = fe_velocity->get_JxW();

  // The velocity shape function gradients at interior
  // quadrature points.
  const std::vector<std::vector<RealGradient> >& dphi =
    fe_velocity->get_dphi();

  // The pressure shape functions at interior
  // quadrature points.
  const std::vector<std::vector<Real> >& psi = fe_pressure->get_phi();

  // The number of local degrees of freedom in each variable
  const unsigned int n_u_dofs = dof_indices_var[u_var].size();
  const unsigned int n_p_dofs = dof_indices_var[p_var].size();

  // The subvectors and submatrices we need to fill:
  const unsigned int dim = this->get_mesh().mesh_dimension();