Hi!
We ran into an issue with our framework here at INL when forming full
jacobian matrix as a preconditioner with AMR. What we hit in our problem
was the different size of the dense matrix (local jacobian contributions)
and the dof indices array. I traced the problem down to the level of
DofMap::constrain_element_matrix. I found out that this method can change
the dof indices array (actually the method that changes the dof indices
array is DofMap::build_constraint_matrix.) and that triggered an assert in
libMesh (attached is a modified ex19 that exhibits this behavior). This
only happens when there are constrains involved. I quess that no one
bumped into this since people usually call this constrain method just
once, but we call it multiple times and trying to reuse those arrays.
So my question: Is this known behavior? And is this used somewhere else in
libMesh, like we call that constrain_element_matrix and we use this
modified array for something later on?
Thanks,
--
David Andrs
/* $Id: ex4.C 2501 2007-11-20 02:33:29Z benkirk $ */
/* The Next Great Finite Element Library. */
/* Copyright (C) 2003 Benjamin S. Kirk */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */
/* License as published by the Free Software Foundation; either */
/* version 2.1 of the License, or (at your option) any later version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
/* Lesser General Public License for more details. */
/* You should have received a copy of the GNU Lesser General Public */
/* License along with this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
// <h1>Example 19 - Solving the 2D Young Laplace Problem using nonlinear
solvers</h1>
//
// This example shows how to use the NonlinearImplicitSystem class
// to efficiently solve nonlinear problems in parallel.
//
// In nonlinear systems, we aim at finding x that satisfy R(x) = 0.
// In nonlinear finite element analysis, the residual is typically
// of the form R(x) = K(x)*x - f, with K(x) the system matrix and f
// the "right-hand-side". The NonlinearImplicitSystem class expects
// two callback functions to compute the residual R and its Jacobian
// for the Newton iterations. Here, we just approximate
// the true Jacobian by K(x).
//
// You can turn on preconditining of the matrix free system using the
// jacobian by passing "-pre" on the command line. Currently this only
// work with Petsc so this isn't used by using "make run"
//
// This example also runs with the experimental Trilinos NOX solvers by
specifying
// the --use-trilinos command line argument.
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <cmath>
// Various include files needed for the mesh & solver functionality.
#include "libmesh.h"
#include "mesh.h"
#include "mesh_refinement.h"
#include "exodusII_io.h"
#include "equation_systems.h"
#include "fe.h"
#include "quadrature_gauss.h"
#include "dof_map.h"
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "elem.h"
#include "string_to_enum.h"
#include "getpot.h"
#include "dense_submatrix.h"
#include "dense_subvector.h"
#include "mesh_generation.h"
#include "coupling_matrix.h"
#include "boundary_info.h"
#include "mesh_refinement.h"
#include "error_vector.h"
#include "kelly_error_estimator.h"
// The nonlinear solver and system we will be using
#include "nonlinear_solver.h"
#include "nonlinear_implicit_system.h"
// Necessary for programmatically setting petsc options
#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;
// Let-s define the physical parameters of the equation
const Real kappa = 1.;
const Real sigma = 0.2;
// 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
NonlinearImplicitSystem& system =
es.get_system<NonlinearImplicitSystem>("Laplace-Young");
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
// 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();
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = dof_map.variable_type(0);
// Build a Finite Element object of the specified type. Since the
// \p FEBase::build() member dynamically creates memory we will
// store the object as an \p AutoPtr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
// A 5th order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, FIFTH);
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
// Here we define some references to cell-specific data that
// will be used to assemble the linear system.
// We begin with the element Jacobian * quadrature weight at each
// integration point.
const std::vector<Real>& JxW = fe->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
// The element shape function gradients evaluated at the quadrature
// points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
// Define data structures to contain the Jacobian element matrix.
// Following basic finite element terminology we will denote these
// "Ke". More detail is in example 3.
// DenseMatrix<Number> Ke;
// DenseMatrix<Number> Kuu(Ke), Kvv(Ke), Kuv(Ke);
DenseMatrix<Number> Kuu, Kvv, Kuv;
// 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_indices;
std::vector<unsigned int> dof_indices_u;
std::vector<unsigned int> dof_indices_v;
// 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)
{
// Store a pointer to the element we are currently
// working on. This allows for nicer syntax later.
const Elem* elem = *el;
// 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);
dof_map.dof_indices (elem, dof_indices_u, u_var);
dof_map.dof_indices (elem, dof_indices_v, v_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();
// 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);
// 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 (dof_indices.size(),
// dof_indices.size());
// Kuu.reposition( 0, 0, n_u_dofs, n_u_dofs);
// Kuv.reposition(n_u_dofs, 0, n_u_dofs, n_v_dofs);
// Kvv.reposition(n_u_dofs, n_u_dofs, n_v_dofs, n_v_dofs);
Kuu.resize(n_u_dofs, n_u_dofs);
Kuv.resize(n_u_dofs, n_v_dofs);
Kvv.resize(n_v_dofs, n_v_dofs);
// Now we will build the element Jacobian. This involves
// a double loop to integrate the test funcions (i) against
// the trial functions (j). Note that the Jacobian depends
// on the current solution x, which we access using the soln
// vector.
//
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Gradient grad_u;
Gradient grad_v;
for (unsigned int i=0; i<phi.size(); i++)
{
grad_u += dphi[i][qp]*soln(dof_indices_u[i]);
grad_v += dphi[i][qp]*soln(dof_indices_v[i]);
}
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
}
}
// dof_map.constrain_element_matrix (Ke, dof_indices);
dof_map.constrain_element_matrix (Kuu, dof_indices_u, dof_indices_u);
dof_map.constrain_element_matrix (Kuv, dof_indices_u, dof_indices_v);
dof_map.constrain_element_matrix (Kvv, dof_indices_v, dof_indices_v);
jacobian.add_matrix (Kuu, dof_indices_u, dof_indices_u);
jacobian.add_matrix (Kuv, dof_indices_u, dof_indices_v);
jacobian.add_matrix (Kvv, dof_indices_v, dof_indices_v);
// Add the element matrix to the system Jacobian.
// jacobian.add_matrix (Ke, dof_indices);
}
jacobian.close();
std::vector<int> rows;
std::vector<unsigned int> nl;
std::vector<short int> il;
mesh.boundary_info->build_node_list(nl, il);
for (int ii = 0; ii < nl.size(); ii++)
{
const Node & nd = mesh.node(nl[ii]);
if (il[ii] == 1 || il[ii] == 3)
{
unsigned int idx = nd.dof_number(system.number(), u_var, 0);
rows.push_back(idx);
}
}
jacobian.zero_rows(rows, 1);
jacobian.close();
// 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 (dim == 2);
// Get a reference to the NonlinearImplicitSystem we are solving
NonlinearImplicitSystem& system =
es.get_system<NonlinearImplicitSystem>("Laplace-Young");
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
// 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();
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = dof_map.variable_type(0);
// Build a Finite Element object of the specified type. Since the
// \p FEBase::build() member dynamically creates memory we will
// store the object as an \p AutoPtr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
// A 5th order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, FIFTH);
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
// Declare a special finite element object for
// boundary integration.
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
// Boundary integration requires one quadraure rule,
// with dimensionality one less than the dimensionality
// of the element.
QGauss qface(dim-1, FIFTH);
// Tell the finte element object to use our
// quadrature rule.
fe_face->attach_quadrature_rule (&qface);
// Here we define some references to cell-specific data that
// will be used to assemble the linear system.
// We begin with the element Jacobian * quadrature weight at each
// integration point.
const std::vector<Real>& JxW = fe->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
// The element shape function gradients evaluated at the quadrature
// points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
// Define data structures to contain the resdual contributions
DenseVector<Number> Re;
DenseSubVector<Number> Re_u(Re), Re_v(Re);
// 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_indices;
std::vector<unsigned int> dof_indices_u, dof_indices_v;
// Now we will loop over all the active elements in the mesh which
// are local to this processor.
// We will compute the element residual.
residual.zero();
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;
// 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);
dof_map.dof_indices (elem, dof_indices_u, u_var);
dof_map.dof_indices (elem, dof_indices_v, v_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();
// 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);
// 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).
Re.resize (dof_indices.size());
Re_u.reposition(u_var * n_u_dofs, n_u_dofs);
Re_v.reposition(v_var * n_u_dofs, n_v_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.
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Number u = 0;
Gradient grad_u;
Number v = 0;
Gradient grad_v;
for (unsigned int j=0; j<phi.size(); j++)
{
u += phi[j][qp]*soln(dof_indices_u[j]);
grad_u += dphi[j][qp]*soln(dof_indices_u[j]);
v += phi[j][qp]*soln(dof_indices_v[j]);
grad_v += dphi[j][qp]*soln(dof_indices_v[j]);
}
for (unsigned int i=0; i<phi.size(); i++)
{
Re_u(i) += JxW[qp]*(dphi[i][qp]*grad_u);
Re_v(i) += JxW[qp]*(dphi[i][qp]*grad_v);
}
}
// At this point the interior element integration has
// been completed. However, we have not yet addressed
// boundary conditions.
// // 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)
// {
// // The value of the shape functions at the quadrature
// // points.
// const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
//
// // The Jacobian * Quadrature Weight at the quadrature
// // points on the face.
// const std::vector<Real>& JxW_face = fe_face->get_JxW();
//
// // Compute the shape function values on the element face.
// fe_face->reinit(elem, side);
//
// // 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++)
// Re_u(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
// }
// }
dof_map.constrain_element_vector (Re, dof_indices);
residual.add_vector (Re, dof_indices);
}
residual.close();
std::vector<unsigned int> nl;
std::vector<short int> il;
mesh.boundary_info->build_node_list(nl, il);
for (int ii = 0; ii < nl.size(); ii++)
{
const Node & nd = mesh.node(nl[ii]);
if (il[ii] == 1)
{
unsigned int idx = nd.dof_number(system.number(), u_var, 0);
residual.set(idx, soln(idx) - 1.);
}
else if (il[ii] == 3)
{
unsigned int idx = nd.dof_number(system.number(), u_var, 0);
residual.set(idx, soln(idx) - 2.);
}
}
residual.close();
// That's it.
}
// Begin the main program.
int main (int argc, char** argv)
{
// Initialize libMesh and any dependent libaries, like in example 2.
LibMeshInit init (argc, argv);
#if !defined(LIBMESH_HAVE_PETSC) && !defined(LIBMESH_HAVE_TRILINOS)
if (libMesh::processor_id() == 0)
std::cerr << "ERROR: This example requires libMesh to be\n"
<< "compiled with nonlinear solver support from\n"
<< "PETSc or Trilinos!"
<< std::endl;
return 0;
#endif
#ifndef LIBMESH_ENABLE_AMR
if (libMesh::processor_id() == 0)
std::cerr << "ERROR: This example requires libMesh to be\n"
<< "compiled with AMR support!"
<< std::endl;
return 0;
#else
// Create a GetPot object to parse the command line
GetPot command_line (argc, argv);
// Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
// Create a mesh from file.
Mesh mesh;
MeshTools::Generation::build_square(mesh,
2, 2,
0, 1, 0, 1, QUAD4);
mesh.boundary_info->build_node_list_from_side_list();
// Print information about the mesh to the screen.
mesh.print_info();
// Create an equation systems object.
EquationSystems equation_systems (mesh);
_equation_system = &equation_systems;
// Declare the system and its variables.
// Creates a system named "Laplace-Young"
NonlinearImplicitSystem& system =
equation_systems.add_system<NonlinearImplicitSystem> ("Laplace-Young");
// 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-12;
equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum
iterations") = 50;
// Adds the variable "u" to "Laplace-Young". "u"
// will be approximated using second-order approximation.
system.add_variable("u", FIRST, LAGRANGE);
system.add_variable("v", FIRST, LAGRANGE);
CouplingMatrix cm(2);
cm(0, 0) = 1;
cm(1, 1) = 1;
cm(0, 1) = 1;
system.get_dof_map()._dof_coupling = &cm;
// 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();
MeshRefinement mesh_refinement (mesh);
const unsigned int max_r_steps = 5;
for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
{
// Solve the system "Laplace-Young", print the number of iterations
// and final residual
equation_systems.get_system("Laplace-Young").solve();
// Print out final convergence information. This duplicates some
// output from during the solve itself, but demonstrates another way
// to get this information after the solve is complete.
std::cout << "Laplace-Young system solved at nonlinear iteration "
<< system.n_nonlinear_iterations()
<< " , final nonlinear residual norm: "
<< system.final_nonlinear_residual()
<< std::endl;
if (r_step+1 != max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
ErrorVector error;
KellyErrorEstimator error_estimator;
error_estimator.estimate_error (system, error);
mesh_refinement.refine_fraction() = 0.01;
mesh_refinement.coarsen_fraction() = 0.0;
mesh_refinement.max_h_level() = 5;
mesh_refinement.flag_elements_by_error_fraction (error);
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit ();
mesh.boundary_info->build_node_list_from_side_list();
}
}
#ifdef LIBMESH_HAVE_EXODUS_API
// After solving the system write the solution
ExodusII_IO (mesh).write_equation_systems ("out.e",
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
#endif // #ifndef LIBMESH_ENABLE_AMR
// All done.
return 0;
}
------------------------------------------------------------------------------
All the data continuously generated in your IT infrastructure contains a
definitive record of customers, application performance, security
threats, fraudulent activity and more. Splunk takes this data and makes
sense of it. Business sense. IT sense. Common sense.
http://p.sf.net/sfu/splunk-d2dcopy1
_______________________________________________
Libmesh-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
