From: Vetter R. <vet...@st...> - 2011-02-09 09:57:40
|
Dear libmesh developers VTK output crashes when run in parallel (mpirun -n 2 ...). The following assertion fails: Assertion `it!=_global_to_local_map.end()' failed. [...]/libmesh/include/numerics/petsc_vector.h, line 1073, [...] The problem is reproducible with the following minimal working example: //START MINIMAL WORKING EXAMPLE CODE #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "equation_systems.h" #include "explicit_system.h" #include "vtk_io.h" using namespace libMesh; int main(int argc, char** argv) { LibMeshInit init(argc, argv); { Mesh mesh; MeshTools::Generation::build_line(mesh, 4, 0, 1, EDGE2); EquationSystems equation_systems(mesh); ExplicitSystem& system = equation_systems.add_system<ExplicitSystem>("sys"); system.add_variable("u"); equation_systems.init(); VTKIO(mesh).write_equation_systems("out.pvtu", equation_systems); // causes assertion failure } return 0; } //END MINIMAL WORKING EXAMPLE CODE The stack trace looks as follows: 0: libMesh::print_trace(std::ostream&) 1: libMesh::PetscVector<double>::map_global_to_local_index(unsigned int) const 2: libMesh::PetscVector<double>::operator()(unsigned int) const 3: libMesh::System::current_solution(unsigned int) const 4: libMesh::VTKIO::solution_to_vtk(libMesh::EquationSystems const&, vtkUnstructuredGrid*&) 5: libMesh::VTKIO::write_equation_systems(std::string const&, libMesh::EquationSystems const&) 6: main 7: __libc_start_main The assertion does NOT fail when running in serial mode (without mpirun) or when using ExodusII_IO. Since I'm bound to VTK and parallel processing, I will need to fix this. The problem seems to be that current_solution.operator() is called for ALL dofs in line 333 of vtk_io.C, instead of only the local ones. Therefore, the global-to-local index map doesn't contain the sought indices, when processor 0 attempts to write data owned by processor 1. Solution output appears to be more sophisticated in ExodusII_IO, where MeshOutput<MT>::_build_variable_names_and_solution_vector is called, which does a EquationSystems::allgather(). What's the expert's ansatz to fixing VTK output? Should I have a go at duplicating the ExodusII output procedure, or does anybody see a simple(r) way of fixing the problem? Perhaps a localize_to_one() call on the solution vector? Would that compromise parallel efficiency? |
From: Roy S. <roy...@ic...> - 2011-02-09 20:59:45
|
On Wed, 9 Feb 2011, Vetter Roman wrote: > VTK output crashes when run in parallel (mpirun -n 2 ...). The following assertion fails: > > Assertion `it!=_global_to_local_map.end()' failed. > [...]/libmesh/include/numerics/petsc_vector.h, line 1073, [...] Yup - the VTK support was a contribution from a non-main developer, one who only needed and only wrote support for output of serial non-adaptive data. We figured that any capability was better than none, but I suppose at the very least we should have tossed a few lines like "if (n_processors>1) libmesh_not_implemented()" to give people more comprehensible warnings. > The assertion does NOT fail when running in serial mode (without > mpirun) or when using ExodusII_IO. Since I'm bound to VTK and > parallel processing, I will need to fix this. We'd appreciate that, thanks! > The problem seems to be that current_solution.operator() is called > for ALL dofs in line 333 of vtk_io.C, instead of only the local > ones. Therefore, the global-to-local index map doesn't contain the > sought indices, when processor 0 attempts to write data owned by > processor 1. Exactly. > Solution output appears to be more sophisticated in ExodusII_IO, > where MeshOutput<MT>::_build_variable_names_and_solution_vector is > called, which does a EquationSystems::allgather(). EquationSystems::allgather() is actually a second step up in sophistication from where VTK output is - it gathers geometry data and degree of freedom numbering in the case of a distributed mesh. The gathering (and sadly, the decimation of higher-order solutions) of vector data happens in EquationSystems::build_vector(); > What's the expert's ansatz to fixing VTK output? Should I have a go > at duplicating the ExodusII output procedure, I'm not an expert on the I/O classes, but that's probably the best way. ExodusII is at least decently maintained by people who are experts. ;-) > or does anybody see a simple(r) way of fixing the problem? Perhaps a > localize_to_one() call on the solution vector? Would that compromise > parallel efficiency? Right now *every* single-file I/O we do compromises parallel efficiency; if that's a serious bottleneck for you then you'll either need to write a support for the parallel multi-file VTK format or switch to Nemesis. --- Roy |
From: Roman V. <vet...@st...> - 2011-02-21 13:19:52
|
> EquationSystems::allgather() is actually a second step up in > sophistication from where VTK output is - it gathers geometry data and > degree of freedom numbering in the case of a distributed mesh. The > gathering (and sadly, the decimation of higher-order solutions) of > vector data happens in EquationSystems::build_vector(); Oh, I see. Well, it would certainly be preferable to have also the higher order components written. From that point of view, I might be better for me not to use the existing procedure from ExodusII, but address the VTK output independently. I'll look into this soon. In the meantime, below is a patch that fixes two rather serious bugs in VTK output (unrelated to parallelism), and removes annoying debugging output leftovers from the original contributor. The first bug fixed is this one: http://sourceforge.net/tracker/?func=detail&aid=3140652&group_id=71130&atid=530254 The second one was that the values ranges were printed only in single precision instead of double, which causes trouble in third-party VTK readers such as in Paraview. Roman Index: src/mesh/vtk_io.C =================================================================== --- src/mesh/vtk_io.C (revision 4198) +++ src/mesh/vtk_io.C (working copy) @@ -317,11 +317,11 @@ std::string name = sys.variable_name(j); - vtkFloatArray *data = vtkFloatArray::New(); + vtkDoubleArray *data = vtkDoubleArray::New(); data->SetName(name.c_str()); - data->SetNumberOfValues(sys.solution->size()); + data->SetNumberOfValues(n_nodes); for(unsigned int k=0;k<n_nodes;++k){ @@ -601,16 +601,12 @@ * we only use Unstructured grids */ _vtk_grid = vtkUnstructuredGrid::New(); - vtkXMLPUnstructuredGridWriter* writer= vtkXMLPUnstructuredGridWriter::New(); - libMesh::out<<"get points "<<std::endl; + vtkXMLPUnstructuredGridWriter* writer= vtkXMLPUnstructuredGridWriter::New(); vtkPoints* pnts = nodes_to_vtk((const MeshBase&)es.get_mesh()); _vtk_grid->SetPoints(pnts); - int * types = new int[es.get_mesh().n_active_elem()]; - libMesh::out<<"get cells"<<std::endl; + int * types = new int[es.get_mesh().n_active_elem()]; vtkCellArray* cells = cells_to_vtk((const MeshBase&)es.get_mesh(), types); - - libMesh::out<<"set cells"<<std::endl; _vtk_grid->SetCells(types,cells); // I'd like to write out meshdata, but this requires some coding, in @@ -618,8 +614,7 @@ // const MeshData& md = es.get_mesh_data(); // if(es.has_mesh_data()) // meshdata_to_vtk(md,_vtk_grid); - // libmesh_assert (soln.size() ==mesh.n_nodes()*names.size()); - libMesh::out<<"write solution"<<std::endl; + // libmesh_assert (soln.size() ==mesh.n_nodes()*names.size()); solution_to_vtk(es,_vtk_grid); //#ifdef DEBUG |
From: Roy S. <roy...@ic...> - 2011-02-21 17:01:48
|
On Mon, 21 Feb 2011, Roman Vetter wrote: >> EquationSystems::allgather() is actually a second step up in >> sophistication from where VTK output is - it gathers geometry data and >> degree of freedom numbering in the case of a distributed mesh. The >> gathering (and sadly, the decimation of higher-order solutions) of >> vector data happens in EquationSystems::build_vector(); > > Oh, I see. Well, it would certainly be preferable to have also the > higher order components written. From that point of view, I might be > better for me not to use the existing procedure from ExodusII, but > address the VTK output independently. I'll look into this soon. Does VTK have native support for higher-order data? If so, adding that to our VTKIO output would be fantastic. If not, then it ought to be possible to implement a more general solution. "Subdivide an order-p element into p^d order-1 elements" is a pretty file-format-independent task, even if we limited the initial testing to VTKIO. > In the meantime, below is a patch that fixes two rather serious bugs > in VTK output (unrelated to parallelism), and removes annoying > debugging output leftovers from the original contributor. > > The first bug fixed is this one: > http://sourceforge.net/tracker/?func=detail&aid=3140652&group_id=71130&atid=530254 > The second one was that the values ranges were printed only in single > precision instead of double, which causes trouble in third-party VTK > readers such as in Paraview. Patch didn't apply cleanly for me (submit as an attachment next time? my email system must have mangled the whitespace or something) but it looks simple enough to apply by hand. I'll separate the bug fixes and the de-verbosifying into two separater commits and put them both in SVN now. --- Roy |
From: Roy S. <roy...@ic...> - 2011-02-21 17:04:40
|
On Mon, 21 Feb 2011, Roy Stogner wrote: > Patch didn't apply cleanly for me (submit as an attachment next time? > my email system must have mangled the whitespace or something) but it > looks simple enough to apply by hand. And I'm glad I looked at it more closely - shouldn't we also be using vtkDoubleArray in system_vectors_to_vtk? It looks as if your patch bumps up the precision on system.solution but not on the other vectors being output. --- Roy |
From: Roman V. <vet...@st...> - 2011-02-22 14:33:21
Attachments:
patch.diff
|
>> Patch didn't apply cleanly for me (submit as an attachment next time? >> my email system must have mangled the whitespace or something) but it >> looks simple enough to apply by hand. Sorry about that. I'll attach patches as files in the future. > > And I'm glad I looked at it more closely - shouldn't we also be using > vtkDoubleArray in system_vectors_to_vtk? It looks as if your patch > bumps up the precision on system.solution but not on the other vectors > being output. Yes, that would seem like a good idea to me. However, system_vectors_to_vtk is not actively being used (the only call is commented and even bracketed by #ifdef DEBUG), and frankly, it is a huge mess. I suggest commenting that method out, with a remark that it doesn't support AMR or higher order variables (it really doesn't -- it relies on n_dof() == n_nodes(), which is fatal). Besides, I don't think that VTKIO::write_equation_systems should write system vectors anyway. Neither of ExodusII_IO or GMVIO do. > Does VTK have native support for higher-order data? If so, adding > that to our VTKIO output would be fantastic. Yes, VTK does, but it's not implemented in VTKIO. I attached a patch that does it. :-) Note: The patch only modifies solution_to_vtk, while leaving system_vectors_to_vtk untouched. See above for the reason why... > If not, then it ought to be possible to implement a more general > solution. "Subdivide an order-p element into p^d order-1 elements" is > a pretty file-format-independent task, even if we limited the initial > testing to VTKIO. Agreed. It's just that VTK output is implemented fundamentally different than others. But it should be relatively easy to modify EquationSystems::build_solution_vector and EquationSystems::build_variable_names, which seem to be only used by ExodusII_IO and GMVIO currently, in the same way, i.e. such that build_variable_names also returns all components with a suffix, and build_solution_vector returns also the higher order components. It seems to me that this would be sufficient already. Roman |
From: Roy S. <roy...@ic...> - 2011-02-22 19:56:47
|
>> Does VTK have native support for higher-order data? If so, adding >> that to our VTKIO output would be fantastic. > > Yes, VTK does, but it's not implemented in VTKIO. I attached a patch > that does it. :-) Does it, but how well? It looks like the higher order coefficients are getting added to the file, fine, but I don't see how they're supposed to end up being accurately visualized. How is a VTK reader going to know whether a given higher order DoF index is the coefficient of a quadratic interpolatory function in a Lagrange basis, a quartic bubble function in a hierarchic basis, or a normal flux function in a cubic macroelement Clough-Tocher basis? --- Roy |
From: Roman V. <vet...@st...> - 2011-02-23 09:06:07
|
2011/2/22, Roy Stogner <roy...@ic...>: > >>> Does VTK have native support for higher-order data? If so, adding >>> that to our VTKIO output would be fantastic. >> >> Yes, VTK does, but it's not implemented in VTKIO. I attached a patch >> that does it. :-) > > Does it, but how well? It looks like the higher order coefficients > are getting added to the file, fine, but I don't see how they're > supposed to end up being accurately visualized. How is a VTK reader > going to know whether a given higher order DoF index is the > coefficient of a quadratic interpolatory function in a Lagrange basis, > a quartic bubble function in a hierarchic basis, or a normal flux > function in a cubic macroelement Clough-Tocher basis? Oh, that's what you mean by "native support". Well, in that case VTK doesn't. It only allows multiple components per node/element ("tuples"), but does not specify what these components stand for. A VTK reader will therefore not be able to visualize the shape functions/solution properly, without the user telling the software how to interpret the data. Specifying interpolation types is just not part of the VTK file format. I suppose it's a matter of opinion/objective. Do the libMesh developers/users want their solution output either a) complete, including higher order terms, in a way that complies with the VTK file format, but not necessarily perfectly visualizable in, say, Paraview without further ado in all cases or b) incomplete, i.e. only first order terms, not using the full capability of the VTK file format, but with a higher chance that VTK readers will be able to deal with at least this first order component properly ? Case a) is implemented with my patch, case b) is the current way of implementation. There is a third possibility, which is a compromise: In the case of higher order variables, one could write each component separately as scalar variables, perhaps with a suffix in the name denoting the component. This way, all information is written, while Paraview will still be able to recognize the first order solution as scalars and display it as it does currently. Roman |
From: Roy S. <roy...@ic...> - 2011-02-23 14:51:46
|
On Wed, 23 Feb 2011, Roman Vetter wrote: > I suppose it's a matter of opinion/objective. Do the libMesh > developers/users want their solution output either > > a) complete, including higher order terms, in a way that complies with > the VTK file format, but not necessarily perfectly visualizable in, > say, Paraview without further ado in all cases At least for me this isn't useful until/unless there's *some* 3rd party support for reading/understanding the higher-order solution terms, and it's actually counterproductive if it makes it less likely for a reader to render the low-order decimated solution properly. I'm not sure how it would be useful for others. As a first step toward eventually having another option for restart files, I assume? Right now (IIRC, Derek please correct me if I'm wrong) we only support higher order restarts via our own xda/xdr based format, with a few more low order solution options. I'd like to hear more opinions, if anyone wants to chime in. I think Ben's busy with NASA stuff and John's busy moving across the country, though... > or > > b) incomplete, i.e. only first order terms, not using the full > capability of the VTK file format, but with a higher chance that VTK > readers will be able to deal with at least this first order component > properly > > ? > > Case a) is implemented with my patch, case b) is the current way of > implementation. I'm definitely not happy with the current implementation either, but the big flaws in it are its failures in parallel and on adapted grids; I won't think of p-decimation as a flaw until I see a viz tool that properly handles high-p output. > There is a third possibility, which is a compromise: > In the case of higher order variables, one could write each component > separately as scalar variables, perhaps with a suffix in the name > denoting the component. This way, all information is written, while > Paraview will still be able to recognize the first order solution as > scalars and display it as it does currently. See my response to option a), but without the "counterproductive" worry. ;-) This actually sounds like a decent improvement, if it works for you. I'd worry about it being incompatible with future 3rd-party support for higher-order bases, but that's a problem with any idea we might have for saving non-decimated solutions. --- Roy |
From: Derek G. <fri...@gm...> - 2011-02-28 17:58:31
Attachments:
Hermites_0lvls_early.png
Hermites_2lvls_early.png
|
Sorry I haven't been more active in this conversation... I was on vacation all of last week. I've sprinkled my thoughts below.... On Feb 23, 2011, at 7:51 AM, Roy Stogner wrote: > > On Wed, 23 Feb 2011, Roman Vetter wrote: > >> I suppose it's a matter of opinion/objective. Do the libMesh >> developers/users want their solution output either >> >> a) complete, including higher order terms, in a way that complies with >> the VTK file format, but not necessarily perfectly visualizable in, >> say, Paraview without further ado in all cases > > At least for me this isn't useful until/unless there's *some* 3rd > party support for reading/understanding the higher-order solution > terms, and it's actually counterproductive if it makes it less likely > for a reader to render the low-order decimated solution properly. > > I'm not sure how it would be useful for others. As a first step > toward eventually having another option for restart files, I assume? > Right now (IIRC, Derek please correct me if I'm wrong) we only support > higher order restarts via our own xda/xdr based format, with a few > more low order solution options. > > I'd like to hear more opinions, if anyone wants to chime in. I think > Ben's busy with NASA stuff and John's busy moving across the country, > though... Agree with Roy. If there is any chance that a user will use Hermites (or any higher order shape functions) and the _default_ result file won't load in Paraview.... that is a FULL STOP. If you want to have an option on the VTK_IO writer to write the extra dofs out.... that's fine with me. But the default should always be compatible with 3rd party viz tools. This is actually a timely discussion for us because we are having trouble visualizing our higher order solutions. We have actually added capability to our code to project higher order solutions onto arbitrarily refined meshes for output.... so that we can actually get a better idea of what the true solution looks like. If _anyone_ knows of a 3rd party viz tool that will visualize higher order solutions.... I would REALLY love to hear about it. The alternative to doing what I just said (projecting to refined meshes) is doing what Roy mentioned a few emails ago: refining _when_ you output. That can be achieved but seemed like a lot more work than just doing a mesh to mesh projection and then using our normal IO code to write the solution out. I wouldn't mind some other opinions here though. As for restart... Roy you are right. We can only _restart_ from xda/xdr in the case of higher order shape functions. We can read "nodal values" out of both Exodus and GMV I believe. But you will obviously not have any of the higher order coefficients. To sum up: 1. Default behavior must work with 3rd party tools 2. If there is a tool that can visualize higher order data.... I want to know about it! 3. There are other alternatives for dealing with higher order data like projecting to finer grids. BTW.... I'm attaching two images of a Cahn-Hilliard solution run with 3rd order hermites.... the first is on the original grid.... the second is the same solution... but projected to a finer grid before being written out. There is no doubt that by visualizing the decimated solution on the original grid we're missing a TON of info! (I believe this visualization was done with Paraview.... but I don't know for sure) Derek |
From: Roman V. <vet...@st...> - 2011-02-24 11:25:47
Attachments:
parallel_vtk.diff
|
>> There is a third possibility, which is a compromise: >> In the case of higher order variables, one could write each component >> separately as scalar variables, perhaps with a suffix in the name >> denoting the component. This way, all information is written, while >> Paraview will still be able to recognize the first order solution as >> scalars and display it as it does currently. > > [...] This actually sounds like a decent improvement, if it > works for you. [...] I have implemented this, and it works fine with Paraview. Would you like me to supply the patch so you (guys) can see what it actually looks like? Getting back to the original motivation for this discussion: Parallelism in VTK output. I have attached a patch which fixes the issue for me. There are certainly many ways of doing it, mine attempts to somewhat minimize the amount of modifications of the current code. Essentially, the solution vectors are localized to one. Roman |
From: Roy S. <roy...@ic...> - 2011-02-28 23:21:24
|
On Thu, 24 Feb 2011, Roman Vetter wrote: >>> There is a third possibility, which is a compromise: >>> In the case of higher order variables, one could write each component >>> separately as scalar variables, perhaps with a suffix in the name >>> denoting the component. This way, all information is written, while >>> Paraview will still be able to recognize the first order solution as >>> scalars and display it as it does currently. >> >> [...] This actually sounds like a decent improvement, if it >> works for you. [...] > > I have implemented this, and it works fine with Paraview. Would you > like me to supply the patch so you (guys) can see what it actually > looks like? Not speaking for Derek, but for me the answer is definitely "yes"! > Getting back to the original motivation for this discussion: > Parallelism in VTK output. I have attached a patch which fixes the > issue for me. There are certainly many ways of doing it, mine attempts > to somewhat minimize the amount of modifications of the current code. > Essentially, the solution vectors are localized to one. This definitely looks like an improvement over the current code. Seems to work on my initial tests, too; I'll commit it to svn. Works on the initial tests where it's expected to work, anyway. Does the newest VTK support triquadratic hexes yet? vtk_io.C seems to think that VTK_BIQUADRATIC_QUAD is relatively new, and it would be nice if we didn't have to die with an error in the 3D case. --- Roy |
From: Roy S. <roy...@ic...> - 2011-03-01 20:30:28
|
On Mon, 28 Feb 2011, Roy Stogner wrote: > Works on the initial tests where it's expected to work, anyway. Does > the newest VTK support triquadratic hexes yet? vtk_io.C seems to > think that VTK_BIQUADRATIC_QUAD is relatively new, and it would be > nice if we didn't have to die with an error in the 3D case. Looks like the answer is yes - if I find time I'll put together the proper Hex27::connectivity and add support. There's another test where we die with an error: mixed-order elements. E.g. ex11 with VTKIO dies screaming as soon as VTKIO::solution_to_vtk hits a side node on a QUAD9, when it expects to see two velocity components plus one pressure component but it only sees one pressure component. I don't know enough about the VTK formats to say what the right way to handle this is. --- Roy |
From: Roman V. <vet...@st...> - 2011-03-03 09:14:46
Attachments:
vtk_high_order.diff
|
> There's another test where we die with an error: mixed-order elements. > E.g. ex11 with VTKIO dies screaming as soon as VTKIO::solution_to_vtk > hits a side node on a QUAD9, when it expects to see two velocity > components plus one pressure component but it only sees one pressure > component. I don't know enough about the VTK formats to say what the > right way to handle this is. Judging from the latest VTK file format specification (4.2) apparently existing on the web, mixed-number-of-component entities seem unsupported. A field can be scalar or vectorial, but has a fixed number of components for all nodes/elements either way. Anyway, here's my patch that adds higher order components as additional scalars with an integer suffix in the name to VTK output files, leaving the first order component untouched. Roman |
From: Roy S. <roy...@ic...> - 2011-03-04 16:45:18
|
On Thu, 3 Mar 2011, Roman Vetter wrote: > Judging from the latest VTK file format specification (4.2) apparently > existing on the web, mixed-number-of-component entities seem > unsupported. A field can be scalar or vectorial, but has a fixed > number of components for all nodes/elements either way. That's a shame, but understandable. > Anyway, here's my patch that adds higher order components as > additional scalars with an integer suffix in the name to VTK output > files, leaving the first order component untouched. Looking over this... would it make sense to interpolate instead when we find nodes with missing Lagrange DoFs, to get an isoparametric result? VTK supports second order elements, and first order data is a subset of second order data, so we could still naturally represent it on those elements. --- Roy |
From: Roman V. <vet...@st...> - 2011-03-09 13:37:45
|
> Looking over this... would it make sense to interpolate instead when > we find nodes with missing Lagrange DoFs, to get an isoparametric > result? VTK supports second order elements, and first order data is a > subset of second order data, so we could still naturally represent it > on those elements. I'm probably not the right person to tell, as I haven't been in the situation of using such elements. From a pure VTK point of view, I'd say why not, but I don't know how people or third party visualizers would expect to find their data in this case. Roman |
From: Roy S. <roy...@ic...> - 2011-03-09 21:02:05
|
On Wed, 9 Mar 2011, Roman Vetter wrote: >> Looking over this... would it make sense to interpolate instead when >> we find nodes with missing Lagrange DoFs, to get an isoparametric >> result? VTK supports second order elements, and first order data is a >> subset of second order data, so we could still naturally represent it >> on those elements. > > I'm probably not the right person to tell, as I haven't been in the > situation of using such elements. From a pure VTK point of view, I'd > say why not, but I don't know how people or third party visualizers > would expect to find their data in this case. It should be the same as for a truely isoparametric result, no? If the visualizer can handle a biquadratic element with three biquadratic variables, then by definition it can handle a biquadratic element with three biquadratic variables even when one of those variables suspiciously only takes on values that could also be expressed in a bilinear basis. --- Roy |