From: <tim@ce...>  20070115 13:43:19
Attachments:
test.cpp

Dear libMesh developer team, Currently, I'm working on some NavierStokes application with solid obstacles in the computational domain. My implementation of the obstacles consists of fixing the velocity to zero on all dofs in the respective cells. Please find attached a test program, which is in fact just a slight modification of ex13; the changes are marked with "qqq". (My real application is more complicated, but this suffices to demsonstrate my problem.) The problem: When using a high resolution grid, no solution is obtained. Either PETSc crashes with a "Detected zero pivot in LU factorization" message, or the solver just does not converge. Everything works well on a lowresolution grid. Also, everything works well on a highresolution grid if the obstacle is disabled. Of course, another way of implementing the obstacle would be to remove the obstacle cells from the computational domain, but unfortunately, in my real application, I have to solve an additional PDE on the whole domain, including the obstacle. Also, the shape and size of the obstacle vary in time. Do you have any ideas what is going wrong and how I can do it better? Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen Universitaetsallee 29 (Office 3.13) tim@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
From: Roy Stogner <roystgnr@ic...>  20070115 15:42:35

On Mon, 15 Jan 2007, Tim Kr=F6ger wrote: > Do you have any ideas what is going wrong and how I can do it better? This might just be the second question in a week whose answer is "we don't constrain away the constant pressure solution in ex13". Unless you've pinned a free pressure node along with all your obstacle velocity nodes, your resulting matrix will be one rank deficient, and it's certainly possible that other changes you make to the problem will render the solver (or much more likely, preconditioner) unable to handle that. You might try playing around with different preconditioners before modifying your system, though.  Roy Stogner 
From: John Peterson <peterson@cf...>  20070115 15:43:04

