The attached code demonstrates what I mean, and below is a summary of some representative results for a case with 150 variables (note that I got the "average bandwidth" data from equations_systems.print_info()):
MONOMIAL with diagonal dof_coupling: average bandwidth = 811 MONOMIAL with all variables coupled: average bandwidth = 960 LAGRANGE with dof_coupling: average bandwidth = 22.3824LAGRANGE with all variables coupled: average bandwidth = [Don't know, my laptop ran out of memory trying to get this value... but it's certainly greater than 1000]
So LAGRANGE behaves the way I would expect, in that the dof_coupling leads to a significant reduction in matrix bandwidth (and hence a significant reduction in memory footprint). Any ideas why we don't get the same behavior with MONOMIAL? Also, suggestions on how to fix this for MONOMIAL would be appreciated!
Thanks, David
######################################################################
#
# Template libMesh application Makefile
LIBMESH_DIR ?= /home/dknez/software/libmesh-real
# include the library options determined by configure
include $(LIBMESH_DIR)/Make.common
target := ./sparsity_test-$(METHOD)
###############################################################################
# File management. This is where the source, header, and object files are
# defined
#
# source files
srcfiles := $(wildcard *.C)
#
# object files
objects := $(patsubst %.C, %.$(obj-suffix), $(srcfiles))
###############################################################################
.PHONY: dust clean distclean
###############################################################################
# Target:
#
all:: $(notdir $(target))
# Production rules: how to make the target - depends on library configuration
$(notdir $(target)): $(objects)
@echo "Linking "$@"..."
@$(libmesh_LIBTOOL) --tag=CXX $(LIBTOOLFLAGS) --mode=link \
$(libmesh_CXX) $(libmesh_CXXFLAGS) $(objects) -o $@ $(libmesh_LIBS)
$(libmesh_LDFLAGS) $(EXTERNAL_FLAGS)
# Useful rules.
dust:
@echo "Deleting old output and runtime files"
@rm -f out*.m job_output.txt output.txt* *.gmv.* *.plt.* out*.xdr*
out*.xda* PI* complete
clean: dust
@rm -f $(objects) *.$(obj-suffix) *.lo
clobber: clean
@rm -f $(target)
distclean: clean
@rm -rf *.o .libs .depend
echo:
@echo srcfiles = $(srcfiles)
@echo objects = $(objects)
@echo target = $(target)
run: complete
complete: $(wildcard *.in)
# @$(MAKE) dust
@$(MAKE) -C $(dir $(target)) $(notdir $(target))
@echo "***************************************************************"
@echo "* Running App " $(notdir $(target))
@echo "***************************************************************"
@echo " "
${LIBMESH_RUN} $(target) ${LIBMESH_OPTIONS} 2>&1 | tee output.txt
@bzip2 -f output.txt
@echo " "
@echo "***************************************************************"
@echo "* Done Running App " $(notdir $(target))
@echo "***************************************************************"
gmv:
@$(MAKE) -C $(LIBMESH_DIR)/roy/meshplot/ meshplot-$(METHOD)
@for file in out.mesh.*; do ${LIBMESH_RUN}
$(LIBMESH_DIR)/roy/meshplot/meshplot-$(METHOD) $$file
out.soln.$${file##out.mesh.} out.gmv.$${file:9:4}; done
# include the dependency list
-include .depend
#
# Dependencies
#
.depend: $(srcfiles) $(LIBMESH_DIR)/include/libmesh/*.h
@$(perl) $(LIBMESH_DIR)/contrib/bin/make_dependencies.pl -I. $(foreach
i, $(LIBMESH_DIR)/include $(wildcard $(LIBMESH_DIR)/include/*), -I$(i))
"-S\$$(obj-suffix)" $(srcfiles) > .depend
###############################################################################
run.sh
Description: run.sh
// Check sparsity pattern with MONOMIAL vs LAGRANGE with a diagonal coupling_matrix
// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>
#include <set>
// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/gnuplot_io.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/coupling_matrix.h"
// Define the Finite Element object.
#include "libmesh/fe.h"
// Define Gauss quadrature rules.
#include "libmesh/quadrature_gauss.h"
// Define the DofMap, which handles degree of freedom
// indexing.
#include "libmesh/dof_map.h"
// Define useful datatypes for finite element
// matrix and vector components.
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
// Define the PerfLog, a performance logging utility.
// It is useful for timing events in a code and giving
// you an idea where bottlenecks lie.
#include "libmesh/perf_log.h"
// The definition of a geometric element
#include "libmesh/elem.h"
// To impose Dirichlet boundary conditions
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/analytic_function.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/getpot.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
void assemble_matrix(EquationSystems& es,
const std::string& system_name);
// Begin the main program.
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
// Create a GetPot object to parse the command line
GetPot command_line (argc, argv);
int dim = 3;
// Create a mesh with user-defined dimension.
// Read number of elements from command line
int n_elem = 15;
if ( command_line.search(1, "-n") )
n_elem = command_line.next(n_elem);
// Read FE order from command line
std::string order = "FIRST";
if ( command_line.search(2, "-Order", "-o") )
order = command_line.next(order);
// Read FE Family from command line
std::string family = "LAGRANGE";
if ( command_line.search(2, "-FEFamily", "-f") )
family = command_line.next(family);
Mesh mesh(init.comm());
MeshTools::Generation::build_cube (mesh,
n_elem, n_elem, n_elem,
-1., 1.,
-1., 1.,
-1., 1.,
HEX8);
mesh.print_info();
EquationSystems equation_systems (mesh);
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("TestSystem");
int n_vars = 15;
if ( command_line.search(1, "-n_vars") )
n_vars = command_line.next(n_vars);
for(unsigned int var=0; var<n_vars; var++)
{
std::stringstream var_name;
var_name << "var_" << var;
system.add_variable(var_name.str(),
Utility::string_to_enum<Order> (order),
Utility::string_to_enum<FEFamily>(family));
}
// Give the system a pointer to the matrix assembly
// function.
system.attach_assemble_function (assemble_matrix);
CouplingMatrix coupling_matrix(system.n_vars());
for(unsigned int var1=0; var1<system.n_vars(); var1++)
for(unsigned int var2=0; var2<system.n_vars(); var2++)
{
unsigned char value = (var1==var2) ? 1 : 0;
coupling_matrix(var1,var2) = value;
}
int use_coupling_matrix = 1;
if ( command_line.search(1, "-use_coupling_matrix") )
use_coupling_matrix = command_line.next(use_coupling_matrix);
if(use_coupling_matrix != 0)
{
system.get_dof_map()._dof_coupling = &coupling_matrix;
}
// Initialize the data structures for the equation system.
equation_systems.init();
// Print information about the system to the screen.
equation_systems.print_info();
return 0;
}
void assemble_matrix(EquationSystems& es,
const std::string& system_name)
{
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("TestSystem");
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<Point>& q_point = fe->get_xyz();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<dof_id_type> 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());
Fe.resize (dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
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]);
}
Fe(i) += JxW[qp]*phi[i][qp];
}
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
------------------------------------------------------------------------------ Managing the Performance of Cloud-Based Applications Take advantage of what the Cloud has to offer - Avoid Common Pitfalls. Read the Whitepaper. http://pubads.g.doubleclick.net/gampad/clk?id=121051231&iu=/4140/ostg.clktrk
_______________________________________________ Libmesh-devel mailing list Libmesh-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-devel
