Roy Stogner <[email protected]> wrote on 06/08/2011 01:28:03 PM:
> > On Wed, 8 Jun 2011, David Andrs wrote: > > > This was actually a serial run with adaptivity and this bug poped out > > after first refinement when calling EquationSystems::reinit(). > > > > I'll try to replicate this behavior with pure libMesh calls and send it to > > you, it might be worth to include it into regression tests. > > Please do! If you mean the bug got triggered by refinement of a > SerialMesh in a parallel run, AFAIK that shouldn't have happened. > > And if you really mean the bug got triggered by refinement in a serial > run, that *definitely* shouldn't have happened. So it took me a while to prepare the test case, but it uncovered the source of the problem. The whole thing is caused by having two systems in the EquationSystems class - one system with some variables in it, but the other with NO variables at all. This caused the older and old vectors to be empty and therefore first_local_dof and last_local_dof were 0 and therefore the fail in the localize. Source code and a little mesh file are attached. The code is hackish, basically copy'n'paste to reproduce the issue. Run it with -snes_mf_operator command line option. My libMesh installation is configured with serial mesh and the run is also serial. Hope this helps, -- David
square.e
Description: Binary data
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <cmath>
// Basic include file needed for the mesh functionality.
#include "libmesh.h"
#include "mesh.h"
#include "mesh_refinement.h"
#include "gmv_io.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 "getpot.h"
#include "o_string_stream.h"
// This example will solve a linear transient system,
// so we need to include the \p TransientLinearImplicitSystem definition.
#include "transient_system.h"
#include "nonlinear_implicit_system.h"
#include "vector_value.h"
#include "nonlinear_solver.h"
#include "boundary_info.h"
// To refine the mesh we need an \p ErrorEstimator
// object to figure out which elements to refine.
#include "error_vector.h"
#include "kelly_error_estimator.h"
// The definition of a geometric element
#include "elem.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
const Real dt = 0.1;
RealVectorValue vel(10, 1);
void compute_jacobian (const NumericVector<Number>& soln,
SparseMatrix<Number>& jacobian,
NonlinearImplicitSystem& sys)
{
EquationSystems &es = sys.get_equation_systems();
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
TransientNonlinearImplicitSystem& system =
es.get_system<TransientNonlinearImplicitSystem>("sys-one");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
std::vector<unsigned int> dof_indices;
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Ke.resize (dof_indices.size(),
dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Gradient grad_u;
Number du_dot_du = 0;
for (unsigned int i=0; i<phi.size(); i++)
{
grad_u += dphi[i][qp]*soln(dof_indices[i]);
du_dot_du += phi[i][qp]*(1./dt);
}
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
Ke(i,j) += JxW[qp]*(
(phi[i][qp]*phi[j][qp]*du_dot_du) +
(dphi[i][qp]*dphi[j][qp]) +
(vel * dphi[j][qp] * phi[i][qp])
);
}
dof_map.constrain_element_matrix (Ke, dof_indices);
jacobian.add_matrix (Ke, dof_indices);
}
}
void compute_residual (const NumericVector<Number>& soln,
NumericVector<Number>& residual,
NonlinearImplicitSystem& sys)
{
EquationSystems &es = sys.get_equation_systems();
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
libmesh_assert (dim == 2);
TransientNonlinearImplicitSystem& system =
es.get_system<TransientNonlinearImplicitSystem>("sys-one");
NumericVector<Number> & old_soln = *system.old_local_solution;
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseVector<Number> Re;
std::vector<unsigned int> dof_indices;
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)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Re.resize (dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Number u = 0;
Gradient grad_u;
Number u_dot = 0;
for (unsigned int j=0; j<phi.size(); j++)
{
u += phi[j][qp]*soln(dof_indices[j]);
grad_u += dphi[j][qp]*soln(dof_indices[j]);
u_dot += phi[j][qp]*(soln(dof_indices[j]) -
old_soln(dof_indices[j])) / dt;
}
for (unsigned int i=0; i<phi.size(); i++)
Re(i) += JxW[qp]*(
(phi[i][qp]*u_dot) +
(dphi[i][qp]*grad_u) +
(vel * grad_u * phi[i][qp])
);
}
dof_map.constrain_element_vector (Re, dof_indices);
residual.add_vector (Re, dof_indices);
}
residual.close();
/// Boundary node list (node ids and corresponding side-set ids, arrays
always have the same length)
std::vector<unsigned int> nodes;
std::vector<short int> ids;
mesh.boundary_info->build_node_list(nodes, ids);
for (unsigned int i = 0; i < nodes.size(); i++)
{
const Node & node = mesh.node(nodes[i]);
short int boundary_id = ids[i];
if(node.processor_id() == libMesh::processor_id())
{
unsigned int dof = node.dof_number(system.number(), 0, 0);
if (boundary_id == 1)
{
residual.set(dof, soln(dof) - 0.);
}
else if (boundary_id == 2)
{
residual.set(dof, soln(dof) - 1.);
}
}
}
residual.close();
}
// Begin the main program. Note that the first
// statement in the program throws an error if
// you are in complex number mode, since this
// example is only intended to work with real
// numbers.
int main (int argc, char** argv)
{
// Initialize libMesh.
LibMeshInit init (argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
// Our Trilinos interface does not yet support adaptive transient
// problems
libmesh_example_assert(libMesh::default_solver_package() != TRILINOS_SOLVERS,
"--enable-petsc");
// This value is also obtained from the commandline and it specifies the
// initial value for the t_step looping variable. We must
// distinguish between the two cases here, whether we read in the
// solution or we started from scratch, so that we do not overwrite the
// gmv output files.
unsigned int init_timestep = 0;
// This value is also obtained from the command line, and specifies
// the number of time steps to take.
unsigned int n_timesteps = 2;
// Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
// Create a new mesh.
Mesh mesh;
// Create an equation systems object.
EquationSystems equation_systems (mesh);
MeshRefinement mesh_refinement (mesh);
// Read the mesh from file.
mesh.read ("square.e");
mesh.prepare_for_use(false);
mesh.boundary_info->build_node_list_from_side_list();
// Print information about the mesh to the screen.
mesh.print_info();
// Declare the system and its variables.
// Begin by creating a transient system
TransientNonlinearImplicitSystem & system =
equation_systems.add_system<TransientNonlinearImplicitSystem>("sys-one");
// Add another empty system
TransientExplicitSystem &sys2 =
equation_systems.add_system<TransientExplicitSystem>("sys-two");
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 "Convection-Diffusion". "u"
// will be approximated using first-order approximation.
system.add_variable ("u", FIRST);
// Give the system a pointer to the matrix assembly
// and initialization functions.
system.nonlinear_solver->residual = compute_residual;
system.nonlinear_solver->jacobian = compute_jacobian;
// Initialize the data structures for the equation system.
equation_systems.init ();
// Prints information about the system to the screen.
equation_systems.print_info();
equation_systems.parameters.set<unsigned int>
("linear solver maximum iterations") = 250;
equation_systems.parameters.set<Real>
("linear solver tolerance") = 10e-10;
// The Convection-Diffusion system requires that we specify
// the flow velocity. We will specify it as a RealVectorValue
// data type and then use the Parameters object to pass it to
// the assemble function.
equation_systems.parameters.set<RealVectorValue>("velocity") =
RealVectorValue (0.8, 0.8);
// The Convection-Diffusion system also requires a specified
// diffusivity. We use an isotropic (hence Real) value.
equation_systems.parameters.set<Real>("diffusivity") = 0.01;
// Solve the system "Convection-Diffusion". This will be done by
// looping over the specified time interval and calling the
// \p solve() member at each time step. This will assemble the
// system and call the linear solver.
Real time = init_timestep*dt;
// We do 25 timesteps both before and after writing out the
// intermediate solution
for(unsigned int t_step=init_timestep;
t_step<(init_timestep+n_timesteps);
t_step++)
{
// Increment the time counter, set the time and the
// time step size as parameters in the EquationSystem.
time += dt;
equation_systems.parameters.set<Real> ("time") = time;
equation_systems.parameters.set<Real> ("dt") = dt;
// A pretty update message
std::cout << " Solving time step ";
// As already seen in example 9, use a work-around
// for missing stream functionality (of older compilers).
// Add a set of scope braces to enforce data locality.
{
OStringStream out;
OSSInt(out,2,t_step);
out << ", time=";
OSSRealzeroleft(out,6,3,time);
out << "...";
std::cout << out.str() << std::endl;
}
// At this point we need to update the old
// solution vector. The old solution vector
// will be the current solution vector from the
// previous time step. We will do this by extracting the
// system from the \p EquationSystems object and using
// vector assignment. Since only \p TransientLinearImplicitSystems
// (and systems derived from them) contain old solutions
// we need to specify the system type when we ask for it.
TransientNonlinearImplicitSystem & system =
equation_systems.get_system<TransientNonlinearImplicitSystem>("sys-one");
*system.old_local_solution = *system.current_local_solution;
TransientExplicitSystem & sys2 =
equation_systems.get_system<TransientExplicitSystem>("sys-two");
*sys2.old_local_solution = *sys2.current_local_solution;
// The number of refinement steps per time step.
const unsigned int max_r_steps = 1;
// A refinement loop.
{
// Assemble & solve the linear system
system.solve();
system.update();
// Print out the H1 norm, for verification purposes:
Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "H1 norm = " << H1norm << std::endl;
char name[256] = { 0 };
sprintf(name, "step-%04d.e", t_step);
ExodusII_IO(mesh).write_equation_systems (name, equation_systems);
// refine the mesh
{
std::cout << " Refining the mesh..." << std::endl;
// The \p ErrorVector is a particular \p StatisticsVector
// for computing error information on a finite element mesh.
ErrorVector error;
// The \p ErrorEstimator class interrogates a finite element
// solution and assigns to each element a positive error value.
// This value is used for deciding which elements to refine
// and which to coarsen.
//ErrorEstimator* error_estimator = new KellyErrorEstimator;
KellyErrorEstimator error_estimator;
// Compute the error for each active element using the provided
// \p flux_jump indicator. Note in general you will need to
// provide an error estimator specifically designed for your
// application.
error_estimator.estimate_error (system, error);
// This takes the error in \p error and decides which elements
// will be coarsened or refined. Any element within 20% of the
// maximum error on any element will be refined, and any
// element within 7% of the minimum error on any element might
// be coarsened. Note that the elements flagged for refinement
// will be refined, but those flagged for coarsening _might_ be
// coarsened.
mesh_refinement.refine_fraction() = 0.2;
mesh_refinement.coarsen_fraction() = 0.3;
mesh_refinement.max_h_level() = 4;
mesh_refinement.flag_elements_by_error_fraction (error);
// This call actually refines and coarsens the flagged
// elements.
mesh_refinement.refine_and_coarsen_elements();
// This call reinitializes the \p EquationSystems object for
// the newly refined mesh. One of the steps in the
// reinitialization is projecting the \p solution,
// \p old_solution, etc... vectors from the old mesh to
// the current one.
equation_systems.reinit ();
mesh.boundary_info->build_node_list_from_side_list();
mesh.print_info();
}
}
}
{
// Print out the H1 norm of the saved solution, for verification purposes:
TransientNonlinearImplicitSystem & system =
equation_systems.get_system<TransientNonlinearImplicitSystem>("sys-one");
Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Final H1 norm = " << H1norm << std::endl << std::endl;
ExodusII_IO(mesh).write_equation_systems ("out.e", equation_systems);
}
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
------------------------------------------------------------------------------ EditLive Enterprise is the world's most technically advanced content authoring tool. Experience the power of Track Changes, Inline Image Editing and ensure content is compliant with Accessibility Checking. http://p.sf.net/sfu/ephox-dev2dev
_______________________________________________ Libmesh-devel mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-devel