One thing that's missing from ex13 is pinning a single value of the pressure. This is a common (at least amongst LibMesh developers) hack for rectifying the nontrivial null space of constant pressures which is associated with the NS equations. Here's a code snippet that can go at the end of your bc routine to fix a single value of the pressure at node "0" (assuming that node 0 is on the boundary.) I will try to get this added to ex13 as well...= John // Pin the pressure to zero const unsigned int pressure=5Fnode=3D0; const Real p=5Fvalue =3D 0.0; for (unsigned int c=3D0; c<elem>n=5Fnodes(); c++) if (elem>node(c) =3D=3D pressure=5Fnode) { =09// here(); =09// std::cout << elem>point(c) << std::endl; =09Kpp(c,c) +=3D penalty; =09Fp(c) +=3D penalty*p=5Fvalue; } Tim Kr=F6ger writes: > Dear libMesh developer team, >=20 > Currently, I'm working on some NavierStokes application with solid=20= > obstacles in the computational domain. My implementation of the=20 > obstacles consists of fixing the velocity to zero on all dofs in the= =20 > respective cells. Please find attached a test program, which is in=20= > fact just a slight modification of ex13; the changes are marked with= =20 > "qqq". (My real application is more complicated, but this suffices = to=20 > demsonstrate my problem.) >=20 > The problem: When using a high resolution grid, no solution is=20 > obtained. Either PETSc crashes with a "Detected zero pivot in LU=20= > factorization" message, or the solver just does not converge. >=20 > Everything works well on a lowresolution grid. Also, everything=20= > works well on a highresolution grid if the obstacle is disabled. >=20 > Of course, another way of implementing the obstacle would be to remo= ve=20 > the obstacle cells from the computational domain, but unfortunately,= =20 > in my real application, I have to solve an additional PDE on the who= le=20 > domain, including the obstacle. Also, the shape and size of the=20 > obstacle vary in time. >=20 > Do you have any ideas what is going wrong and how I can do it better= =3F >=20 > Best Regards, >=20 > Tim >=20 > =20 > Dr. Tim Kroeger > CeVis  Center of Complex Systems and Visualization > University of Bremen > Universitaetsallee 29 (Office 3.13) tim@... > D28359 Bremen Phone +494212187710 > Germany Fax +494212184236/* = $Id: ex13.C,v 1.16 2006/10/23 18:46:06 jwpeterson Exp $ */ >=20 > /* The Next Great Finite Element Library. */ > /* Copyright (C) 2003 Benjamin S. Kirk */ >=20 > /* 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= . */ >=20 > /* 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. */ >=20 > /* You should have received a copy of the GNU Lesser General Public = */ > /* License along with this library; if not, write to the Free Softwa= re */ > /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 021111= 307 USA */ >=20 > // <h1>Example 13  Unsteady NavierStokes Equations  Unsteady Non= linear Systems of Equations</h1> > // > // This example shows how a simple, unsteady, nonlinear system of e= quations > // can be solved in parallel. The system of equations are the fami= liar > // NavierStokes equations for lowspeed incompressible fluid flow.= This > // example introduces the concept of the inner nonlinear loop for e= ach > // timestep, and requires a good deal of linear algebra numbercrun= ching > // at each step. If you have the General Mesh Viewer (GMV) install= ed, > // 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 i= f 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=5Ftimesteps and dt variables below. >=20 > // C++ include files that we need > #include <iostream> > #include <algorithm> > #include <math.h> >=20 > // Basic include file needed for the mesh functionality. > #include "libmesh.h" > #include "mesh.h" > #include "mesh=5Fgeneration.h" > #include "gmv=5Fio.h" > #include "equation=5Fsystems.h" > #include "fe.h" > #include "quadrature=5Fgauss.h" > #include "dof=5Fmap.h" > #include "sparse=5Fmatrix.h" > #include "numeric=5Fvector.h" > #include "dense=5Fmatrix.h" > #include "dense=5Fvector.h" > #include "linear=5Fimplicit=5Fsystem.h" > #include "transient=5Fsystem.h" > #include "perf=5Flog.h" > #include "boundary=5Finfo.h" >=20 > // Some (older) compilers do not offer full stream=20 > // functionality, OStringStream works around this. > #include "o=5Fstring=5Fstream.h" >=20 > // For systems of equations the \p DenseSubMatrix > // and \p DenseSubVector provide convenient ways for > // assembling the element matrix and vector on a > // componentbycomponent basis. > #include "dense=5Fsubmatrix.h" > #include "dense=5Fsubvector.h" >=20 > // The definition of a geometric element > #include "elem.h" >=20 > // Function prototype. This function will assemble the system > // matrix and righthandside. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname); >=20 > // The main program. > int main (int argc, char** argv) > { > // Initialize libMesh. > libMesh::init (argc, argv); > { =20 > // Set the dimensionality of the mesh =3D 2 > const unsigned int dim =3D 2; =20 > =20 > // Create a twodimensional mesh. > Mesh mesh (dim); > =20 > // Use the MeshTools::Generation mesh generator to create a unif= orm > // grid on the square [1,1]^D. We instruct the mesh generator > // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 > // elements in 3D. Building these higherorder elements allows > // us to use higherorder approximation, as in example 3. > MeshTools::Generation::build=5Fsquare (mesh, > =09=09=09=09=09 //qqq=09=09=09=09=09 20, 20, > =09=09=09=09=09 100, 100,//qqq > =09=09=09=09=09 0., 1., > =09=09=09=09=09 0., 1., > =09=09=09=09=09 QUAD9); > =20 > // Print information about the mesh to the screen. > mesh.print=5Finfo(); > =20 > // Create an equation systems object. > EquationSystems equation=5Fsystems (mesh); > =20 > // Declare the system and its variables. > { > // Creates a transient system named "NavierStokes" > TransientLinearImplicitSystem & system =3D=20 > =09equation=5Fsystems.add=5Fsystem<TransientLinearImplicitSystem> ("= NavierStokes"); > =20 > // Add the variables "u" & "v" to "NavierStokes". They > // will be approximated using secondorder approximation. > system.add=5Fvariable ("u", SECOND); > system.add=5Fvariable ("v", SECOND); >=20 > // Add the variable "p" to "NavierStokes". This will > // be approximated with a firstorder basis, > // providing an LBBstable pressurevelocity pair. > system.add=5Fvariable ("p", FIRST); >=20 > // Give the system a pointer to the matrix assembly > // function. > system.attach=5Fassemble=5Ffunction (assemble=5Fstokes); > =20 > // Initialize the data structures for the equation system. > equation=5Fsystems.init (); >=20 > equation=5Fsystems.parameters.set<unsigned int>("linear solver= maximum iterations") =3D 250; > equation=5Fsystems.parameters.set<Real> ("linear solver= tolerance") =3D std::sqrt(TOLERANCE); > =20 > // Prints information about the system to the screen. > equation=5Fsystems.print=5Finfo(); > } >=20 > // Create a performancelogging object for this example > PerfLog perf=5Flog("Example 13"); > =20 > // Now we begin the timestep loop to compute the timeaccurate > // solution of the equations. > const Real dt =3D 0.005; > Real time =3D 0.0; > const unsigned int n=5Ftimesteps =3D 15; >=20 > // The number of steps and the stopping criterion are also requi= red > // for the nonlinear iterations. > const unsigned int n=5Fnonlinear=5Fsteps =3D 15; > const Real nonlinear=5Ftolerance =3D 1.e2; > =20 > // Tell the system of equations what the timestep is by using > // the set=5Fparameter function. The matrix assembly routine ca= n > // then reference this parameter. > equation=5Fsystems.parameters.set<Real> ("dt") =3D dt; >=20 > // Get a reference to the Stokes system to use later. > TransientLinearImplicitSystem& navier=5Fstokes=5Fsystem =3D > =09 equation=5Fsystems.get=5Fsystem<TransientLinearImplicitSystem>(= "NavierStokes"); >=20 > // 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=5Fnonlinear=5Fsoln (navier=5Fstokes=5Fsystem.solution>cl= one()); >=20 > for (unsigned int t=5Fstep=3D0; t=5Fstep<n=5Ftimesteps; ++t=5Fst= ep) > { > =09// Incremenet the time counter, set the time and the > =09// time step size as parameters in the EquationSystem. > =09time +=3D dt; >=20 > =09// Let the system of equations know the current time. > =09// This might be necessary for a timedependent forcing > =09// function for example. > =09equation=5Fsystems.parameters.set<Real> ("time") =3D time; >=20 > =09// A pretty update message > =09std::cout << " Solving time step " << t=5Fstep << ", time =3D " <= < time << std::endl; >=20 > =09// Now we need to update the solution vector from the > =09// previous time step. This is done directly through > =09// the reference to the Stokes system. > =09*navier=5Fstokes=5Fsystem.old=5Flocal=5Fsolution =3D *navier=5Fst= okes=5Fsystem.current=5Flocal=5Fsolution; >=20 > =09// Now we begin the nonlinear loop > =09for (unsigned int l=3D0; l<n=5Fnonlinear=5Fsteps; ++l) > =09 { > =09 // Update the nonlinear solution. > =09 last=5Fnonlinear=5Fsoln>zero(); > =09 last=5Fnonlinear=5Fsoln>add(*navier=5Fstokes=5Fsystem.soluti= on); > =09 =20 > =09 // Assemble & solve the linear system. > =09 perf=5Flog.start=5Fevent("linear solve"); > =09 equation=5Fsystems.get=5Fsystem("NavierStokes").solve(); > =09 perf=5Flog.stop=5Fevent("linear solve"); >=20 > =09 // Compute the difference between this solution and the last > =09 // nonlinear iterate. > =09 last=5Fnonlinear=5Fsoln>add (1., *navier=5Fstokes=5Fsystem.= solution); >=20 > =09 // Close the vector before computing its norm > =09 last=5Fnonlinear=5Fsoln>close(); >=20 > =09 // Compute the l2 norm of the difference > =09 const Real norm=5Fdelta =3D last=5Fnonlinear=5Fsoln>l2=5Fnor= m(); >=20 > =09 // Print out convergence information for the linear and > =09 // nonlinear iterations. > =09 std::cout << "Linear solver converged at step: " > =09=09 << navier=5Fstokes=5Fsystem.n=5Flinear=5Fiterations() > =09=09 << ", final residual: " > =09=09 << navier=5Fstokes=5Fsystem.final=5Flinear=5Fresidual() > =09=09 << " Nonlinear convergence: u  u=5Fold =3D " > =09=09 << norm=5Fdelta > =09=09 << std::endl; >=20 > =09 // Terminate the solution iteration if the difference between= > =09 // this iteration and the last is sufficiently small. > =09 if (norm=5Fdelta < nonlinear=5Ftolerance) > =09 { > =09=09std::cout << " Nonlinear solver converged at step " > =09=09=09 << l > =09=09=09 << std::endl; > =09=09break; > =09 } > =09 =20 > =09 } // end nonlinear loop > =09 > =09// Write out every nth timestep to file. > =09const unsigned int write=5Finterval =3D 1; > =09 > =09if ((t=5Fstep+1)%write=5Finterval =3D=3D 0) > =09 { > =09 OStringStream file=5Fname; >=20 > =09 // We write the file name in the gmv autoread format. > =09 file=5Fname << "out.gmv."; > =09 OSSRealzeroright(file=5Fname,3,0, t=5Fstep + 1); > =09 =20 > =09 GMVIO(mesh).write=5Fequation=5Fsystems (file=5Fname.str(), > =09=09=09=09=09=09equation=5Fsystems); > =09 } > } // end timestep loop. > } > =20 > // All done. =20 > return libMesh::close (); > } >=20 >=20 >=20 >=20 >=20 >=20 > // The matrix assembly function to be called at each time step to > // prepare for the linear solve. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname) > { > // It is a good idea to make sure we are assembling > // the proper system. > assert (system=5Fname =3D=3D "NavierStokes"); > =20 > // Get a constant reference to the mesh object. > const Mesh& mesh =3D es.get=5Fmesh(); > =20 > // The dimension that we are running > const unsigned int dim =3D mesh.mesh=5Fdimension(); > =20 > // Get a reference to the Stokes system object. > TransientLinearImplicitSystem & navier=5Fstokes=5Fsystem =3D > es.get=5Fsystem<TransientLinearImplicitSystem> ("NavierStokes")= ; >=20 > // Numeric ids corresponding to each variable in the system > const unsigned int u=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("u"); > const unsigned int v=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("v"); > const unsigned int p=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("p"); > =20 > // Get the Finite Element type for "u". Note this will be > // the same as the type for "v". > FEType fe=5Fvel=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Ftyp= e(u=5Fvar); > =20 > // Get the Finite Element type for "p". > FEType fe=5Fpres=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Fty= pe(p=5Fvar); >=20 > // Build a Finite Element object of the specified type for > // the velocity variables. > AutoPtr<FEBase> fe=5Fvel (FEBase::build(dim, fe=5Fvel=5Ftype)); > =20 > // Build a Finite Element object of the specified type for > // the pressure variables. > AutoPtr<FEBase> fe=5Fpres (FEBase::build(dim, fe=5Fpres=5Ftype)); > =20 > // A Gauss quadrature rule for numerical integration. > // Let the \p FEType object decide what order rule is appropriate.= > QGauss qrule (dim, fe=5Fvel=5Ftype.default=5Fquadrature=5Forder())= ; >=20 > // Tell the finite element objects to use our quadrature rule. > fe=5Fvel>attach=5Fquadrature=5Frule (&qrule); > fe=5Fpres>attach=5Fquadrature=5Frule (&qrule); > =20 > // Here we define some references to cellspecific data that > // will be used to assemble the linear system. > // > // The element Jacobian * quadrature weight at each integration po= int. =20 > const std::vector<Real>& JxW =3D fe=5Fvel>get=5FJxW(); >=20 > // The element shape functions evaluated at the quadrature points.= > const std::vector<std::vector<Real> >& phi =3D fe=5Fvel>get=5Fphi= (); >=20 > // The element shape function gradients for the velocity > // variables evaluated at the quadrature points. > const std::vector<std::vector<RealGradient> >& dphi =3D fe=5Fvel>= get=5Fdphi(); >=20 > // The element shape functions for the pressure variable > // evaluated at the quadrature points. > const std::vector<std::vector<Real> >& psi =3D fe=5Fpres>get=5Fph= i(); >=20 > // The value of the linear shape function gradients at the quadrat= ure points > // const std::vector<std::vector<RealGradient> >& dpsi =3D fe=5Fpr= es>get=5Fdphi(); > =20 > // A reference to the \p DofMap object for this system. The \p Do= fMap > // object handles the index translation from node and element numb= ers > // to degree of freedom numbers. We will talk more about the \p D= ofMap > // in future examples. > const DofMap & dof=5Fmap =3D navier=5Fstokes=5Fsystem.get=5Fdof=5F= map(); >=20 > // Define data structures to contain the element matrix > // and righthandside vector contribution. Following > // basic finite element terminology we will denote these > // "Ke" and "Fe". > DenseMatrix<Number> Ke; > DenseVector<Number> Fe; >=20 > DenseSubMatrix<Number> > Kuu(Ke), Kuv(Ke), Kup(Ke), > Kvu(Ke), Kvv(Ke), Kvp(Ke), > Kpu(Ke), Kpv(Ke), Kpp(Ke); >=20 > DenseSubVector<Number> > Fu(Fe), > Fv(Fe), > Fp(Fe); >=20 > // 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<unsigned int> dof=5Findices; > std::vector<unsigned int> dof=5Findices=5Fu; > std::vector<unsigned int> dof=5Findices=5Fv; > std::vector<unsigned int> dof=5Findices=5Fp; >=20 > // Find out what the timestep size parameter is from the system, a= nd > // the value of theta for the theta method. We use implicit Euler= (theta=3D1) > // for this simulation even though it is only firstorder accurate= in time. > // The reason for this decision is that the secondorder CrankNic= olson > // method is notoriously oscillatory for problems with discontinuo= us > // initial data such as the liddriven cavity. Therefore, > // we sacrifice accuracy in time for stability, but since the solu= tion > // reaches steady state relatively quickly we can afford to take s= mall > // timesteps. If you monitor the initial nonlinear residual for t= his > // simulation, you should see that it is monotonically decreasing = in time. > const Real dt =3D es.parameters.get<Real>("dt"); > // const Real time =3D es.parameters.get<Real>("time"); > const Real theta =3D 1.; > =20 > // Now we will loop over all the elements in the mesh that > // live on the local processor. We will compute the element > // matrix and righthandside contribution. Since the mesh > // will be refined we want to only consider the ACTIVE elements, > // hence we use a variant of the \p active=5Felem=5Fiterator. > // const=5Factive=5Flocal=5Felem=5Fiterator el (mesh.ele= ments=5Fbegin()); > // const const=5Factive=5Flocal=5Felem=5Fiterator end=5Fel (mesh.e= lements=5Fend()); >=20 > MeshBase::const=5Felement=5Fiterator el =3D mesh.active=5F= local=5Felements=5Fbegin(); > const MeshBase::const=5Felement=5Fiterator end=5Fel =3D mesh.activ= e=5Flocal=5Felements=5Fend();=20 > =20 > for ( ; el !=3D end=5Fel; ++el) > { =20 > // Store a pointer to the element we are currently > // working on. This allows for nicer syntax later. > const Elem* elem =3D *el; > =20 > // Get the degree of freedom indices for the > // current element. These define where in the global > // matrix and righthandside this element will > // contribute to. > dof=5Fmap.dof=5Findices (elem, dof=5Findices); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fu, u=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fv, v=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fp, p=5Fvar); >=20 > const unsigned int n=5Fdofs =3D dof=5Findices.size(); > const unsigned int n=5Fu=5Fdofs =3D dof=5Findices=5Fu.size();=20= > const unsigned int n=5Fv=5Fdofs =3D dof=5Findices=5Fv.size();=20= > const unsigned int n=5Fp=5Fdofs =3D dof=5Findices=5Fp.size(); > =20 > // Compute the elementspecific data for the current > // element. This involves computing the location of the > // quadrature points (q=5Fpoint) and the shape functions > // (phi, dphi) for the current element. > fe=5Fvel>reinit (elem); > fe=5Fpres>reinit (elem); >=20 > // Zero the element matrix and righthand side 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=5Fdofs, n=5Fdofs); > Fe.resize (n=5Fdofs); >=20 > // Reposition the submatrices... The idea is this: > // > //     > //  Kuu Kuv Kup   Fu  > // Ke =3D  Kvu Kvv Kvp ; Fe =3D  Fv  > //  Kpu Kpv Kpp   Fp  > //     > // > // The \p DenseSubMatrix.repostition () member takes the > // (row=5Foffset, column=5Foffset, row=5Fsize, column=5Fsize).= > //=20 > // Similarly, the \p DenseSubVector.reposition () member > // takes the (row=5Foffset, row=5Fsize) > Kuu.reposition (u=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fu=5Fdofs); > Kuv.reposition (u=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fv=5Fdofs); > Kup.reposition (u=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fp=5Fdofs); > =20 > Kvu.reposition (v=5Fvar*n=5Fv=5Fdofs, u=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fu=5Fdofs); > Kvv.reposition (v=5Fvar*n=5Fv=5Fdofs, v=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fv=5Fdofs); > Kvp.reposition (v=5Fvar*n=5Fv=5Fdofs, p=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fp=5Fdofs); > =20 > Kpu.reposition (p=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fu=5Fdofs); > Kpv.reposition (p=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fv=5Fdofs); > Kpp.reposition (p=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fp=5Fdofs); >=20 > Fu.reposition (u=5Fvar*n=5Fu=5Fdofs, n=5Fu=5Fdofs); > Fv.reposition (v=5Fvar*n=5Fu=5Fdofs, n=5Fv=5Fdofs); > Fp.reposition (p=5Fvar*n=5Fu=5Fdofs, n=5Fp=5Fdofs); >=20 > // Now we will build the element matrix and righthandside. > // Constructing the RHS requires the solution and its > // gradient from the previous timestep. This must be > // calculated at each quadrature point by summing the > // solution degreeoffreedom values by the appropriate > // weight functions. > for (unsigned int qp=3D0; qp<qrule.n=5Fpoints(); qp++) > =09{ > =09 // Values to hold the solution & its gradient at the previous t= imestep. > =09 Number u =3D 0., u=5Fold =3D 0.; > =09 Number v =3D 0., v=5Fold =3D 0.; > =09 Number p=5Fold =3D 0.; > =09 Gradient grad=5Fu, grad=5Fu=5Fold; > =09 Gradient grad=5Fv, grad=5Fv=5Fold; > =09 =20 > =09 // Compute the velocity & its gradient from the previous timest= ep > =09 // and the old Newton iterate. > =09 for (unsigned int l=3D0; l<n=5Fu=5Fdofs; l++) > =09 { > =09 // From the old timestep: > =09 u=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fu[l]); > =09 v=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fv[l]); > =09 grad=5Fu=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fv[l])); >=20 > =09 // From the previous Newton iterate: > =09 u +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fu[l]);=20 > =09 v +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fv[l]); > =09 grad=5Fu.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fv[l])); > =09 } >=20 > =09 // Compute the old pressure value at this quadrature point. > =09 for (unsigned int l=3D0; l<n=5Fp=5Fdofs; l++) > =09 { > =09 p=5Fold +=3D psi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fp[l]); > =09 } >=20 > =09 // Definitions for convenience. It is sometimes simpler to do = a > =09 // dot product if you have the full vector at your disposal. > =09 const NumberVectorValue U=5Fold (u=5Fold, v=5Fold); > =09 const NumberVectorValue U (u, v); > =09 const Number u=5Fx =3D grad=5Fu(0); > =09 const Number u=5Fy =3D grad=5Fu(1); > =09 const Number v=5Fx =3D grad=5Fv(0); > =09 const Number v=5Fy =3D grad=5Fv(1); > =09 =20 > =09 // First, an iloop over the velocity degrees of freedom. > =09 // We know that n=5Fu=5Fdofs =3D=3D n=5Fv=5Fdofs so we can comp= ute contributions > =09 // for both at the same time. > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Fu(i) +=3D JxW[qp]*(u=5Fold*phi[i][qp]  = // massmatrix term=20 > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fu=5Fold)*phi[i][qp] + // c= onvection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](0)  // press= ure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fu=5Fold*dphi[i][qp]) + // dif= fusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fu)*phi[i][qp]); // Newto= n term >=20 > =09=09 > =09 Fv(i) +=3D JxW[qp]*(v=5Fold*phi[i][qp]  = // massmatrix term > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fv=5Fold)*phi[i][qp] + // = convection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](1)  // pres= sure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fv=5Fold*dphi[i][qp]) + // di= ffusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fv)*phi[i][qp]); // Newt= on term > =09 =20 >=20 > =09 // Note that the Fp block is identically zero unless we are= using > =09 // some kind of artificial compressibility scheme... >=20 > =09 // Matrix contributions for the uu and vv couplings. > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kuu(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*u=5Fx*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kuv(i,j) +=3D JxW[qp]*theta*dt*u=5Fy*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09 =20 > =09=09 Kvv(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*v=5Fy*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kvu(i,j) +=3D JxW[qp]*theta*dt*v=5Fx*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09} >=20 > =09 // Matrix contributions for the up and vp couplings. > =09 for (unsigned int j=3D0; j<n=5Fp=5Fdofs; j++) > =09=09{ > =09=09 Kup(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](0));= > =09=09 Kvp(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](1));= > =09=09} > =09 } >=20 > =09 // Now an iloop over the pressure degrees of freedom. This co= de computes > =09 // the matrix entries due to the continuity equation. Note: To= maintain a > =09 // symmetric matrix, we may (or may not) multiply the continuit= y equation by > =09 // negative one. Here we do not. > =09 for (unsigned int i=3D0; i<n=5Fp=5Fdofs; i++) > =09 { > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kpu(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](0); > =09=09 Kpv(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](1); > =09=09} > =09 } > =09} // end of the quadrature point qploop >=20 > =20 > // At this point the interior element integration has > // been completed. However, we have not yet addressed > // boundary conditions. For this example we will only > // consider simple Dirichlet boundary conditions imposed > // via the penalty method. The penalty method used here > // is equivalent (for Lagrange basis functions) to lumping > // the matrix resulting from the L2 projection penalty > // approach introduced in example 3. > { > =09// The following loops over the sides of the element. > =09// If the element has no neighbor on a side then that > =09// side MUST live on a boundary of the domain. > =09for (unsigned int s=3D0; s<elem>n=5Fsides(); s++) > =09 if (elem>neighbor(s) =3D=3D NULL) > =09 { > =09 // Get the boundary ID for side 's'. > =09 // These are set internally by build=5Fsquare(). > =09 // 0=3Dbottom > =09 // 1=3Dright > =09 // 2=3Dtop > =09 // 3=3Dleft > =09 short int bc=5Fid =3D mesh.boundary=5Finfo>boundary=5Fid (= elem,s); > =09 if (bc=5Fid=3D=3DBoundaryInfo::invalid=5Fid) > =09=09 error(); >=20 > =09 =20 > =09 AutoPtr<Elem> side (elem>build=5Fside(s)); > =09 =09 =20 > =09 // Loop over the nodes on the side. > =09 for (unsigned int ns=3D0; ns<side>n=5Fnodes(); ns++) > =09=09{ > =09=09 // The penalty value. \f$ \frac{1}{\epsilon \f$ > =09=09 const Real penalty =3D 1.e10; > =09=09 =20 > =09=09 // Get the boundary values. > =09=09 =20 > =09=09 // Set u =3D 1 on the top boundary, 0 everywhere else > =09=09 const Real u=5Fvalue =3D (bc=5Fid=3D=3D2) =3F 1. : 0.; > =09=09 =20 > =09=09 // Set v =3D 0 everywhere > =09=09 const Real v=5Fvalue =3D 0.; > =09=09 =20 > =09=09 // Find the node on the element matching this node on > =09=09 // the side. That defined where in the element matrix > =09=09 // the boundary condition will be applied. > =09=09 for (unsigned int n=3D0; n<elem>n=5Fnodes(); n++) > =09=09 if (elem>node(n) =3D=3D side>node(ns)) > =09=09 { > =09=09=09// Matrix contribution. > =09=09=09Kuu(n,n) +=3D penalty; > =09=09=09Kvv(n,n) +=3D penalty; > =09=09 =09=09 =20 > =09=09=09// Righthandside contribution. > =09=09=09Fu(n) +=3D penalty*u=5Fvalue; > =09=09=09Fv(n) +=3D penalty*v=5Fvalue; > =09=09 } > =09=09} // end face node loop=09 =20 > =09 } // end if (elem>neighbor(side) =3D=3D NULL) > } // end boundary condition section=09 =20 > =20 > /* qqq begin */ > if(fabs(elem>centroid()(0)0.5)<=3D0.4 && fabs(elem>centroid= ()(1)0.5)<=3D0.4) > =09{ > =09 const Real penalty =3D 1.e7; > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Kuu(i,i) +=3D penalty; > =09 Kvv(i,i) +=3D penalty; > =09 } > =09} > /* qqq end */ >=20 > // The element matrix and righthandside are now built > // for this element. Add them to the global matrix and > // righthandside vector. The \p PetscMatrix::add=5Fmatrix()= > // and \p PetscVector::add=5Fvector() members do this for us. > navier=5Fstokes=5Fsystem.matrix>add=5Fmatrix (Ke, dof=5Findic= es); > navier=5Fstokes=5Fsystem.rhs>add=5Fvector (Fe, dof=5Findic= es); > } // end of element loop > =20 > // That's it. > return; > } > =  > Take Surveys. Earn Cash. Influence the Future of IT > Join SourceForge.net's Techsay panel and you'll get the chance to sh= are your > opinions on IT & business topics through brief surveys  and earn ca= sh > http://www.techsay.com/default.php=3Fpage=3Djoin.php&p=3Dsourceforge= &CID=3DDEVDEV=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: <tim@ce...>  20070116 08:14:30
Attachments:
test.cpp

