Dear Roy,
On Wed, 25 Mar 2009, Tim Kroeger wrote:
Unfortunately, the residuals still don't coincide. What suprises me
even more is that they even differ between identical runs, i.e. ghost
dofs are disabled both times. In other words, my application gets
results that are not reproducible, and this has got nothing to do with
ghost dofs. I will try to track this down.
I tracked it down in the sense that I can reproduce it fast, but I
have still no idea why it happens. The grid is now uniform and small.
In the matrix assembly function, I wrote the dof_indices, Ke, and Fe
to files (one file for each processor) directly before they are added
to the system matrix/rhs.
You might want to check whether you can reproduce the
non-reproducibility. To do this, please download
www.mevis.de/~tim/a.tar.gz and unpack it (small this time). Then run
the attached test.cpp on 8 CPUs, which creates a simple grid,
assembles a matrix as given by the files, writes the complete matrix
to another file (in PETSc binary format), and then solves. On solving
the system, it will crash, but that's not important in this case.
(The grid is so coarse that important geometric structures are not
"seen", so that the system becomes singular.) Run that program
several times and rename the resulting matrix files. Using "diff" you
will find that they differ. You can use the attached test2.cpp to get
an ASCII representation of the matrices.
Best Regards,
Tim
--
Dr. Tim Kroeger
tim.kroe...@mevis.fraunhofer.de Phone +49-421-218-7710
tim.kroe...@cevis.uni-bremen.de Fax +49-421-218-4236
Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
template<typename T>
static void read(std::istream& f, T& t)
{
T t0 = t;
f.read(reinterpret_cast<char*>(&t0),sizeof(t));
for(unsigned int i=0; i<sizeof(t); i++)
{
reinterpret_cast<char*>(&t)[i] =
reinterpret_cast<char*>(&t0)[sizeof(t)-1-i];
}
}
int main(int argc,char**argv)
{
if(argc!=2)
{
std::cout << "Usage: " << argv[0] << " infilename" << std::endl;
exit(1);
}
std::cout.precision(20);
std::ifstream f(argv[1]);
unsigned int m = 0;
read(f,m);//dummy
read(f,m);//this is really m
unsigned int n = 0;
read(f,n);
unsigned int numEntriesTotal = 0;
read(f,numEntriesTotal);
std::vector<std::vector<unsigned int> > nonZeros(m,std::vector<unsigned
int>());
unsigned int i = 0;
unsigned int j = 0;
unsigned int k = 0;
for(i=0; i<m; i++)
{
read(f,k);
nonZeros[i].resize(k,0);
}
for(i=0; i<m; i++)
{
for(j=0; j<nonZeros[i].size(); j++)
{
read(f,nonZeros[i][j]);
}
}
double value = 0.0;
for(i=0; i<m; i++)
{
for(j=0; j<nonZeros[i].size(); j++)
{
read(f,value);
std::cout << "[" << i << "," << j << "] = " << value << std::endl;
}
}
return 0;
}
#include <iostream>
#include <math.h>
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "mesh_refinement.h"
#include "elem.h"
#include "equation_systems.h"
#include "linear_implicit_system.h"
#include "dof_map.h"
#include "fe_base.h"
#include "quadrature_gauss.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "sparse_matrix.h"
#include "numeric_vector.h"
#include "linear_solver.h"
#include "petsc_linear_solver.h"
static void writeMatrix(std::string filename, SparseMatrix<Number>& matrix)
{
matrix.close();
Mat mat = dynamic_cast<PetscMatrix<Number>*>(&matrix)->mat();
int ierr = 0;
PetscViewer petsc_viewer;
ierr = PetscViewerCreate(libMesh::COMM_WORLD, &petsc_viewer);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = PetscViewerBinaryOpen(libMesh::COMM_WORLD, filename.c_str(),
FILE_MODE_WRITE,&petsc_viewer);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = MatView(mat, petsc_viewer);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = PetscViewerDestroy(petsc_viewer);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
static void assemble(EquationSystems& es,
const std::string& /*system_name*/)
{
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("s");
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<unsigned int> dof_indices;
std::stringstream ss;
ss << "assemble-" << libMesh::processor_id();
std::ifstream f(ss.str().c_str());
unsigned int n = 0;
f >> n;
while(!f.eof())
{
dof_indices.resize(n,0);
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
for(unsigned int i=0; i<dof_indices.size(); i++) {f >> dof_indices[i];}
for(unsigned int i=0; i<dof_indices.size(); i++)
for(unsigned int j=0; j<dof_indices.size(); j++) {f >> Ke(i,j);}
for(unsigned int i=0; i<dof_indices.size(); i++) {f >> Fe(i);}
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
f >> n;
}
writeMatrix("matrix",*system.matrix);
}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
Mesh mesh (3);
MeshTools::Generation::build_cube(mesh,
7,6,4,
-32.0,24.01,
-24.0,24.0,
-32.0,32.0,
HEX8);
EquationSystems es(mesh);
LinearImplicitSystem& s = es.add_system<LinearImplicitSystem>("s");
s.add_variable("v");
s.attach_assemble_function (assemble);
std::cout << "Init EquationSystems" << std::endl;
es.init();
std::cout << "Solve" << std::endl;
s.solve();
return 0;
}
------------------------------------------------------------------------------
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
