On Thu, 5 Mar 2009, Tim Kroeger wrote:
Well, actually the crash was due to the missing a.close() between a=b
and a.scale(). After adding that, but before my latest patch, it didn't
crash any more, but it produced totally wrong results. That led me to the
idea that PetscVector::scale() (and hence other methods) should use the local
form.
This is surprising. What were you calling scale() on? As of now the
only ghosted vectors are *supposed* to be current_local_solution and kin.
However, I how have another problem: After the first ~20 time steps of my
application, it starts to refine the grid and -- crashes. I have no idea
why, but I will try to track that down.
Thanks. The ghosted vectors seem to do fine with our AMR/C examples
and with my apps, and without a test case demonstrating the bug I
don't know where I'd begin to fix it.
However, (a) the results do not *exactly* coincide with the previous ones;
that it, I get e.g. a residual norm of 1.16562e-08 instead of 1.16561e-08
after 608 linear iterations (although using the same number of
processors),
Hmm... maybe we still have some inaccuracy bugs in there? There
aren't any other confounding changes? If you configure using
--disable-ghosted you get back the old results?
I haven't tested this yet. I will try to find the cause of the crash first.
Right, first things first.
(Re-configuring takes a lot of time, you know.)
Actually, I do all these big tests on a cluster with distcc set up.
"make -j 50" is *fantastic*. ;-)
and (b) I had to add a.close() between a=b and a.scale(), where a and b
are *System::solution of two different systems.
The systems are identical? Same mesh, same variable types added in
the same order?
The mesh is the same, and also the variables. One is an ExplicitSystem, the
other a LinearImplicitSystem, but that shouldn't be important, should it?
No, they should have the same vector index ordering.
I don't know what to think yet... but I actually can't replicate
problem (b) myself with the attached vectest.C -
I can't either. Very strange.
Odd. How about with this vectest2.C? I can't even get problem (b) to
trigger with ghosted vectors!
On the other hand, here's something even more odd: I just realized
that I'm not properly matching the ghosted vector init() function
signature. The compiler must be deciding to convert GHOSTED to true
and then using the default AUTOMATIC for the fifth argument? That
wouldn't explain any crashes, but if we're fast-initializing vectors
we should be clearing, that might explain your slightly-off residuals.
I'll fix it in SVN.
---
Roy
// C++ include files that we need
#include <iostream>
#include <cmath>
// Basic include file needed for the mesh functionality.
#include "libmesh.h"
#include "numeric_vector.h"
const unsigned int dofsperproc = 1;
int main (int argc, char** argv)
{
// Initialize libMesh.
LibMeshInit init (argc, argv);
std::vector<unsigned int> send_list;
if (libMesh::processor_id() > 0)
send_list.push_back(libMesh::processor_id() * dofsperproc - 1);
if (libMesh::processor_id() < libMesh::n_processors() - 1)
send_list.push_back(libMesh::processor_id() * dofsperproc + dofsperproc +
1);
AutoPtr<NumericVector<Number> > vec1 = NumericVector<Number>::build();
AutoPtr<NumericVector<Number> > vec2 = NumericVector<Number>::build();
// vec1->init(libMesh::n_processors() * dofsperproc, dofsperproc, false,
PARALLEL);
// vec2->init(libMesh::n_processors() * dofsperproc, dofsperproc, false,
PARALLEL);
vec1->init(libMesh::n_processors() * dofsperproc, dofsperproc, send_list,
GHOSTED);
vec2->init(libMesh::n_processors() * dofsperproc, dofsperproc, send_list,
GHOSTED);
*vec1 = *vec2;
vec1->scale(1.0);
return 0;
}
------------------------------------------------------------------------------
Open Source Business Conference (OSBC), March 24-25, 2009, San Francisco, CA
-OSBC tackles the biggest issue in open source: Open Sourcing the Enterprise
-Strategies to boost innovation and cut costs with open source participation
-Receive a $600 discount off the registration fee with the source code: SFAD
http://p.sf.net/sfu/XcvMzF8H
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