Dear Roy and John, thank you for your nearly simultaneous answer. (Did you use a stopwatch to get that right?) (: Actually, constraining away the onedimensional null space did not cure the problem in my case. Please find attached a modified version of my test program that includes John's code snippet. Again, all the (very little) places where the program differs from ex13 are marked with "qqq". The behaviour is the same as before, i.e. the program works well without obstacle (rough and fine grid), and it works well with obstacle on the rough grid (20x20), but it crashes with obstacle on the fine grid (100x100). I pinned the pressure to 100.0 at node 0, but pinning it to 0.0 instead does not change the problem. It would be nice if you at least tested whether this phenomenon also happens with your installation, in order to exclude version conflicts. Thank you very much in advance, Tim 
From: Roy Stogner <roystgnr@ic...>  20070116 17:48:32

On Tue, 16 Jan 2007, Tim Kr=F6ger wrote: > It would be nice if you at least tested whether this phenomenon also ha= ppens=20 > with your installation, in order to exclude version conflicts. It breaks for me unless I set pc_factor_zeropivot 0, but that happens with a lot of libMesh problems, where our penalty methods confuse PETSc. Try that option and see how it works for you. That may not be the only thing you need to tweak, though. For one thing, you may also want to try adding more linear solves. 250 isn't a general rule, it's just something that worked well for that one problem and default mesh size.  Roy 
From: <tim@ce...>  20070130 10:55:24
Attachments:
test.cpp

Dear Roy, On Tue, 16 Jan 2007, Roy Stogner wrote: > It breaks for me unless I set pc_factor_zeropivot 0, but that happens > with a lot of libMesh problems, where our penalty methods confuse > PETSc. Try that option and see how it works for you. Thank you very much; this seems to help. However, it seems to work only in oneprocessor mode. Please find attached a test program that solves the potential equation with Dirichlet boundaries in 2D. Since I didn't manage to create a simple test program that breaks unless pc_factor_zeropivot is used, I tested it the over way around: I ran the program * without this option (which works), and * with pc_factor_zeropivot 100 (which of course breaks). But when I use two processors, the program always works, and instead I get the message: WARNING! There are options you set that were not used! WARNING! could be spelling mistake, etc! Option left: name:pc_factor_zeropivot value: 100 Similar things happen if I try to hardcode this option into libMesh: In src/numerics/petsc_linear_solver.C at line 291, I inserted ierr = PCFactorSetZeroPivot(_pc,100); CHKERRABORT(libMesh::COMM_WORLD,ierr); Then, the program (of course) always breaks in oneprocessor mode. But with two processors, it still works fine, i.e. this option is again not respected. So what can I do to modify this tolerance also in multiprocessor mode? Best Regards, Tim 
From: Roy Stogner <roystgnr@ic...>  20070130 12:46:52

On Tue, 30 Jan 2007, Tim Kr=F6ger wrote: > On Tue, 16 Jan 2007, Roy Stogner wrote: > >> It breaks for me unless I set pc_factor_zeropivot 0, but that happens >> with a lot of libMesh problems, where our penalty methods confuse >> PETSc. Try that option and see how it works for you. > > Thank you very much; this seems to help. However, it seems to work onl= y in=20 > oneprocessor mode. In multiprocessor mode, ILU defaults to being a subpreconditioner, inside Block Jacobi I believe. Try sub_pc_factor_zeropivot 0.  Roy 
From: <tim@ce...>  20070201 10:17:17
Attachments:
test.cpp

Dear Roy, On Tue, 30 Jan 2007, Roy Stogner wrote: > On Tue, 30 Jan 2007, Tim Kr=F6ger wrote: > >> Thank you very much; this seems to help. However, it seems to work only= in=20 >> oneprocessor mode. > > In multiprocessor mode, ILU defaults to being a subpreconditioner, > inside Block Jacobi I believe. Try sub_pc_factor_zeropivot 0. Ah, thank you very much again. This makes it better. However, I=20 still face problems with my computation. I again attached you my test program, which differs from ex13.C only=20 in three points: (1) The grid is 100x100 (instead of 20x20), (2) the=20 linear tolarance is computed slightly different, and (3) a large fixed=20 obstacle is located in the computational domain (marked with "qqq" in=20 the source code). I tried to run this program on 1, 2, 4, 9, 10, 16, 18, 19, 20, 21, and=20 22 processors, with "sub_pc_factor_zeropivot 0" (on one processor=20 with "pc_factor_zeropivot 0"). In all these cases except with 20 processors, the program converged=20 towards a sensiblelooking solution (although the stopping criterion=20 was not reached, but that's not the point). For 20 processors,=20 however, the solution had strange oscillations. This was=20 reproducible. Could you please try whether this happens on your system as well? Do you have any idea what the reason is and what I can do to overcome=20 this problem? (In my real application, the nonconvergence just occures "sometimes",=20 without reproducibility, but much too often to be negligible. This=20 happens although the application does not contain any=20 nondeterministic part, and altough the number of processors is always=20 the same.) Best Regards, Tim 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20070201 12:14:44

As you point out, in parallel the default preconditioner is blockJacobi = ILU0. This preconditioner is particularly sensitive to nuances in the = partitioning... In elasticity problems if you have two plates connected = by a hinge, and the partitioning puts each plate on a processor with the = hinge as the interface (a perfectly natural thing to do...) the = processor local matrices look like two disjoint plates and the = preconditioner quality is severely degraded. (sorry for the aside into = structures, but I don't know how to conceptualize the problem in a = flow.) Anyway, as a matter of course I tend to use an additive schwarz = preconditioner with two levels of overlap. (one may be sufficient... I = have not tried to optimize the choice.) I have found that for my = applications this gives much more reliable convergence characteristics = on an arbitrary number of processors. Of course, your results may = vary... A typical commandline for my application code looks something = like this: mpirun np $np ./fins ./$case.in node_major_dofs \ ksp_gmres_restart 10 ksp_monitor \ sub_pc_type asm sub_pc_asm_overlap 2 \ sub_sub_pc_ilu_zeropivot 1.e20 \ log_summary \ > $case.log 2>&1 (note that is petsc2.2.x syntax) In your real application, is the mesh identical on the same number of = processors and this only occurs some of the time? =20 Original Message From: libmeshusersbounces@... = [mailto:libmeshusersbounces@...] On Behalf Of Tim = Kr=F6ger Sent: Thursday, February 01, 2007 4:11 AM To: Roy Stogner Cc: libmeshusers@... Subject: Re: [Libmeshusers] NavierStokes with obstacles Dear Roy, On Tue, 30 Jan 2007, Roy Stogner wrote: > On Tue, 30 Jan 2007, Tim Kr=F6ger wrote: > >> Thank you very much; this seems to help. However, it seems to work=20 >> only in oneprocessor mode. > > In multiprocessor mode, ILU defaults to being a subpreconditioner,=20 > inside Block Jacobi I believe. Try sub_pc_factor_zeropivot 0. Ah, thank you very much again. This makes it better. However, I still = face problems with my computation. I again attached you my test program, which differs from ex13.C only in = three points: (1) The grid is 100x100 (instead of 20x20), (2) the linear = tolarance is computed slightly different, and (3) a large fixed obstacle = is located in the computational domain (marked with "qqq" in the source = code). I tried to run this program on 1, 2, 4, 9, 10, 16, 18, 19, 20, 21, and 22 processors, with "sub_pc_factor_zeropivot 0" (on one processor with = "pc_factor_zeropivot 0"). In all these cases except with 20 processors, the program converged = towards a sensiblelooking solution (although the stopping criterion was = not reached, but that's not the point). For 20 processors, however, the = solution had strange oscillations. This was reproducible. Could you please try whether this happens on your system as well? Do you have any idea what the reason is and what I can do to overcome = this problem? (In my real application, the nonconvergence just occures "sometimes", = without reproducibility, but much too often to be negligible. This = happens although the application does not contain any nondeterministic = part, and altough the number of processors is always the same.) Best Regards, Tim 
From: <tim@ce...>  20070206 13:05:47

Dear Ben and Roy, On Thu, 1 Feb 2007, Kirk, Benjamin (JSCEG) wrote: > mpirun np $np ./fins ./$case.in node_major_dofs \ > ksp_gmres_restart 10 ksp_monitor \ > sub_pc_type asm sub_pc_asm_overlap 2 \ > sub_sub_pc_ilu_zeropivot 1.e20 \ > log_summary \ > > $case.log 2>&1 > > (note that is petsc2.2.x syntax) Thank you for your help. Indeed, this helps a lot. The code now runs much more stable (although I had to remove the "node_major_dofs" option). And if it breaks, this is now reproducible. A minor problem is now that I cannot use command line arguments with my programs (except for testing). The main reason is that the code will be embedded into a large software package where I have no access to the overall command line arguments. Hence, I have to do things with calls to PetSC functions. I tried to do this before the System::solve() call (by getting the KSP and PC context from the PetscLinearSolver class), but this failed since I have to call functions like PCBJacobiGetSubKSP(), which requires to call KSPSetUp() first, which requires to call KSPSetOperators() first, which I cannot do before matrix assembly. As a dirty hack, I implemented everything inside PetscLinearSolver::solve(). My idea to overcome this is to implement the possibility to set a callback function that is called in PetscLinearSolver::solve() immediately before the call to KSPSolve(). I will implement this and send you a patch, unless some of you has a better idea. > In your real application, is the mesh identical on the same number of processors and this only occurs some of the time? Yes, that's what happened before I applied your hint. Best Regards, Tim 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20070207 11:40:22

In the case that you do not have access to the command line you can = still 'pretend' you do to set arguments. The functions = 'PetscOptionSetValue' and 'PetscOptionClearValue' can be used. This is = useful because it prevents us from having to wrap *every* API call that = petsc has available. (yes, there are a lot.) For example, here is a = snippet from my app: if (Control::have_variable ("Parameters/matrix_free")) PetscOptionsSetValue("snes_mf_operator",PETSC_NULL); else PetscOptionsClearValue("snes_mf_operator"); It would be a straightforward thing to set in your input file something = like Petsc_ops "foo bar baz..." And then loop over the options and pass them to PETSc. Seems weird, but will this do the trick? =20 Original Message From: Tim Kr=F6ger [mailto:tim@...]=20 Sent: Tuesday, February 06, 2007 7:04 AM To: Kirk, Benjamin (JSCEG) Cc: Roy Stogner; libmeshusers@... Subject: RE: [Libmeshusers] NavierStokes with obstacles Dear Ben and Roy, On Thu, 1 Feb 2007, Kirk, Benjamin (JSCEG) wrote: > mpirun np $np ./fins ./$case.in node_major_dofs \ > ksp_gmres_restart 10 ksp_monitor \ > sub_pc_type asm sub_pc_asm_overlap 2 \ > sub_sub_pc_ilu_zeropivot 1.e20 \ > log_summary \ > > $case.log 2>&1 > > (note that is petsc2.2.x syntax) Thank you for your help. Indeed, this helps a lot. The code now runs = much more stable (although I had to remove the "node_major_dofs"=20 option). And if it breaks, this is now reproducible. A minor problem is now that I cannot use command line arguments with my = programs (except for testing). The main reason is that the code will be = embedded into a large software package where I have no access to the = overall command line arguments. Hence, I have to do things with calls = to PetSC functions. I tried to do this before the System::solve() call (by getting the KSP and PC context from the = PetscLinearSolver class), but this failed since I have to call functions = like PCBJacobiGetSubKSP(), which requires to call KSPSetUp() first, = which requires to call KSPSetOperators() first, which I cannot do before = matrix assembly. As a dirty hack, I implemented everything inside = PetscLinearSolver::solve(). My idea to overcome this is to implement = the possibility to set a callback function that is called in PetscLinearSolver::solve() immediately before the call to KSPSolve().=20 I will implement this and send you a patch, unless some of you has a = better idea. > In your real application, is the mesh identical on the same number of = processors and this only occurs some of the time? Yes, that's what happened before I applied your hint. Best Regards, Tim 
From: <tim@ce...>  20070208 13:45:53

Dear Ben, On Wed, 7 Feb 2007, Kirk, Benjamin (JSCEG) wrote: > In the case that you do not have access to the command line you can > still 'pretend' you do to set arguments. The functions > 'PetscOptionSetValue' and 'PetscOptionClearValue' can be used. Thank you very much. This solves the problem. (At least until the next problem arises...) Best Regards, Tim 
From: John Peterson <peterson@cf...>  20070116 08:32:15

Hi Tim, Upon closer look at your code I finally noticed this part: What does it do exactly=3F It can be quite tricky trying to work with two different "penalty" values, in case they accidentily end up in the same row. /* qqq begin */ if(fabs(elem>centroid()(0)0.5)<=3D0.4 && fabs(elem>centroid()(= 1)0.5)<=3D0.4) =09{ =09 const Real penalty =3D 1.e7; =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) =09 { =09 Kuu(i,i) +=3D penalty; =09 Kvv(i,i) +=3D penalty; =09 } =09} /* qqq end */ Tim Kr=F6ger writes: > Dear Roy and John, >=20 > thank you for your nearly simultaneous answer. (Did you use a=20 > stopwatch to get that right=3F) (: >=20 > Actually, constraining away the onedimensional null space did not=20= > cure the problem in my case. Please find attached a modified versio= n=20 > of my test program that includes John's code snippet. Again, all th= e=20 > (very little) places where the program differs from ex13 are marked=20= > with "qqq". The behaviour is the same as before, i.e. the program=20= > works well without obstacle (rough and fine grid), and it works well= =20 > with obstacle on the rough grid (20x20), but it crashes with obstacl= e=20 > on the fine grid (100x100). I pinned the pressure to 100.0 at node = 0,=20 > but pinning it to 0.0 instead does not change the problem. >=20 > It would be nice if you at least tested whether this phenomenon also= =20 > happens with your installation, in order to exclude version conflict= s. >=20 > Thank you very much in advance, >=20 > Tim > /* $Id: ex13.C,v 1.16 2006/10/23 18:46:06 jwpeterson Exp $ */ >=20 > /* The Next Great Finite Element Library. */ > /* Copyright (C) 2003 Benjamin S. Kirk */ >=20 > /* 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= . */ >=20 > /* 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. */ >=20 > /* You should have received a copy of the GNU Lesser General Public = */ > /* License along with this library; if not, write to the Free Softwa= re */ > /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 021111= 307 USA */ >=20 > // <h1>Example 13  Unsteady NavierStokes Equations  Unsteady Non= linear Systems of Equations</h1> > // > // This example shows how a simple, unsteady, nonlinear system of e= quations > // can be solved in parallel. The system of equations are the fami= liar > // NavierStokes equations for lowspeed incompressible fluid flow.= This > // example introduces the concept of the inner nonlinear loop for e= ach > // timestep, and requires a good deal of linear algebra numbercrun= ching > // at each step. If you have the General Mesh Viewer (GMV) install= ed, > // 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 i= f 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=5Ftimesteps and dt variables below. >=20 > // C++ include files that we need > #include <iostream> > #include <algorithm> > #include <math.h> >=20 > // Basic include file needed for the mesh functionality. > #include "libmesh.h" > #include "mesh.h" > #include "mesh=5Fgeneration.h" > #include "gmv=5Fio.h" > #include "equation=5Fsystems.h" > #include "fe.h" > #include "quadrature=5Fgauss.h" > #include "dof=5Fmap.h" > #include "sparse=5Fmatrix.h" > #include "numeric=5Fvector.h" > #include "dense=5Fmatrix.h" > #include "dense=5Fvector.h" > #include "linear=5Fimplicit=5Fsystem.h" > #include "transient=5Fsystem.h" > #include "perf=5Flog.h" > #include "boundary=5Finfo.h" >=20 > // Some (older) compilers do not offer full stream=20 > // functionality, OStringStream works around this. > #include "o=5Fstring=5Fstream.h" >=20 > // For systems of equations the \p DenseSubMatrix > // and \p DenseSubVector provide convenient ways for > // assembling the element matrix and vector on a > // componentbycomponent basis. > #include "dense=5Fsubmatrix.h" > #include "dense=5Fsubvector.h" >=20 > // The definition of a geometric element > #include "elem.h" >=20 > // Function prototype. This function will assemble the system > // matrix and righthandside. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname); >=20 > // The main program. > int main (int argc, char** argv) > { > // Initialize libMesh. > libMesh::init (argc, argv); > { =20 > // Set the dimensionality of the mesh =3D 2 > const unsigned int dim =3D 2; =20 > =20 > // Create a twodimensional mesh. > Mesh mesh (dim); > =20 > // Use the MeshTools::Generation mesh generator to create a unif= orm > // grid on the square [1,1]^D. We instruct the mesh generator > // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 > // elements in 3D. Building these higherorder elements allows > // us to use higherorder approximation, as in example 3. > MeshTools::Generation::build=5Fsquare (mesh, > =09=09=09=09=09 //qqq=09=09=09=09=09 20, 20, > =09=09=09=09=09 100, 100,//qqq > =09=09=09=09=09 0., 1., > =09=09=09=09=09 0., 1., > =09=09=09=09=09 QUAD9); > =20 > // Print information about the mesh to the screen. > mesh.print=5Finfo(); > =20 > // Create an equation systems object. > EquationSystems equation=5Fsystems (mesh); > =20 > // Declare the system and its variables. > { > // Creates a transient system named "NavierStokes" > TransientLinearImplicitSystem & system =3D=20 > =09equation=5Fsystems.add=5Fsystem<TransientLinearImplicitSystem> ("= NavierStokes"); > =20 > // Add the variables "u" & "v" to "NavierStokes". They > // will be approximated using secondorder approximation. > system.add=5Fvariable ("u", SECOND); > system.add=5Fvariable ("v", SECOND); >=20 > // Add the variable "p" to "NavierStokes". This will > // be approximated with a firstorder basis, > // providing an LBBstable pressurevelocity pair. > system.add=5Fvariable ("p", FIRST); >=20 > // Give the system a pointer to the matrix assembly > // function. > system.attach=5Fassemble=5Ffunction (assemble=5Fstokes); > =20 > // Initialize the data structures for the equation system. > equation=5Fsystems.init (); >=20 > equation=5Fsystems.parameters.set<unsigned int>("linear solver= maximum iterations") =3D 250; > equation=5Fsystems.parameters.set<Real> ("linear solver= tolerance") =3D std::sqrt(TOLERANCE); > =20 > // Prints information about the system to the screen. > equation=5Fsystems.print=5Finfo(); > } >=20 > // Create a performancelogging object for this example > PerfLog perf=5Flog("Example 13"); > =20 > // Now we begin the timestep loop to compute the timeaccurate > // solution of the equations. > const Real dt =3D 0.005; > Real time =3D 0.0; > const unsigned int n=5Ftimesteps =3D 15; >=20 > // The number of steps and the stopping criterion are also requi= red > // for the nonlinear iterations. > const unsigned int n=5Fnonlinear=5Fsteps =3D 15; > const Real nonlinear=5Ftolerance =3D 1.e2; > =20 > // Tell the system of equations what the timestep is by using > // the set=5Fparameter function. The matrix assembly routine ca= n > // then reference this parameter. > equation=5Fsystems.parameters.set<Real> ("dt") =3D dt; >=20 > // Get a reference to the Stokes system to use later. > TransientLinearImplicitSystem& navier=5Fstokes=5Fsystem =3D > =09 equation=5Fsystems.get=5Fsystem<TransientLinearImplicitSystem>(= "NavierStokes"); >=20 > // 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=5Fnonlinear=5Fsoln (navier=5Fstokes=5Fsystem.solution>cl= one()); >=20 > for (unsigned int t=5Fstep=3D0; t=5Fstep<n=5Ftimesteps; ++t=5Fst= ep) > { > =09// Incremenet the time counter, set the time and the > =09// time step size as parameters in the EquationSystem. > =09time +=3D dt; >=20 > =09// Let the system of equations know the current time. > =09// This might be necessary for a timedependent forcing > =09// function for example. > =09equation=5Fsystems.parameters.set<Real> ("time") =3D time; >=20 > =09// A pretty update message > =09std::cout << " Solving time step " << t=5Fstep << ", time =3D " <= < time << std::endl; >=20 > =09// Now we need to update the solution vector from the > =09// previous time step. This is done directly through > =09// the reference to the Stokes system. > =09*navier=5Fstokes=5Fsystem.old=5Flocal=5Fsolution =3D *navier=5Fst= okes=5Fsystem.current=5Flocal=5Fsolution; >=20 > =09// Now we begin the nonlinear loop > =09for (unsigned int l=3D0; l<n=5Fnonlinear=5Fsteps; ++l) > =09 { > =09 // Update the nonlinear solution. > =09 last=5Fnonlinear=5Fsoln>zero(); > =09 last=5Fnonlinear=5Fsoln>add(*navier=5Fstokes=5Fsystem.soluti= on); > =09 =20 > =09 // Assemble & solve the linear system. > =09 perf=5Flog.start=5Fevent("linear solve"); > =09 equation=5Fsystems.get=5Fsystem("NavierStokes").solve(); > =09 perf=5Flog.stop=5Fevent("linear solve"); >=20 > =09 // Compute the difference between this solution and the last > =09 // nonlinear iterate. > =09 last=5Fnonlinear=5Fsoln>add (1., *navier=5Fstokes=5Fsystem.= solution); >=20 > =09 // Close the vector before computing its norm > =09 last=5Fnonlinear=5Fsoln>close(); >=20 > =09 // Compute the l2 norm of the difference > =09 const Real norm=5Fdelta =3D last=5Fnonlinear=5Fsoln>l2=5Fnor= m(); >=20 > =09 // Print out convergence information for the linear and > =09 // nonlinear iterations. > =09 std::cout << "Linear solver converged at step: " > =09=09 << navier=5Fstokes=5Fsystem.n=5Flinear=5Fiterations() > =09=09 << ", final residual: " > =09=09 << navier=5Fstokes=5Fsystem.final=5Flinear=5Fresidual() > =09=09 << " Nonlinear convergence: u  u=5Fold =3D " > =09=09 << norm=5Fdelta > =09=09 << std::endl; >=20 > =09 // Terminate the solution iteration if the difference between= > =09 // this iteration and the last is sufficiently small. > =09 if (norm=5Fdelta < nonlinear=5Ftolerance) > =09 { > =09=09std::cout << " Nonlinear solver converged at step " > =09=09=09 << l > =09=09=09 << std::endl; > =09=09break; > =09 } > =09 =20 > =09 } // end nonlinear loop > =09 > =09// Write out every nth timestep to file. > =09const unsigned int write=5Finterval =3D 1; > =09 > =09if ((t=5Fstep+1)%write=5Finterval =3D=3D 0) > =09 { > =09 OStringStream file=5Fname; >=20 > =09 // We write the file name in the gmv autoread format. > =09 file=5Fname << "out.gmv."; > =09 OSSRealzeroright(file=5Fname,3,0, t=5Fstep + 1); > =09 =20 > =09 GMVIO(mesh).write=5Fequation=5Fsystems (file=5Fname.str(), > =09=09=09=09=09=09equation=5Fsystems); > =09 } > } // end timestep loop. > } > =20 > // All done. =20 > return libMesh::close (); > } >=20 >=20 >=20 >=20 >=20 >=20 > // The matrix assembly function to be called at each time step to > // prepare for the linear solve. > void assemble=5Fstokes (EquationSystems& es, > =09=09 const std::string& system=5Fname) > { > // It is a good idea to make sure we are assembling > // the proper system. > assert (system=5Fname =3D=3D "NavierStokes"); > =20 > // Get a constant reference to the mesh object. > const Mesh& mesh =3D es.get=5Fmesh(); > =20 > // The dimension that we are running > const unsigned int dim =3D mesh.mesh=5Fdimension(); > =20 > // Get a reference to the Stokes system object. > TransientLinearImplicitSystem & navier=5Fstokes=5Fsystem =3D > es.get=5Fsystem<TransientLinearImplicitSystem> ("NavierStokes")= ; >=20 > // Numeric ids corresponding to each variable in the system > const unsigned int u=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("u"); > const unsigned int v=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("v"); > const unsigned int p=5Fvar =3D navier=5Fstokes=5Fsystem.variable=5F= number ("p"); > =20 > // Get the Finite Element type for "u". Note this will be > // the same as the type for "v". > FEType fe=5Fvel=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Ftyp= e(u=5Fvar); > =20 > // Get the Finite Element type for "p". > FEType fe=5Fpres=5Ftype =3D navier=5Fstokes=5Fsystem.variable=5Fty= pe(p=5Fvar); >=20 > // Build a Finite Element object of the specified type for > // the velocity variables. > AutoPtr<FEBase> fe=5Fvel (FEBase::build(dim, fe=5Fvel=5Ftype)); > =20 > // Build a Finite Element object of the specified type for > // the pressure variables. > AutoPtr<FEBase> fe=5Fpres (FEBase::build(dim, fe=5Fpres=5Ftype)); > =20 > // A Gauss quadrature rule for numerical integration. > // Let the \p FEType object decide what order rule is appropriate.= > QGauss qrule (dim, fe=5Fvel=5Ftype.default=5Fquadrature=5Forder())= ; >=20 > // Tell the finite element objects to use our quadrature rule. > fe=5Fvel>attach=5Fquadrature=5Frule (&qrule); > fe=5Fpres>attach=5Fquadrature=5Frule (&qrule); > =20 > // Here we define some references to cellspecific data that > // will be used to assemble the linear system. > // > // The element Jacobian * quadrature weight at each integration po= int. =20 > const std::vector<Real>& JxW =3D fe=5Fvel>get=5FJxW(); >=20 > // The element shape functions evaluated at the quadrature points.= > const std::vector<std::vector<Real> >& phi =3D fe=5Fvel>get=5Fphi= (); >=20 > // The element shape function gradients for the velocity > // variables evaluated at the quadrature points. > const std::vector<std::vector<RealGradient> >& dphi =3D fe=5Fvel>= get=5Fdphi(); >=20 > // The element shape functions for the pressure variable > // evaluated at the quadrature points. > const std::vector<std::vector<Real> >& psi =3D fe=5Fpres>get=5Fph= i(); >=20 > // The value of the linear shape function gradients at the quadrat= ure points > // const std::vector<std::vector<RealGradient> >& dpsi =3D fe=5Fpr= es>get=5Fdphi(); > =20 > // A reference to the \p DofMap object for this system. The \p Do= fMap > // object handles the index translation from node and element numb= ers > // to degree of freedom numbers. We will talk more about the \p D= ofMap > // in future examples. > const DofMap & dof=5Fmap =3D navier=5Fstokes=5Fsystem.get=5Fdof=5F= map(); >=20 > // Define data structures to contain the element matrix > // and righthandside vector contribution. Following > // basic finite element terminology we will denote these > // "Ke" and "Fe". > DenseMatrix<Number> Ke; > DenseVector<Number> Fe; >=20 > DenseSubMatrix<Number> > Kuu(Ke), Kuv(Ke), Kup(Ke), > Kvu(Ke), Kvv(Ke), Kvp(Ke), > Kpu(Ke), Kpv(Ke), Kpp(Ke); >=20 > DenseSubVector<Number> > Fu(Fe), > Fv(Fe), > Fp(Fe); >=20 > // 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<unsigned int> dof=5Findices; > std::vector<unsigned int> dof=5Findices=5Fu; > std::vector<unsigned int> dof=5Findices=5Fv; > std::vector<unsigned int> dof=5Findices=5Fp; >=20 > // Find out what the timestep size parameter is from the system, a= nd > // the value of theta for the theta method. We use implicit Euler= (theta=3D1) > // for this simulation even though it is only firstorder accurate= in time. > // The reason for this decision is that the secondorder CrankNic= olson > // method is notoriously oscillatory for problems with discontinuo= us > // initial data such as the liddriven cavity. Therefore, > // we sacrifice accuracy in time for stability, but since the solu= tion > // reaches steady state relatively quickly we can afford to take s= mall > // timesteps. If you monitor the initial nonlinear residual for t= his > // simulation, you should see that it is monotonically decreasing = in time. > const Real dt =3D es.parameters.get<Real>("dt"); > // const Real time =3D es.parameters.get<Real>("time"); > const Real theta =3D 1.; > =20 > // Now we will loop over all the elements in the mesh that > // live on the local processor. We will compute the element > // matrix and righthandside contribution. Since the mesh > // will be refined we want to only consider the ACTIVE elements, > // hence we use a variant of the \p active=5Felem=5Fiterator. > // const=5Factive=5Flocal=5Felem=5Fiterator el (mesh.ele= ments=5Fbegin()); > // const const=5Factive=5Flocal=5Felem=5Fiterator end=5Fel (mesh.e= lements=5Fend()); >=20 > MeshBase::const=5Felement=5Fiterator el =3D mesh.active=5F= local=5Felements=5Fbegin(); > const MeshBase::const=5Felement=5Fiterator end=5Fel =3D mesh.activ= e=5Flocal=5Felements=5Fend();=20 > =20 > for ( ; el !=3D end=5Fel; ++el) > { =20 > // Store a pointer to the element we are currently > // working on. This allows for nicer syntax later. > const Elem* elem =3D *el; > =20 > // Get the degree of freedom indices for the > // current element. These define where in the global > // matrix and righthandside this element will > // contribute to. > dof=5Fmap.dof=5Findices (elem, dof=5Findices); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fu, u=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fv, v=5Fvar); > dof=5Fmap.dof=5Findices (elem, dof=5Findices=5Fp, p=5Fvar); >=20 > const unsigned int n=5Fdofs =3D dof=5Findices.size(); > const unsigned int n=5Fu=5Fdofs =3D dof=5Findices=5Fu.size();=20= > const unsigned int n=5Fv=5Fdofs =3D dof=5Findices=5Fv.size();=20= > const unsigned int n=5Fp=5Fdofs =3D dof=5Findices=5Fp.size(); > =20 > // Compute the elementspecific data for the current > // element. This involves computing the location of the > // quadrature points (q=5Fpoint) and the shape functions > // (phi, dphi) for the current element. > fe=5Fvel>reinit (elem); > fe=5Fpres>reinit (elem); >=20 > // Zero the element matrix and righthand side 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=5Fdofs, n=5Fdofs); > Fe.resize (n=5Fdofs); >=20 > // Reposition the submatrices... The idea is this: > // > //     > //  Kuu Kuv Kup   Fu  > // Ke =3D  Kvu Kvv Kvp ; Fe =3D  Fv  > //  Kpu Kpv Kpp   Fp  > //     > // > // The \p DenseSubMatrix.repostition () member takes the > // (row=5Foffset, column=5Foffset, row=5Fsize, column=5Fsize).= > //=20 > // Similarly, the \p DenseSubVector.reposition () member > // takes the (row=5Foffset, row=5Fsize) > Kuu.reposition (u=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fu=5Fdofs); > Kuv.reposition (u=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fv=5Fdofs); > Kup.reposition (u=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= u=5Fdofs, n=5Fp=5Fdofs); > =20 > Kvu.reposition (v=5Fvar*n=5Fv=5Fdofs, u=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fu=5Fdofs); > Kvv.reposition (v=5Fvar*n=5Fv=5Fdofs, v=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fv=5Fdofs); > Kvp.reposition (v=5Fvar*n=5Fv=5Fdofs, p=5Fvar*n=5Fv=5Fdofs, n=5F= v=5Fdofs, n=5Fp=5Fdofs); > =20 > Kpu.reposition (p=5Fvar*n=5Fu=5Fdofs, u=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fu=5Fdofs); > Kpv.reposition (p=5Fvar*n=5Fu=5Fdofs, v=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fv=5Fdofs); > Kpp.reposition (p=5Fvar*n=5Fu=5Fdofs, p=5Fvar*n=5Fu=5Fdofs, n=5F= p=5Fdofs, n=5Fp=5Fdofs); >=20 > Fu.reposition (u=5Fvar*n=5Fu=5Fdofs, n=5Fu=5Fdofs); > Fv.reposition (v=5Fvar*n=5Fu=5Fdofs, n=5Fv=5Fdofs); > Fp.reposition (p=5Fvar*n=5Fu=5Fdofs, n=5Fp=5Fdofs); >=20 > // Now we will build the element matrix and righthandside. > // Constructing the RHS requires the solution and its > // gradient from the previous timestep. This must be > // calculated at each quadrature point by summing the > // solution degreeoffreedom values by the appropriate > // weight functions. > for (unsigned int qp=3D0; qp<qrule.n=5Fpoints(); qp++) > =09{ > =09 // Values to hold the solution & its gradient at the previous t= imestep. > =09 Number u =3D 0., u=5Fold =3D 0.; > =09 Number v =3D 0., v=5Fold =3D 0.; > =09 Number p=5Fold =3D 0.; > =09 Gradient grad=5Fu, grad=5Fu=5Fold; > =09 Gradient grad=5Fv, grad=5Fv=5Fold; > =09 =20 > =09 // Compute the velocity & its gradient from the previous timest= ep > =09 // and the old Newton iterate. > =09 for (unsigned int l=3D0; l<n=5Fu=5Fdofs; l++) > =09 { > =09 // From the old timestep: > =09 u=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fu[l]); > =09 v=5Fold +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fv[l]); > =09 grad=5Fu=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv=5Fold.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5F= system.old=5Fsolution (dof=5Findices=5Fv[l])); >=20 > =09 // From the previous Newton iterate: > =09 u +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fu[l]);=20 > =09 v +=3D phi[l][qp]*navier=5Fstokes=5Fsystem.current=5Fsoluti= on (dof=5Findices=5Fv[l]); > =09 grad=5Fu.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fu[l])); > =09 grad=5Fv.add=5Fscaled (dphi[l][qp],navier=5Fstokes=5Fsystem= .current=5Fsolution (dof=5Findices=5Fv[l])); > =09 } >=20 > =09 // Compute the old pressure value at this quadrature point. > =09 for (unsigned int l=3D0; l<n=5Fp=5Fdofs; l++) > =09 { > =09 p=5Fold +=3D psi[l][qp]*navier=5Fstokes=5Fsystem.old=5Fsolu= tion (dof=5Findices=5Fp[l]); > =09 } >=20 > =09 // Definitions for convenience. It is sometimes simpler to do = a > =09 // dot product if you have the full vector at your disposal. > =09 const NumberVectorValue U=5Fold (u=5Fold, v=5Fold); > =09 const NumberVectorValue U (u, v); > =09 const Number u=5Fx =3D grad=5Fu(0); > =09 const Number u=5Fy =3D grad=5Fu(1); > =09 const Number v=5Fx =3D grad=5Fv(0); > =09 const Number v=5Fy =3D grad=5Fv(1); > =09 =20 > =09 // First, an iloop over the velocity degrees of freedom. > =09 // We know that n=5Fu=5Fdofs =3D=3D n=5Fv=5Fdofs so we can comp= ute contributions > =09 // for both at the same time. > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Fu(i) +=3D JxW[qp]*(u=5Fold*phi[i][qp]  = // massmatrix term=20 > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fu=5Fold)*phi[i][qp] + // c= onvection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](0)  // press= ure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fu=5Fold*dphi[i][qp]) + // dif= fusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fu)*phi[i][qp]); // Newto= n term >=20 > =09=09 > =09 Fv(i) +=3D JxW[qp]*(v=5Fold*phi[i][qp]  = // massmatrix term > =09=09=09=09(1.theta)*dt*(U=5Fold*grad=5Fv=5Fold)*phi[i][qp] + // = convection term > =09=09=09=09(1.theta)*dt*p=5Fold*dphi[i][qp](1)  // pres= sure term on rhs > =09=09=09=09(1.theta)*dt*(grad=5Fv=5Fold*dphi[i][qp]) + // di= ffusion term on rhs > =09=09=09=09theta*dt*(U*grad=5Fv)*phi[i][qp]); // Newt= on term > =09 =20 >=20 > =09 // Note that the Fp block is identically zero unless we are= using > =09 // some kind of artificial compressibility scheme... >=20 > =09 // Matrix contributions for the uu and vv couplings. > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kuu(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*u=5Fx*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kuv(i,j) +=3D JxW[qp]*theta*dt*u=5Fy*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09 =20 > =09=09 Kvv(i,j) +=3D JxW[qp]*(phi[i][qp]*phi[j][qp] + = // mass matrix term > =09=09=09=09 theta*dt*(dphi[i][qp]*dphi[j][qp]) + // diffusi= on term > =09=09=09=09 theta*dt*(U*dphi[j][qp])*phi[i][qp] + // convect= ion term > =09=09=09=09 theta*dt*v=5Fy*phi[i][qp]*phi[j][qp]); // Newto= n term >=20 > =09=09 Kvu(i,j) +=3D JxW[qp]*theta*dt*v=5Fx*phi[i][qp]*phi[j][qp]; = // Newton term > =09=09} >=20 > =09 // Matrix contributions for the up and vp couplings. > =09 for (unsigned int j=3D0; j<n=5Fp=5Fdofs; j++) > =09=09{ > =09=09 Kup(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](0));= > =09=09 Kvp(i,j) +=3D JxW[qp]*(theta*dt*psi[j][qp]*dphi[i][qp](1));= > =09=09} > =09 } >=20 > =09 // Now an iloop over the pressure degrees of freedom. This co= de computes > =09 // the matrix entries due to the continuity equation. Note: To= maintain a > =09 // symmetric matrix, we may (or may not) multiply the continuit= y equation by > =09 // negative one. Here we do not. > =09 for (unsigned int i=3D0; i<n=5Fp=5Fdofs; i++) > =09 { > =09 for (unsigned int j=3D0; j<n=5Fu=5Fdofs; j++) > =09=09{ > =09=09 Kpu(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](0); > =09=09 Kpv(i,j) +=3D JxW[qp]*psi[i][qp]*dphi[j][qp](1); > =09=09} > =09 } > =09} // end of the quadrature point qploop >=20 > =20 > // At this point the interior element integration has > // been completed. However, we have not yet addressed > // boundary conditions. For this example we will only > // consider simple Dirichlet boundary conditions imposed > // via the penalty method. The penalty method used here > // is equivalent (for Lagrange basis functions) to lumping > // the matrix resulting from the L2 projection penalty > // approach introduced in example 3. > { > =09// The following loops over the sides of the element. > =09// If the element has no neighbor on a side then that > =09// side MUST live on a boundary of the domain. > =09for (unsigned int s=3D0; s<elem>n=5Fsides(); s++) > =09 if (elem>neighbor(s) =3D=3D NULL) > =09 { > =09 // Get the boundary ID for side 's'. > =09 // These are set internally by build=5Fsquare(). > =09 // 0=3Dbottom > =09 // 1=3Dright > =09 // 2=3Dtop > =09 // 3=3Dleft > =09 short int bc=5Fid =3D mesh.boundary=5Finfo>boundary=5Fid (= elem,s); > =09 if (bc=5Fid=3D=3DBoundaryInfo::invalid=5Fid) > =09=09 error(); >=20 > =09 =20 > =09 AutoPtr<Elem> side (elem>build=5Fside(s)); > =09 =09 =20 > =09 // Loop over the nodes on the side. > =09 for (unsigned int ns=3D0; ns<side>n=5Fnodes(); ns++) > =09=09{ > =09=09 // The penalty value. \f$ \frac{1}{\epsilon \f$ > =09=09 const Real penalty =3D 1.e10; > =09=09 =20 > =09=09 // Get the boundary values. > =09=09 =20 > =09=09 // Set u =3D 1 on the top boundary, 0 everywhere else > =09=09 const Real u=5Fvalue =3D (bc=5Fid=3D=3D2) =3F 1. : 0.; > =09=09 =20 > =09=09 // Set v =3D 0 everywhere > =09=09 const Real v=5Fvalue =3D 0.; > =09=09 =20 > =09=09 // Find the node on the element matching this node on > =09=09 // the side. That defined where in the element matrix > =09=09 // the boundary condition will be applied. > =09=09 for (unsigned int n=3D0; n<elem>n=5Fnodes(); n++) > =09=09 if (elem>node(n) =3D=3D side>node(ns)) > =09=09 { > =09=09=09// Matrix contribution. > =09=09=09Kuu(n,n) +=3D penalty; > =09=09=09Kvv(n,n) +=3D penalty; > =09=09 =09=09 =20 > =09=09=09// Righthandside contribution. > =09=09=09Fu(n) +=3D penalty*u=5Fvalue; > =09=09=09Fv(n) +=3D penalty*v=5Fvalue; > =09=09 } > =09=09} // end face node loop=09 =20 > =09 } // end if (elem>neighbor(side) =3D=3D NULL) >=20 > =09/* qqq begin */ >=20 > =09for (unsigned int c=3D0; c<elem>n=5Fnodes(); c++) > =09 if (elem>node(c) =3D=3D 0) > =09 { > =09 const Real penalty =3D 1.e10; > =09 Kpp(c,c) +=3D penalty; > =09 Fp(c) +=3D penalty*100.0; > =09 } > =09 > =09/* qqq end */ >=20 > } // end boundary condition section=09 =20 > =20 > /* qqq begin */ > if(fabs(elem>centroid()(0)0.5)<=3D0.4 && fabs(elem>centroid= ()(1)0.5)<=3D0.4) > =09{ > =09 const Real penalty =3D 1.e7; > =09 for (unsigned int i=3D0; i<n=5Fu=5Fdofs; i++) > =09 { > =09 Kuu(i,i) +=3D penalty; > =09 Kvv(i,i) +=3D penalty; > =09 } > =09} > /* qqq end */ >=20 > // The element matrix and righthandside are now built > // for this element. Add them to the global matrix and > // righthandside vector. The \p PetscMatrix::add=5Fmatrix()= > // and \p PetscVector::add=5Fvector() members do this for us. > navier=5Fstokes=5Fsystem.matrix>add=5Fmatrix (Ke, dof=5Findic= es); > navier=5Fstokes=5Fsystem.rhs>add=5Fvector (Fe, dof=5Findic= es); > } // end of element loop > =20 > // That's it. > return; > } 
From: <tim@ce...>  20070116 13:23:26

Dear John, On Tue, 16 Jan 2007, John Peterson wrote: > Upon closer look at your code I finally noticed this part: > What does it do exactly? Well, it's the obstacle. I.e., in this case, a (large) square area is considered to be solid and fixed, i.e. u=0 in this area. The idea is that this acts as a noslip boundary condition at the boundary of the obstacle. Within the obstacle, I'm not interested in the NavierStokes quantities anyway. However, removing these cells from the grid might be problematic in my real application for two reasons: First, I have a second equation that has to be solved in the whole domain, including the obstacles. Second, the size, number, and shape of obstacles changes in time. > It can be quite tricky trying to > work with two different "penalty" values, in case they accidentily > end up in the same row. In my opinion, they can't happen to be in the same row in this case (except for extremly rough grids) since the obstacle is away from the boundary of the domain. Also, the behaviour of the program does not change if 1e7 is used for *all* penalties. However, if 1e10 is used for all penalties, then the computation with obstacle fails already on the 20x20 grid. Best Regards, Tim 