Dear Libmesher,
I tried to write an incompressible Naiver-Stokes solver based Libmesh.
Following features includes:
coupled pressure-velocity method
TransientNonlinearImplicitSystem used with
compute_jacobian(..) and compute_residual(..) run separatedly
Lid-driven cavity case for initial tests
BDF1 time marching
QUAD4 element,
vel_order=2, pres-order=1, both LAGRANGE elements
Pressure pinned at node_id=0
I tried DirichletBoundary class to impose the velocity BC, the code runs OK.
Now the problem is :
I coded myself the common penalty method for term
\int_{\Gamma_D}dt*\gamma/h*u*v,
and applied the forms in boundary integral part of jacobian function
Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] *phi_face_[i][qp_face]
and
Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(u-u_bc)
at counterparts of residual function
the results show the boundary velocity not fully applied ,esp. at the side
walls near
the two top corners. Actually there are nonzero u normal to the local side wall.
So I wonder anybody here will kindly show me some hints.
Furthermore, I tried applied the weak Dirichlet conditions by add futher
boundary integrals like
dt\int_{\partial\Omega} (-\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} +
dt\int_{\partial\Omega} (-\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u}
at left and
dt\int_{\partial\Omega} (-\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u_bc}
The case is even worse.
I also attached the code, and I am looking forward to any suggestions. Many
thanks.
Zhenyu
/* The libMesh 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 */
// <h1>Systems Example 2 - Unsteady Nonlinear Navier-Stokes</h1>
//
// This example shows how a simple, unsteady, nonlinear system of equations
// can be solved in parallel. The system of equations are the familiar
// Navier-Stokes equations for low-speed incompressible fluid flow. This
// example introduces the concept of the inner nonlinear loop for each
// timestep, and requires a good deal of linear algebra number-crunching
// at each step. If you have a ExodusII viewer such as ParaView installed,
// the script movie.sh in this directory will also take appropriate screen
// shots of each of the solution files in the time sequence. These rgb files
// can then be animated with the "animate" utility of ImageMagick if it is
// installed on your system. On a PIII 1GHz machine in debug mode, this
// example takes a little over a minute to run. If you would like to see
// a more detailed time history, or compute more timesteps, that is certainly
// possible by changing the n_timesteps and dt variables below.
//c.f.
//miscellaneous_ex5.C
//
// Required libraries:
// GetOpt: parameters input
// GINAC: symbol processing
// Petsc/Trilinos: Equation solver
// Exodus/VTK: postprocessing
// Paraview:
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <sstream>
#include <cmath>
// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/perf_log.h"
#include "libmesh/boundary_info.h"
#include "libmesh/utility.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
// For systems of equations the \p DenseSubMatrix
// and \p DenseSubVector provide convenient ways for
// assembling the element matrix and vector on a
// component-by-component basis.
#include "libmesh/dense_submatrix.h"
#include "libmesh/dense_subvector.h"
// The definition of a geometric element
#include "libmesh/elem.h"
#include "libmesh/nonlinear_solver.h"
#include "libmesh/nonlinear_implicit_system.h"
#ifdef LIBMESH_HAVE_PETSC
#include <petsc.h>
#endif
// Bring in everything from the libMesh namespace
using namespace libMesh;
// A reference to our equation system
EquationSystems *_equation_system = NULL;
struct InputParameters
{
Real Reynolds;
Real dt;
unsigned int nstep;
Real penalty; // for inner panalty method of Dirichlet weak BC
Order vel_order;
Order pres_order;
bool pin_pressure;
unsigned int pin_pressure_node;
Real penalty_pin; //penalty of pin_pressrure
Real pref;
};
InputParameters input;
void compute_residual (const NumericVector<Number>& soln,
NumericVector<Number>& residual,
NonlinearImplicitSystem& /*sys*/);
void compute_jacobian (const NumericVector<Number>& soln,
SparseMatrix<Number>& jacobian,
NonlinearImplicitSystem& /*sys*/);
void ReadInputParameters(std::string& filename, InputParameters& input)
{
GetPot input_file(filename);
//Read in parameters from the input file
input.Reynolds=input_file("Reynolds", 100.);
input.dt=input_file("delta_t", 0.01);
input.pref=input_file("pressure_reference", 0.);
input.nstep=input_file("n_timestep", 25);
input.penalty=input_file("penalty", 10.);
input.vel_order=static_cast<Order>(input_file("vel_order",2));
input.pres_order=static_cast<Order>(input_file("pres_order",1));
input.pin_pressure=input_file("pin_pressure", true);
input.penalty_pin=input_file("penalty_pin", 1.e10);
input.pin_pressure_node=input_file("pin_pressure_node", 0);
input.pref=input_file("pref", 0.0);
}
// The main program.
int main (int argc, char** argv)
{
// Initialize libMesh.
LibMeshInit init (argc, argv);
//------------------------------------------------------------------
// Create a GetPot object to parse the command line
GetPot command_line (argc, argv);
// Read number of refinements
int nr = 2;
if ( command_line.search(1, "-r") )
nr = command_line.next(nr);
// Read FE order from command line
std::string order = "SECOND";
if ( command_line.search(2, "-Order", "-o") )
order = command_line.next(order);
// Read FE Family from command line
std::string family = "LAGRANGE";
if ( command_line.search(2, "-FEFamily", "-f") )
family = command_line.next(family);
// Cannot use dicontinuous basis.
if ((family == "MONOMIAL") || (family == "XYZ"))
libmesh_error_msg("This example requires a C^0 (or higher) FE basis.");
if ( command_line.search(1, "-pre") )
{
#ifdef LIBMESH_HAVE_PETSC
//Use the jacobian for preconditioning.
PetscOptionsSetValue("-snes_mf_operator",PETSC_NULL);
#else
std::cerr<<"Must be using PetsC to use jacobian based
preconditioning"<<std::endl;
//returning zero so that "make run" won't fail if we ever enable this
capability there.
return 0;
#endif //LIBMESH_HAVE_PETSC
}
std::cout<<"OK 0-1"<<std::endl;
std::string paramfile="laminar2dnl.param";
// std::cout<<"paramfile: "<<paramfile<<std::endl;
// if ( command_line.search(3, "-param", "-par") ){
// paramfile = command_line.next(paramfile);
// std::cout<<"paramfile: "<<paramfile<<std::endl;
// }
std::cout<<"OK 0-1-1"<<std::endl;
ReadInputParameters(paramfile,input);
std::cout<<"Re: "<<input.Reynolds<<"\n";
std::cout<<"OK 0-2"<<std::endl;
// Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
// std::cout<<"OK 0-3"<<std::endl;
// This example NaNs with the Eigen sparse linear solvers and
// Trilinos solvers, but should work OK with either PETSc or
// Laspack.
libmesh_example_requires(libMesh::default_solver_package() !=
EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");
libmesh_example_requires(libMesh::default_solver_package() !=
TRILINOS_SOLVERS, "--enable-petsc or --enable-laspack");
// std::cout<<"OK 1"<<std::endl;
//------------------------------------------------------------------
// Create a mesh, with dimension to be overridden later, distributed
// across the default MPI communicator.
Mesh mesh(init.comm());
// Use the MeshTools::Generation mesh generator to create a uniform
// 2D grid on the square [-1,1]^2. We instruct the mesh generator
// to build a mesh of 8x8 \p Quad9 elements in 2D. Building these
// higher-order elements allows us to use higher-order
// approximation, as in example 3.
MeshTools::Generation::build_square (mesh,
20, 20,
0., 1.,
0., 1.,
QUAD4);
mesh.all_second_order();
// Print information about the mesh to the screen.
mesh.print_info();
EquationSystems equation_systems (mesh);
_equation_system = &equation_systems;
// Declare the system and its variables.
// Creates a transient system named "Navier-Stokes"
TransientNonlinearImplicitSystem & system =
equation_systems.add_system<TransientNonlinearImplicitSystem>
("Navier-Stokes");
// Here we specify the tolerance for the nonlinear solver and
// the maximum of nonlinear iterations.
equation_systems.parameters.set<Real> ("nonlinear solver
tolerance") = 1.e-6;
equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum
iterations") = 50;
// Add the variables "u" & "v" to "Navier-Stokes". They
// will be approximated using second-order approximation.
//~ system.add_variable ("u", SECOND);
//~ system.add_variable ("v", SECOND);
Order vel_order=Utility::string_to_enum<Order>(order);
system.add_variable
("u",vel_order,Utility::string_to_enum<FEFamily>(family));
system.add_variable ("v",
Utility::string_to_enum<Order>(order),Utility::string_to_enum<FEFamily>(family));
// Add the variable "p" to "Navier-Stokes". This will
// be approximated with a first-order basis,
// providing an LBB-stable pressure-velocity pair.
//~ Order pres_order=static_cast<Order>(vel_order-1);
Order pres_order=Order(vel_order-1);
std::cout<<"pres_order:"<<pres_order<<"\n";
//~ pres_order;
system.add_variable ("p", pres_order);
//~ system.add_variable ("p", Utility::string_to_enum<Order>(order)-1);
// Give the system a pointer to the matrix assembly
// function.
//~ system.attach_assemble_function (assemble_stokes);
// Give the system a pointer to the functions that update
// the residual and Jacobian.
system.nonlinear_solver->residual = compute_residual;
system.nonlinear_solver->jacobian = compute_jacobian;
// Initialize the data structures for the equation system.
equation_systems.init ();
// Prints information about the system to the screen.
equation_systems.print_info();
//------------------------------------------------------------------
equation_systems.parameters.set<Real> ("viscosity") = 1./input.Reynolds;
equation_systems.parameters.set<Real> ("pref") = input.pref;
equation_systems.parameters.set<Real> ("penalty") = input.penalty;
std::cout<<"penalty: "<<input.penalty<<"\n";
equation_systems.parameters.set<Real> ("penalty_pin") = input.penalty_pin;
equation_systems.parameters.set<bool> ("pin_pressure") =
input.pin_pressure;
equation_systems.parameters.set<unsigned int> ("pin_pressure_node") =
input.pin_pressure_node;
// Create a performance-logging object for this example
PerfLog perf_log("Systems Example 2");
// Get a reference to the Stokes system to use later.
TransientNonlinearImplicitSystem& navier_stokes_system =
equation_systems.get_system<TransientNonlinearImplicitSystem>("Navier-Stokes");
// Tell the system of equations what the timestep is by using
// the set_parameter function. The matrix assembly routine can
// then reference this parameter.
equation_systems.parameters.set<Real> ("dt") = input.dt;
const unsigned int n_timesteps = input.nstep;
navier_stokes_system.time = 0.0;
// The number of steps and the stopping criterion are also required
// for the nonlinear iterations.
//~ const unsigned int n_nonlinear_steps = 20;
// const Real nonlinear_tolerance = 1.e-3;
// We also set a standard linear solver flag in the EquationSystems object
// which controls the maxiumum number of linear solver iterations allowed.
//~ equation_systems.parameters.set<unsigned int>("linear solver maximum
iterations") = 250;
// The first thing to do is to get a copy of the solution at
// the current nonlinear iteration. This value will be used to
// determine if we can exit the nonlinear loop.
AutoPtr<NumericVector<Number> >
last_nonlinear_soln (navier_stokes_system.solution->clone());
for (unsigned int t_step=0; t_step<n_timesteps; ++t_step)
{
last_nonlinear_soln->zero();
last_nonlinear_soln->add(*navier_stokes_system.solution);
navier_stokes_system.time += input.dt;
// A pretty update message
std::cout << "\n\n*** Solving time step " << t_step <<
", time = " << navier_stokes_system.time <<
" ***" << std::endl;
// Now we need to update the solution vector from the
// previous time step.
*navier_stokes_system.old_local_solution =
*navier_stokes_system.current_local_solution;
// At the beginning of each solve, reset the linear solver tolerance
// to a "reasonable" starting value.
//~ const Real initial_linear_solver_tol = 1.e-6;
//~ equation_systems.parameters.set<Real> ("linear solver tolerance") =
initial_linear_solver_tol;
perf_log.push("nonlinear solve");
equation_systems.get_system("Navier-Stokes").solve();
perf_log.pop("nonlinear solve");
last_nonlinear_soln->add (-1., *navier_stokes_system.solution);
last_nonlinear_soln->close();
const Real norm_delta = last_nonlinear_soln->l2_norm();
std::cout<<" Nonlinear convergence: ||u - u_old|| = "
<< norm_delta
<< std::endl;
#ifdef LIBMESH_HAVE_EXODUS_API
// Write out every nth timestep to file.
const unsigned int write_interval = 2;
if ((t_step+1)%write_interval == 0)
{
std::ostringstream file_name;
// We write the file in the ExodusII format.
file_name << "out_"
<< "Re_"<<int(input.Reynolds)<<"_"
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step + 1
<< ".e";
ExodusII_IO(mesh).write_equation_systems (file_name.str(),
equation_systems);
}
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
} // end timestep loop.
return 0;
}
// This function computes the Jacobian K(x)
void compute_jacobian (const NumericVector<Number>& soln,
SparseMatrix<Number>& jacobian,
NonlinearImplicitSystem& /*sys*/)
{
// Get a reference to the equation system.
EquationSystems &es = *_equation_system;
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
// Get a reference to the NonlinearImplicitSystem we are solving
TransientNonlinearImplicitSystem& system =
es.get_system<TransientNonlinearImplicitSystem>("Navier-Stokes");
// Numeric ids corresponding to each variable in the system
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int p_var = system.variable_number ("p");
FEType fe_vel_type = system.variable_type(u_var);
FEType fe_pres_type = system.variable_type(p_var);
// Finite Element objects of for vel and pres
AutoPtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
// A reference to the \p DofMap object for this system. The \p DofMap
// object handles the index translation from node and element numbers
// to degree of freedom numbers. We will talk more about the \p DofMap
// in future examples.
const DofMap& dof_map = system.get_dof_map();
QGauss qrule (dim, fe_vel_type.default_quadrature_order());
fe_vel->attach_quadrature_rule (&qrule);
fe_pres->attach_quadrature_rule (&qrule);
const std::vector<Real>& JxW = fe_vel->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe_vel->get_phi();
// The element shape function gradients for the velocity
// variables evaluated at the quadrature points.
const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();
// The element shape functions for the pressure variable
// evaluated at the quadrature points.
const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();
#if 1
// Declare a special finite element object for
// boundary integration.
AutoPtr<FEBase> fe_vel_face (FEBase::build(dim, fe_vel_type));
AutoPtr<FEBase> fe_pres_face (FEBase::build(dim, fe_pres_type));
// Boundary integration requires one quadraure rule,
// with dimensionality one less than the dimensionality
// of the element.
QGauss qface(dim-1, fe_vel_type.default_quadrature_order());
QGauss qface_pres(dim-1, fe_pres_type.default_quadrature_order());
// Tell the finte element object to use our
// quadrature rule.
fe_vel_face->attach_quadrature_rule (&qface);
fe_pres_face->attach_quadrature_rule (&qface_pres);
//c.f. miscellaneous_ex5
const std::vector<std::vector<Real> >& phi_face = fe_vel_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face =
fe_vel_face->get_dphi();
const std::vector<Real>& JxW_face = fe_vel_face->get_JxW();
const std::vector<Point>& qface_normals = fe_vel_face->get_normals();
const std::vector<Point>& qface_points = fe_vel_face->get_xyz();
//pressure terms
const std::vector<std::vector<Real> >& psi_face = fe_pres_face->get_phi();
const std::vector<std::vector<RealGradient> >& dpsi_face =
fe_pres_face->get_dphi();
const std::vector<Real>& JxW_face_pres = fe_pres_face->get_JxW();
const std::vector<Point>& qface_normals_pres = fe_pres_face->get_normals();
const std::vector<Point>& qface_points_pres = fe_pres_face->get_xyz();
#endif
// Define data structures to contain the Jacobian element matrix.
DenseMatrix<Number> Ke;
DenseSubMatrix<Number>
Kuu(Ke), Kuv(Ke), Kup(Ke),
Kvu(Ke), Kvv(Ke), Kvp(Ke),
Kpu(Ke), Kpv(Ke), Kpp(Ke);
std::vector<dof_id_type> dof_indices;
std::vector<dof_id_type> dof_indices_u;
std::vector<dof_id_type> dof_indices_v;
std::vector<dof_id_type> dof_indices_p;
// ???
const Real dt = es.parameters.get<Real>("dt");
const Real visc = es.parameters.get<Real>("viscosity");
const Real pref = es.parameters.get<Real>("pref");
const Real penalty = es.parameters.get<Real>("penalty");
const Real penalty_pin = es.parameters.get<Real>("penalty_pin");
const bool pin_pressure=es.parameters.get<bool>("pin_pressure");
const unsigned int pin_pressure_node=es.parameters.get<unsigned
int>("pin_pressure_node");
//~ const Real theta = 1.;
// Now we will loop over all the active elements in the mesh which
// are local to this processor.
// We will compute the element Jacobian contribution.
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
dof_map.dof_indices (elem, dof_indices_u, u_var);
dof_map.dof_indices (elem, dof_indices_v, v_var);
dof_map.dof_indices (elem, dof_indices_p, p_var);
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_u_dofs = dof_indices_u.size();
const unsigned int n_v_dofs = dof_indices_v.size();
const unsigned int n_p_dofs = dof_indices_p.size();
// Get the degree of freedom indices for the
// current element. These define where in the global
// matrix and right-hand-side this element will
// contribute to.
//~ dof_map.dof_indices (elem, dof_indices);
// Compute the element-specific data for the current
// element. This involves computing the location of the
// quadrature points (q_point) and the shape functions
// (phi, dphi) for the current element.
//~ fe->reinit (elem);
fe_vel->reinit (elem);
fe_pres->reinit (elem);
// Zero the element Jacobian before
// summing them. We use the resize member here because
// the number of degrees of freedom might have changed from
// the last element. Note that this will be the case if the
// element type is different (i.e. the last element was a
// triangle, now we are on a quadrilateral).
Ke.resize (n_dofs, n_dofs);
// Reposition the submatrices... The idea is this:
//
// - - - -
// | Kuu Kuv Kup | | Fu |
// Ke = | Kvu Kvv Kvp |; Fe = | Fv |
// | Kpu Kpv Kpp | | Fp |
// - - - -
//
// The \p DenseSubMatrix.repostition () member takes the
// (row_offset, column_offset, row_size, column_size).
//
// Similarly, the \p DenseSubVector.reposition () member
// takes the (row_offset, row_size)
Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
Number u,v;
Number temp;
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
u=0.;
v=0.;
Gradient grad_u,grad_v;
//~ for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int l=0; l<n_u_dofs; l++)
{
//~ grad_u += dphi[i][qp]*soln(dof_indices[i]);
//~ u += phi[l][qp]*system.current_solution(dof_indices_u[l]);
//~ v += phi[l][qp]*system.current_solution(dof_indices_v[l]);
// From the previous Newton iterate:
u += phi[l][qp]*soln(dof_indices_u[l]);
v += phi[l][qp]*soln(dof_indices_v[l]);
grad_u.add_scaled (dphi[l][qp],soln (dof_indices_u[l]));
grad_v.add_scaled (dphi[l][qp],soln (dof_indices_v[l]));
}
//~ const Number K = 1./std::sqrt(1. + grad_u*grad_u);
const NumberVectorValue U(u,v);
const Number u_x = grad_u(0);
const Number u_y = grad_u(1);
const Number v_x = grad_v(0);
const Number v_y = grad_v(1);
//~ for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int i=0; i<n_u_dofs; i++)
{
for (unsigned int j=0; j<n_u_dofs; j++)
{
Kuu(i,j) += JxW[qp]*(phi[j][qp]*phi[i][qp]
// mass matrix term
+dt*(visc*dphi[j][qp]*dphi[i][qp])
// diffusion stiffness term
+dt*(U*dphi[j][qp])*phi[i][qp] //
convection term
+dt*u_x*phi[j][qp]*phi[i][qp] //
Newton term
);
Kuv(i,j) += JxW[qp]*dt*u_y*phi[j][qp]*phi[i][qp]; //
Newton term
Kvu(i,j) += JxW[qp]*dt*v_x*phi[i][qp]*phi[j][qp]; //
Newton term
Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
// mass matrix term
dt*(visc*dphi[i][qp]*dphi[j][qp]) +
// diffusion term
dt*(U*dphi[j][qp])*phi[i][qp] + //
convection term
dt*v_y*phi[i][qp]*phi[j][qp]); //
Newton term
} // j loop of vel
// loop for Kup , Kvp
for (unsigned int j=0; j<n_p_dofs; j++)
{
Kup(i,j) += JxW[qp]*(-dt*psi[j][qp]*dphi[i][qp](0)); //
dphi[i][qp](0): \frac{\partial v_1}{\partial x}
Kvp(i,j) += JxW[qp]*(-dt*psi[j][qp]*dphi[i][qp](1)); //
dphi[i][qp](1): \frac{\partial v_2}{\partial y}
} // j loop of pres
} // i loop
// loop for Kpu, Kpv
for (unsigned int i=0; i<n_p_dofs; i++)
{
for (unsigned int j=0; j<n_u_dofs; j++)
{
temp=-JxW[qp]*dt*psi[i][qp];
Kpu(i,j) += temp*dphi[j][qp](0);
Kpv(i,j) += temp*dphi[j][qp](1);
}
}
} // qp loop of qrule
for (unsigned int side=0; side<elem->n_sides(); side++)
{
if(elem->neighbor(side) == NULL) //boundary
{
boundary_id_type
bid=mesh.get_boundary_info().boundary_id(elem,side);
AutoPtr<Elem> _side (elem->build_side(side));
fe_vel_face->reinit(elem, side);
fe_pres_face->reinit(elem, side);
const unsigned int elem_b_order = static_cast<unsigned int>
(fe_vel_face->get_order());
const double h_elem = elem->volume()/_side->volume() *
1./pow(elem_b_order, 2.);
// only for lid-driven cavity , apply Diri BC
if (bid<4 && bid>=0)
{
// Boundary ids are set internally by
// build_square().
// 0=bottom // 1=right // 2=top // 3=left
// Set u = 1 on the top boundary, 0 everywhere else
const Real u_value
=(mesh.get_boundary_info().has_boundary_id(elem,side,2))
? 1. : 0.;
const Real v_value = 0.;
Number temp;
// ******** Dirichlet BC by penalty
***************************
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
//~ p=0.;
//~ for (unsigned int l=0; l<n_p_dofs; l++)
//~ {
//~ p+=psi[l][qp]*soln(dof_indices_p[l]);
//~ }
//~ Number bc_value = exact_solution(qface_points[qp],
es.parameters,"null","void");
for (unsigned int i=0; i<phi_face.size(); i++)
{
for (unsigned int j=0; j<phi_face.size(); j++)
{
//~ (-idt(p)*N()+mu*gradt(u)*N())
//~ tau_n_t=
//~ (-id(p)*N()+mu*grad(v)*N())
//~ tau_n=-psi_face[i][qp]*
Kuu(i,j) += JxW_face[qp] * dt* penalty/h_elem *
phi_face[i][qp] * phi_face[j][qp]; // stability , term E_u
//Kuu(i,j) -= JxW_face[qp] * dt*
visc*(//phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) // term A_u,
//??? be zero here since as Dirichlet BC ?
//+
// phi_face[j][qp] *
(dphi_face[i][qp]*qface_normals[qp])); // term C_u, both for consistency
Kvv(i,j) += JxW_face[qp] * dt* penalty/h_elem *
phi_face[i][qp] * phi_face[j][qp]; // stability, E_v
//Kvv(i,j) -= JxW_face[qp] * dt*
visc*(//phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) // term A_v
//+
// phi_face[j][qp] *
(dphi_face[i][qp]*qface_normals[qp])); // term C_v, both for consistency
} // j loop of vel
// RHS contribution
//~ Fe(i) += JxW_face[qp] * bc_value *
penalty/h_elem * phi_face[i][qp]; // stability
//~ Fe(i) -= JxW_face[qp] * dphi_face[i][qp] *
(bc_value*qface_normals[qp]); // consistency
#if 0 //Dirichlet BC this may be zero
for (unsigned int j=0; j<psi_face.size(); j++)
{
temp=JxW_face[qp] * dt*
psi_face[j][qp]*phi_face[i][qp];
Kup(i,j)+= temp*qface_normals[qp](0); // B_u
Kvp(i,j)+= temp*qface_normals[qp](1); // B_v
} //j loop of pres
#endif
} // i loop of vel
} //qface.qp loop
#if 0 // D_u,D_v terms
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
for (unsigned int i=0; i<psi_face.size(); i++)
{
for (unsigned int j=0; j<phi_face.size(); j++)
{
temp=JxW_face_pres[qp] * dt*
psi_face[i][qp]*phi_face[j][qp];
Kpu(i,j)+=temp*qface_normals_pres[qp](0); // D_u
Kpv(i,j)+=temp*qface_normals_pres[qp](1); // D_u
} // i loop for vel
} // j loop for pres
} // qface_pres loop
#endif
} // if (bid<4 && bid>=0)
} //if(elem->neighbor(side) == NULL)
else // inner nodes , face integration for DG procedure, see in
miscellaneous_ex5.C
{
} // inner nodes
} // side loop
if(pin_pressure)
{
// const unsigned int pressure_node = 0;
for (unsigned int c=0; c<elem->n_nodes(); c++)
{
if (elem->node(c) == pin_pressure_node)
{
Kpp(c,c) += penalty_pin;
}
}
} // pin_pressure
#if 0
{
// old version of
//beginning of boudary conditions
// The penalty value. \f$ \frac{1}{\epsilon} \f$
const Real penalty = 1.e10;
// The following loops over the sides of the element.
// If the element has no neighbor on a side then that
// side MUST live on a boundary of the domain.
for (unsigned int s=0; s<elem->n_sides(); s++)
{
if (elem->neighbor(s) == NULL)
{
AutoPtr<Elem> side (elem->build_side(s));
// Loop over the nodes on the side.
for (unsigned int ns=0; ns<side->n_nodes(); ns++)
{
// Boundary ids are set internally by
// build_square().
// 0=bottom // 1=right // 2=top // 3=left
// Set u = 1 on the top boundary, 0 everywhere else
// const Real u_value =
//
(mesh.get_boundary_info().has_boundary_id(elem,s,2))
// ? 1. : 0.;
// Set v = 0 everywhere
// const Real v_value = 0.;
// Find the node on the element matching this node on
// the side. That defined where in the element matrix
// the boundary condition will be applied.
for (unsigned int n=0; n<elem->n_nodes(); n++)
if (elem->node(n) == side->node(ns))
{
// Matrix contribution.
Kuu(n,n) += penalty;
Kvv(n,n) += penalty;
//~ // Right-hand-side contribution.
//~ Fu(n) += penalty*u_value;
//~ Fv(n) += penalty*v_value;
}
} // end face node loop
} // end if (elem->neighbor(side) == NULL)
} // side loop
// Pin the pressure to zero at global node number "pressure_node".
// This effectively removes the non-trivial null space of constant
// pressure solutions.
const bool pin_pressure = true;
if (pin_pressure)
{
const unsigned int pressure_node = 0;
// const Real p_value = pref;
for (unsigned int c=0; c<elem->n_nodes(); c++)
if (elem->node(c) == pressure_node)
{
Kpp(c,c) += penalty;
//~ Fp(c) += penalty*p_value;
}
} //if(pin pressure)
} // end boundary condition section
#endif
dof_map.constrain_element_matrix (Ke, dof_indices);
// dof_map.constrain_element_matrix (Ke, dof_indices, dof_indices);
// Add the element matrix to the system Jacobian.
jacobian.add_matrix (Ke, dof_indices);
} // for ( ; el != end_el; ++el)
// That's it.
}
// Here we compute the residual R(x) = K(x)*x - f. The current solution
// x is passed in the soln vector
void compute_residual (const NumericVector<Number>& soln,
NumericVector<Number>& residual,
NonlinearImplicitSystem& /*sys*/)
{
EquationSystems &es = *_equation_system;
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
//~ libmesh_assert_equal_to (dim, 2);
// Get a reference to the NonlinearImplicitSystem we are solving
TransientNonlinearImplicitSystem& system =
es.get_system<TransientNonlinearImplicitSystem>("Navier-Stokes");
// Numeric ids corresponding to each variable in the system
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int p_var = system.variable_number ("p");
FEType fe_vel_type = system.variable_type(u_var);
FEType fe_pres_type = system.variable_type(p_var);
// Build a Finite Element object of the specified type for
// the velocity variables.
AutoPtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type));
AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
// A reference to the \p DofMap object for this system. The \p DofMap
// object handles the index translation from node and element numbers
// to degree of freedom numbers. We will talk more about the \p DofMap
// in future examples.
const DofMap& dof_map = system.get_dof_map();
// A Gauss quadrature rule for numerical integration.
// Let the \p FEType object decide what order rule is appropriate.
QGauss qrule (dim, fe_vel_type.default_quadrature_order());
// Tell the finite element objects to use our quadrature rule.
fe_vel->attach_quadrature_rule (&qrule);
fe_pres->attach_quadrature_rule (&qrule);
// Here we define some references to cell-specific data that
// will be used to assemble the linear system.
// The element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe_vel->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe_vel->get_phi();
// The element shape function gradients for the velocity
// variables evaluated at the quadrature points.
const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();
// The element shape functions for the pressure variable
// evaluated at the quadrature points.
const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();
// The value of the linear shape function gradients at the quadrature points
// const std::vector<std::vector<RealGradient> >& dpsi =
fe_pres->get_dphi();
#if 1
// Declare a special finite element object for
// boundary integration.
AutoPtr<FEBase> fe_vel_face (FEBase::build(dim, fe_vel_type));
AutoPtr<FEBase> fe_pres_face (FEBase::build(dim, fe_pres_type));
// Boundary integration requires one quadraure rule,
// with dimensionality one less than the dimensionality
// of the element.
QGauss qface(dim-1, fe_vel_type.default_quadrature_order());
QGauss qface_pres(dim-1, fe_pres_type.default_quadrature_order());
// Tell the finte element object to use our
// quadrature rule.
fe_vel_face->attach_quadrature_rule (&qface);
fe_pres_face->attach_quadrature_rule (&qface_pres);
//c.f. miscellaneous_ex5
const std::vector<std::vector<Real> >& phi_face = fe_vel_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face =
fe_vel_face->get_dphi();
const std::vector<Real>& JxW_face = fe_vel_face->get_JxW();
const std::vector<Point>& qface_normals = fe_vel_face->get_normals();
const std::vector<Point>& qface_points = fe_vel_face->get_xyz();
//pressure terms
const std::vector<std::vector<Real> >& psi_face = fe_pres_face->get_phi();
const std::vector<std::vector<RealGradient> >& dpsi_face =
fe_pres_face->get_dphi();
const std::vector<Real>& JxW_face_pres = fe_pres_face->get_JxW();
const std::vector<Point>& qface_normals_pres = fe_pres_face->get_normals();
const std::vector<Point>& qface_points_pres = fe_pres_face->get_xyz();
#endif
// Define data structures to contain the resdual contributions
DenseVector<Number> Fe;
DenseSubVector<Number>
Fu(Fe),
Fv(Fe),
Fp(Fe);
// This vector will hold the degree of freedom indices for
// the element. These define where in the global system
// the element degrees of freedom get mapped.
std::vector<dof_id_type> dof_indices;
std::vector<dof_id_type> dof_indices_u;
std::vector<dof_id_type> dof_indices_v;
std::vector<dof_id_type> dof_indices_p;
// ???
const Real dt = es.parameters.get<Real>("dt");
const Real visc = es.parameters.get<Real>("viscosity");
const Real pref = es.parameters.get<Real>("pref");
const Real penalty = es.parameters.get<Real>("penalty");
const Real penalty_pin = es.parameters.get<Real>("penalty_pin");
const bool pin_pressure=es.parameters.get<bool>("pin_pressure");
const unsigned int pin_pressure_node=es.parameters.get<unsigned
int>("pin_pressure_node");
residual.zero();
// Old solution
//~ const NumericVector<Number>& u_old = *old_local_solution;
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
// Store a pointer to the element we are currently
// working on. This allows for nicer syntax later.
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
dof_map.dof_indices (elem, dof_indices_u, u_var);
dof_map.dof_indices (elem, dof_indices_v, v_var);
dof_map.dof_indices (elem, dof_indices_p, p_var);
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_u_dofs = dof_indices_u.size();
const unsigned int n_v_dofs = dof_indices_v.size();
const unsigned int n_p_dofs = dof_indices_p.size();
// Compute the element-specific data for the current
// element. This involves computing the location of the
// quadrature points (q_point) and the shape functions
// (phi, dphi) for the current element.
fe_vel->reinit (elem);
fe_pres->reinit (elem);
// We use the resize member here because
// the number of degrees of freedom might have changed from
// the last element. Note that this will be the case if the
// element type is different (i.e. the last element was a
// triangle, now we are on a quadrilateral).
Fe.resize (n_dofs);
Fu.reposition (u_var*n_u_dofs, n_u_dofs);
Fv.reposition (v_var*n_u_dofs, n_v_dofs);
Fp.reposition (p_var*n_u_dofs, n_p_dofs);
// Now we will build the residual. This involves
// the construction of the matrix K and multiplication of it
// with the current solution x. We rearrange this into two loops:
// In the first, we calculate only the contribution of
// K_ij*x_j which is independent of the row i. In the second loops,
// we multiply with the row-dependent part and add it to the element
// residual.
Number u_old,v_old,p_old;
Number u,v,p,u_x,v_y,w_z;
Number jxw;
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
jxw=JxW[qp];
u_old=0.;
v_old=0.;
p_old=0.;
u=0.;
v=0.;
p=0.;
// Gradient grad_u_old,grad_v_old;
Gradient gradu,gradv;
for (unsigned int l=0; l<n_u_dofs; l++)
{
u_old+= phi[l][qp]*system.old_solution(dof_indices_u[l]);
v_old+= phi[l][qp]*system.old_solution(dof_indices_v[l]);
//~ grad_u_old +=
dphi[l][qp]*system.old_solution(dof_indices_u[l]);
//~ grad_v_old +=
dphi[l][qp]*system.old_solution(dof_indices_v[l]);
u += phi[l][qp]*soln(dof_indices_u[l]);
v += phi[l][qp]*soln(dof_indices_v[l]);
gradu.add_scaled (dphi[l][qp],soln (dof_indices_u[l]));
gradv.add_scaled (dphi[l][qp],soln (dof_indices_v[l]));
}
for (unsigned int l=0; l<n_p_dofs; l++)
{
p_old += psi[l][qp]*system.old_solution (dof_indices_p[l]);
p += psi[l][qp]*soln (dof_indices_p[l]);
}
// const NumberVectorValue U_old (u_old, v_old);
const NumberVectorValue U (u, v);
for (unsigned int i=0; i<n_u_dofs; i++)
{
Fu(i)+=jxw*((u-u_old)*phi[i][qp] //mass matrix
+dt*(visc*dphi[i][qp]*gradu // diffusion
+U*gradu*phi[i][qp]
-p*dphi[i][qp](0)
)); //convection
Fv(i)+=jxw*((v-v_old)*phi[i][qp] //mass matrix
+dt*(visc*dphi[i][qp]*gradv //diffusion
+U*gradv*phi[i][qp]
-p*dphi[i][qp](1)
)); //convection
} // i loop n_u_dofs
for (unsigned int i=0; i<n_p_dofs; i++)
{
Fp(i) -=jxw*dt*psi[i][qp]*(gradu(0)+gradv(1));
}
} //qprule loop in element
// boundary integrals
// The following loops over the sides of the element.
// If the element has no neighbor on a side then that
// side MUST live on a boundary of the domain.
for (unsigned int side=0; side<elem->n_sides(); side++)
{
if (elem->neighbor(side) == NULL)
{
boundary_id_type
bid=mesh.get_boundary_info().boundary_id(elem,side);
AutoPtr<Elem> _side (elem->build_side(side));
// Compute the shape function values on the element face.
fe_vel_face->reinit(elem, side);
fe_pres_face->reinit(elem, side);
// const unsigned int elem_b_order = static_cast<unsigned int>
(fe_vel_face->get_order());
// const double h_elem = elem->volume()/_side->volume() *
1./pow(elem_b_order, 2.);
// ******** Dirichlet BC by penalty ***************************
// only for lid-driven cavity , apply Diri BC
if (bid<4 && bid>=0)
{
// std::cout<<"phi_face (nqp x
nphi_face="<<qface.n_points()<<"x"<<phi_face.size()<<"):"<<std::endl;
// std::cout<<"psi_face_pres (nqp x
nphi_face="<<qface_pres.n_points()<<"x"<<psi_face.size()<<"):"<<std::endl;
// Boundary ids are set internally by
// build_square().
// 0=bottom // 1=right // 2=top // 3=left
// Set u = 1 on the top boundary, 0 everywhere else
const Real u_value =
(mesh.get_boundary_info().has_boundary_id(elem,side,2))
? 1. : 0.;
// Set v = 0 everywhere
const Real v_value = 0.;
const unsigned int elem_b_order = static_cast<unsigned int>
(fe_vel_face->get_order());
const double h_elem = elem->volume()/_side->volume() *
1./pow(elem_b_order, 2.);
// Loop over the nodes on the side.
Number jxw;
Number temp;
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Gradient gradu,gradv;
jxw=JxW_face[qp];
u=0.;
v=0.;
p=0.;
for (unsigned int l=0; l<phi_face.size(); l++)//???
{
//~ u_old+=
phi[l][qp]*system.old_solution(dof_indices_u[l]);
//~ v_old+=
phi[l][qp]*system.old_solution(dof_indices_v[l]);
//~ grad_u_old +=
dphi[l][qp]*system.old_solution(dof_indices_u[l]);
//~ grad_v_old +=
dphi[l][qp]*system.old_solution(dof_indices_v[l]);
u += phi_face[l][qp]*soln(dof_indices_u[l]);
v += phi_face[l][qp]*soln(dof_indices_v[l]);
gradu.add_scaled (dphi_face[l][qp],soln
(dof_indices_u[l]));
gradv.add_scaled (dphi_face[l][qp],soln
(dof_indices_v[l]));
}
// for (unsigned int l=0; l<n_p_dofs; l++)
for (unsigned int l=0; l<psi_face.size(); l++)//???
{
//~ p_old += psi[l][qp]*system.old_solution
(dof_indices_p[l]);
p += psi_face[l][qp]*soln(dof_indices_p[l]);
}
// const NumberVectorValue U_old (u_old, v_old);
const NumberVectorValue U (u, v);
const NumberVectorValue U_bc (u_value, v_value);
for (unsigned int i=0; i<phi_face.size(); i++)
{
Fu(i)+=jxw*dt* penalty/h_elem * phi_face[i][qp] *
(u-u_value);
//Fu(i)-=jxw*dt* visc*(//phi_face[i][qp] *
(gradu*qface_normals[qp]) // term A_u
//+
// (u-u_value) *
(dphi_face[i][qp]*qface_normals[qp])); // C_u-F_u
Fv(i)+=jxw*dt* penalty/h_elem * phi_face[i][qp] *
(v-v_value);
//Fv(i)-=jxw*dt* visc*(//phi_face[i][qp] *
(gradv*qface_normals[qp]) // term A_u
//+
// (v-v_value) *
(dphi_face[i][qp]*qface_normals[qp])); // term C_v-F_v, both for consistency
}
#if 0 //B Block
for (unsigned int i=0; i<phi_face.size(); i++)
{
temp=jxw * dt* p*phi_face[i][qp];
Fu(i)+= temp*qface_normals[qp](0); // B_u
Fv(i)+= temp*qface_normals[qp](1); // B_v
}
#endif
#if 0 //D-G block
for (unsigned int i=0; i<psi_face.size(); i++)
{
temp=jxw * dt* psi_face[i][qp];
Fp(i)+=temp*(U-U_bc)*qface_normals[qp]; //D-G
}
#endif
}// qp loop in qface
#if 0
for (unsigned int ns=0; ns<_side->n_nodes(); ns++)
{
u=0.;
v=0.;
// Find the node on the element matching this node on
// the side. That defined where in the element matrix
// the boundary condition will be applied.
for (unsigned int n=0; n<elem->n_nodes(); n++)
if (elem->node(n) == _side->node(ns))
{
//~ // Right-hand-side contribution.
Fu(n) += penalty*(u-u_value);
Fv(n) += penalty*(v-v_value);
}
} // end face node loop
// ******** Neumann BC *******************************
// Loop over the face quadrature points for integration.
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
// This is the right-hand-side contribution (f),
// which has to be subtracted from the current residual
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
} // qp loop in qface
#endif
} // if (bid<4 && bid>=0)
} //if (elem->neighbor(side) == NULL)
} // side loop , if (elem->neighbor(side) == NULL)
// Pin the pressure to zero at global node number "pressure_node".
// This effectively removes the non-trivial null space of constant
// pressure solutions.
if (pin_pressure)
{
for (unsigned int c=0; c<elem->n_nodes(); c++)
{
if (elem->node(c) == pin_pressure_node)
{
Fp(c) += penalty_pin*(soln(dof_indices_p[c])-pref);
}
} // c loop
} // pin pressure
dof_map.constrain_element_vector (Fe, dof_indices);
residual.add_vector (Fe, dof_indices);
} // end of element loop
// That's it.
}
------------------------------------------------------------------------------
Dive into the World of Parallel Programming The Go Parallel Website, sponsored
by Intel and developed in partnership with Slashdot Media, is your hub for all
things parallel software development, from weekly thought leadership blogs to
news, videos, case studies, tutorials and more. Take a look and join the
conversation now. http://goparallel.sourceforge.net/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
