From: Roy Stogner <roystgnr@ic...>  20070219 19:44:57

My current work with another libMesh user has made me realize that we need: 1. Support for taking .xda/.xdr files which have been written with one setting of ENABLE_INFINITE_ELEMENTS and reading them with the opposite setting. Currently this seems to break. However, my recent experiences fixing the "scrambled DoF indices with nonLagrange AMR restarts" bug and the "not writing .xda files with enough precision" bug have made me realize the importance of: 2. Me never having to read through the xdr* source code again. I'd love to dive right in and implement my wishlist item 1, but unfortunately that conflicts with wishlist item 2. My hands are tied.  Roy 
From: Roy Stogner <roystgnr@ic...>  20070219 19:44:57

My current work with another libMesh user has made me realize that we need: 1. Support for taking .xda/.xdr files which have been written with one setting of ENABLE_INFINITE_ELEMENTS and reading them with the opposite setting. Currently this seems to break. However, my recent experiences fixing the "scrambled DoF indices with nonLagrange AMR restarts" bug and the "not writing .xda files with enough precision" bug have made me realize the importance of: 2. Me never having to read through the xdr* source code again. I'd love to dive right in and implement my wishlist item 1, but unfortunately that conflicts with wishlist item 2. My hands are tied.  Roy 
From: <tim@ce...>  20070220 08:13:49

Dear Roy, On Mon, 19 Feb 2007, Roy Stogner wrote: > My current work with another libMesh user has made me realize that we > need: > > 1. Support for taking .xda/.xdr files which have been written with > one setting of ENABLE_INFINITE_ELEMENTS and reading them with the > opposite setting. Currently this seems to break. > > However, my recent experiences fixing the "scrambled DoF indices with > nonLagrange AMR restarts" bug and the "not writing .xda files with > enough precision" bug have made me realize the importance of: > > 2. Me never having to read through the xdr* source code again. > > I'd love to dive right in and implement my wishlist item 1, but > unfortunately that conflicts with wishlist item 2. My hands are tied. Hmm. Just to confuse you, I wrote the attached test program. It creates a grid, solves a simple PDE on it, solves it twice again, then writes everything to xdr files, reads them in and solves again. Obviously (from the attached output, produced with ksp_monitor), reading the files in does not exactly restore the original state. However, I guess that changing this behaviour, even if it could become an item 3 of your wishlist, would still conflict with item 2. Best Regards, Tim 
From: John Peterson <peterson@cf...>  20070220 13:51:50

Hi Tim, If we write out only 12 digits of the solution (which I think is true for xda) could this explain the discrepancy? The residual norm computed from the readin solution looks the same to within about 12 digits. John Tim Kr,bv(Bger writes: > Dear Roy, > > On Mon, 19 Feb 2007, Roy Stogner wrote: > > > My current work with another libMesh user has made me realize that we > > need: > > > > 1. Support for taking .xda/.xdr files which have been written with > > one setting of ENABLE_INFINITE_ELEMENTS and reading them with the > > opposite setting. Currently this seems to break. > > > > However, my recent experiences fixing the "scrambled DoF indices with > > nonLagrange AMR restarts" bug and the "not writing .xda files with > > enough precision" bug have made me realize the importance of: > > > > 2. Me never having to read through the xdr* source code again. > > > > I'd love to dive right in and implement my wishlist item 1, but > > unfortunately that conflicts with wishlist item 2. My hands are tied. > > Hmm. Just to confuse you, I wrote the attached test program. It > creates a grid, solves a simple PDE on it, solves it twice again, then > writes everything to xdr files, reads them in and solves again. > Obviously (from the attached output, produced with ksp_monitor), > reading the files in does not exactly restore the original state. > However, I guess that changing this behaviour, even if it could become > an item 3 of your wishlist, would still conflict with item 2. > > Best Regards, > > Tim./test ksp_monitor > 0 KSP Residual norm 1.733328971787e+01 > 1 KSP Residual norm 1.185778941653e+01 > 2 KSP Residual norm 1.139134743178e+01 > 3 KSP Residual norm 9.913079902632e+00 > 4 KSP Residual norm 7.191145025915e+00 > 5 KSP Residual norm 5.172029883955e+00 > 6 KSP Residual norm 2.147425476520e+00 > 7 KSP Residual norm 8.687575748560e01 > 8 KSP Residual norm 4.164173667107e01 > 9 KSP Residual norm 2.443014766559e01 > 10 KSP Residual norm 1.888013690466e01 > 11 KSP Residual norm 1.545632671036e01 > 12 KSP Residual norm 9.547251214559e02 > 13 KSP Residual norm 3.558309380772e02 > 14 KSP Residual norm 2.119986479829e02 > 15 KSP Residual norm 1.472055648826e02 > 16 KSP Residual norm 7.844202526193e03 > 17 KSP Residual norm 4.324782781627e03 > 18 KSP Residual norm 2.537192913599e03 > 19 KSP Residual norm 1.372628313957e03 > 20 KSP Residual norm 7.072639404731e04 > 21 KSP Residual norm 2.981094657572e04 > 22 KSP Residual norm 1.807779624842e04 > 23 KSP Residual norm 1.340784951662e04 > 24 KSP Residual norm 1.088577505522e04 > 25 KSP Residual norm 7.207500151591e05 > 26 KSP Residual norm 3.887867139016e05 > 27 KSP Residual norm 2.190821106975e05 > 28 KSP Residual norm 1.040052415721e05 > 29 KSP Residual norm 3.856840518033e06 > 30 KSP Residual norm 1.989373498593e06 > 31 KSP Residual norm 1.442089144671e06 > 32 KSP Residual norm 8.398949810780e07 > 33 KSP Residual norm 4.613052610992e07 > 34 KSP Residual norm 2.446011262135e07 > 35 KSP Residual norm 1.402013920665e07 > 36 KSP Residual norm 9.288528753379e08 > 37 KSP Residual norm 5.569095323071e08 > 38 KSP Residual norm 3.080619009371e08 > 39 KSP Residual norm 2.149324393284e08 > 40 KSP Residual norm 1.668599996685e08 > 41 KSP Residual norm 1.297303487430e08 > 42 KSP Residual norm 9.127468731142e09 > 43 KSP Residual norm 6.183347158799e09 > 44 KSP Residual norm 4.688102496189e09 > 45 KSP Residual norm 4.418529129665e09 > 46 KSP Residual norm 4.266320304576e09 > 47 KSP Residual norm 4.097419313482e09 > 48 KSP Residual norm 3.918923888652e09 > 49 KSP Residual norm 3.395719439036e09 > 50 KSP Residual norm 2.538813185638e09 > 51 KSP Residual norm 1.415015667448e09 > 52 KSP Residual norm 6.308264312374e10 > 53 KSP Residual norm 3.843345277332e10 > 54 KSP Residual norm 2.514734506456e10 > 55 KSP Residual norm 1.710479832166e10 > 56 KSP Residual norm 1.151796274421e10 > 57 KSP Residual norm 7.714264551552e11 > 58 KSP Residual norm 4.742723602833e11 > 59 KSP Residual norm 3.379626986680e11 > 60 KSP Residual norm 2.548253614890e11 > 61 KSP Residual norm 2.095190091567e11 > 62 KSP Residual norm 1.512617662327e11 >  > 0 KSP Residual norm 1.510050418779e11 >  > 0 KSP Residual norm 1.510050418779e11 >  > 0 KSP Residual norm 1.611365184359e11 > #include <iostream> > #include <algorithm> > #include <math.h> > > #include "libmesh.h" > #include "mesh.h" > #include "mesh_generation.h" > #include "gmv_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 "linear_implicit_system.h" > #include "transient_system.h" > #include "perf_log.h" > #include "boundary_info.h" > #include "utility.h" > #include "elem.h" > #include "mesh_refinement.h" > > void assemble(EquationSystems& es, const std::string& system_name) > { > const Mesh& mesh = es.get_mesh(); > const unsigned int dim = mesh.mesh_dimension(); > LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("e"); > 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; > DenseVector<Number> Fe; > std::vector<unsigned int> dof_indices; > > MeshBase::const_element_iterator el = mesh.elements_begin(); > const MeshBase::const_element_iterator end_el = mesh.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]*(dphi[i][qp]*dphi[j][qp]+phi[i][qp]*phi[j][qp]); > } > for (unsigned int i=0; i<phi.size(); i++) > { > Fe(i) += JxW[qp]*1.0*phi[i][qp]; > } > } > > system.matrix>add_matrix (Ke, dof_indices); > system.rhs>add_vector (Fe, dof_indices); > } > } > > int main (int argc, char** argv) > { > libMesh::init (argc, argv); > { > const unsigned int dim = 3; > Mesh mesh (dim); > > MeshTools::Generation::build_cube (mesh, > 4, 4, 4, > 0., 1., > 0., 1., > 0., 1., > HEX27); > MeshRefinement mesh_refinement(mesh); > mesh_refinement.uniformly_refine(1); > EquationSystems equation_systems (mesh); > equation_systems.add_system<LinearImplicitSystem> ("e"); > equation_systems.get_system("e").add_variable("u", SECOND); > equation_systems.get_system("e").attach_assemble_function (assemble); > equation_systems.init(); > equation_systems.get_system("e").solve(); > > std::cout << "" << std::endl; > > equation_systems.get_system("e").solve(); > > std::cout << "" << std::endl; > > equation_systems.get_system("e").solve(); > > std::cout << "" << std::endl; > > mesh.write("mesh.xdr"); > equation_systems.write("systems.xdr", > libMeshEnums::WRITE, > EquationSystems::WRITE_DATAEquationSystems::WRITE_ADDITIONAL_DATA); > > Mesh mesh2(dim); > mesh2.read("mesh.xdr"); > EquationSystems equation_systems2(mesh2); > equation_systems2.read("systems.xdr", > libMeshEnums::READ, > EquationSystems::READ_DATAEquationSystems::READ_ADDITIONAL_DATAEquationSystems::READ_HEADER); > equation_systems2.reinit(); > equation_systems2.get_system("e").attach_assemble_function (assemble); > equation_systems2.get_system("e").solve(); > } > > return libMesh::close (); > } > >  > Take Surveys. Earn Cash. Influence the Future of IT > Join SourceForge.net's Techsay panel and you'll get the chance to share your > opinions on IT & business topics through brief surveysand earn cash > http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV_______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20070220 13:58:04

I am too lazy to do this myself  what if you write/read an XDR file instead? It should store a double as a true double and hopefully bypass this issue. Original Message From: libmeshdevelbounces@... [mailto:libmeshdevelbounces@...] On Behalf Of John Peterson Sent: Tuesday, February 20, 2007 7:51 AM To: Tim Kroger Cc: libmeshdevel@...; Roy Stogner Subject: Re: [Libmeshdevel] Two items on my wishlist Hi Tim, If we write out only 12 digits of the solution (which I think is true for xda) could this explain the discrepancy? The residual norm computed from the readin solution looks the same to within about 12 digits. John Tim Kr=1B,bvger writes: > Dear Roy, > > On Mon, 19 Feb 2007, Roy Stogner wrote: > > > My current work with another libMesh user has made me realize that we > > need: > > > > 1. Support for taking .xda/.xdr files which have been written with > > one setting of ENABLE_INFINITE_ELEMENTS and reading them with the > > opposite setting. Currently this seems to break. > > > > However, my recent experiences fixing the "scrambled DoF indices with > > nonLagrange AMR restarts" bug and the "not writing .xda files with > > enough precision" bug have made me realize the importance of: > > > > 2. Me never having to read through the xdr* source code again. > > > > I'd love to dive right in and implement my wishlist item 1, but > > unfortunately that conflicts with wishlist item 2. My hands are tied. > > Hmm. Just to confuse you, I wrote the attached test program. It > creates a grid, solves a simple PDE on it, solves it twice again, then > writes everything to xdr files, reads them in and solves again.=20 > Obviously (from the attached output, produced with ksp_monitor), > reading the files in does not exactly restore the original state.=20 > However, I guess that changing this behaviour, even if it could become > an item 3 of your wishlist, would still conflict with item 2. > > Best Regards, > > Tim./test ksp_monitor > 0 KSP Residual norm 1.733328971787e+01=20 > 1 KSP Residual norm 1.185778941653e+01=20 > 2 KSP Residual norm 1.139134743178e+01=20 > 3 KSP Residual norm 9.913079902632e+00=20 > 4 KSP Residual norm 7.191145025915e+00=20 > 5 KSP Residual norm 5.172029883955e+00=20 > 6 KSP Residual norm 2.147425476520e+00=20 > 7 KSP Residual norm 8.687575748560e01=20 > 8 KSP Residual norm 4.164173667107e01=20 > 9 KSP Residual norm 2.443014766559e01=20 > 10 KSP Residual norm 1.888013690466e01 > 11 KSP Residual norm 1.545632671036e01 > 12 KSP Residual norm 9.547251214559e02 > 13 KSP Residual norm 3.558309380772e02 > 14 KSP Residual norm 2.119986479829e02 > 15 KSP Residual norm 1.472055648826e02 > 16 KSP Residual norm 7.844202526193e03 > 17 KSP Residual norm 4.324782781627e03 > 18 KSP Residual norm 2.537192913599e03 > 19 KSP Residual norm 1.372628313957e03 > 20 KSP Residual norm 7.072639404731e04 > 21 KSP Residual norm 2.981094657572e04 > 22 KSP Residual norm 1.807779624842e04 > 23 KSP Residual norm 1.340784951662e04 > 24 KSP Residual norm 1.088577505522e04 > 25 KSP Residual norm 7.207500151591e05 > 26 KSP Residual norm 3.887867139016e05 > 27 KSP Residual norm 2.190821106975e05 > 28 KSP Residual norm 1.040052415721e05 > 29 KSP Residual norm 3.856840518033e06 > 30 KSP Residual norm 1.989373498593e06 > 31 KSP Residual norm 1.442089144671e06 > 32 KSP Residual norm 8.398949810780e07 > 33 KSP Residual norm 4.613052610992e07 > 34 KSP Residual norm 2.446011262135e07 > 35 KSP Residual norm 1.402013920665e07 > 36 KSP Residual norm 9.288528753379e08 > 37 KSP Residual norm 5.569095323071e08 > 38 KSP Residual norm 3.080619009371e08 > 39 KSP Residual norm 2.149324393284e08 > 40 KSP Residual norm 1.668599996685e08 > 41 KSP Residual norm 1.297303487430e08 > 42 KSP Residual norm 9.127468731142e09 > 43 KSP Residual norm 6.183347158799e09 > 44 KSP Residual norm 4.688102496189e09 > 45 KSP Residual norm 4.418529129665e09 > 46 KSP Residual norm 4.266320304576e09 > 47 KSP Residual norm 4.097419313482e09 > 48 KSP Residual norm 3.918923888652e09 > 49 KSP Residual norm 3.395719439036e09 > 50 KSP Residual norm 2.538813185638e09 > 51 KSP Residual norm 1.415015667448e09 > 52 KSP Residual norm 6.308264312374e10 > 53 KSP Residual norm 3.843345277332e10 > 54 KSP Residual norm 2.514734506456e10 > 55 KSP Residual norm 1.710479832166e10 > 56 KSP Residual norm 1.151796274421e10 > 57 KSP Residual norm 7.714264551552e11 > 58 KSP Residual norm 4.742723602833e11 > 59 KSP Residual norm 3.379626986680e11 > 60 KSP Residual norm 2.548253614890e11 > 61 KSP Residual norm 2.095190091567e11 > 62 KSP Residual norm 1.512617662327e11 >  > 0 KSP Residual norm 1.510050418779e11=20 >  > 0 KSP Residual norm 1.510050418779e11=20 >  > 0 KSP Residual norm 1.611365184359e11=20 > #include <iostream> > #include <algorithm> > #include <math.h> > > #include "libmesh.h" > #include "mesh.h" > #include "mesh_generation.h" > #include "gmv_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 "linear_implicit_system.h" > #include "transient_system.h" > #include "perf_log.h" > #include "boundary_info.h" > #include "utility.h" > #include "elem.h" > #include "mesh_refinement.h" > > void assemble(EquationSystems& es, const std::string& system_name) > { > const Mesh& mesh =3D es.get_mesh(); > const unsigned int dim =3D mesh.mesh_dimension(); > LinearImplicitSystem& system =3D = es.get_system<LinearImplicitSystem> ("e"); > const DofMap& dof_map =3D system.get_dof_map(); > FEType fe_type =3D 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 =3D fe>get_JxW(); > const std::vector<std::vector<Real> >& phi =3D fe>get_phi(); > const std::vector<std::vector<RealGradient> >& dphi =3D fe>get_dphi(); > DenseMatrix<Number> Ke; > DenseVector<Number> Fe; > std::vector<unsigned int> dof_indices; >=20 > MeshBase::const_element_iterator el =3D mesh.elements_begin(); > const MeshBase::const_element_iterator end_el =3D mesh.elements_end(); > for ( ; el !=3D end_el ; ++el) > { > const Elem* elem =3D *el; > dof_map.dof_indices (elem, dof_indices); > fe>reinit (elem); > Ke.resize (dof_indices.size(), > dof_indices.size()); > Fe.resize (dof_indices.size()); >=20 > for (unsigned int qp=3D0; qp<qrule.n_points(); qp++) > for (unsigned int i=3D0; i<phi.size(); i++) > { > for (unsigned int j=3D0; j<phi.size(); j++) > { > Ke(i,j) +=3D JxW[qp]*(dphi[i][qp]*dphi[j][qp]+phi[i][qp]*phi[j][qp]); > } > for (unsigned int i=3D0; i<phi.size(); i++) > { > Fe(i) +=3D JxW[qp]*1.0*phi[i][qp]; > } > }=20 > =20 > system.matrix>add_matrix (Ke, dof_indices); > system.rhs>add_vector (Fe, dof_indices); > } > } > > int main (int argc, char** argv) > { > libMesh::init (argc, argv); > { =20 > const unsigned int dim =3D 3; =20 > Mesh mesh (dim); > =20 > MeshTools::Generation::build_cube (mesh, > 4, 4, 4, > 0., 1., > 0., 1., > 0., 1., > HEX27); > MeshRefinement mesh_refinement(mesh); > mesh_refinement.uniformly_refine(1); > EquationSystems equation_systems (mesh); > equation_systems.add_system<LinearImplicitSystem> ("e"); > equation_systems.get_system("e").add_variable("u", SECOND); > equation_systems.get_system("e").attach_assemble_function (assemble); > equation_systems.init(); > equation_systems.get_system("e").solve(); >=20 > std::cout << "" << std::endl; >=20 > equation_systems.get_system("e").solve(); >=20 > std::cout << "" << std::endl; >=20 > equation_systems.get_system("e").solve(); >=20 > std::cout << "" << std::endl; >=20 > mesh.write("mesh.xdr"); > equation_systems.write("systems.xdr", > libMeshEnums::WRITE, > EquationSystems::WRITE_DATAEquationSystems::WRITE_ADDITIONAL_DATA); >=20 > Mesh mesh2(dim); > mesh2.read("mesh.xdr"); > EquationSystems equation_systems2(mesh2); > equation_systems2.read("systems.xdr", > libMeshEnums::READ, > EquationSystems::READ_DATAEquationSystems::READ_ADDITIONAL_DATAEquatio nSystems::READ_HEADER); > equation_systems2.reinit(); > equation_systems2.get_system("e").attach_assemble_function (assemble); > equation_systems2.get_system("e").solve(); > } > =20 > return libMesh::close (); > } > >   > Take Surveys. Earn Cash. Influence the Future of IT > Join SourceForge.net's Techsay panel and you'll get the chance to share your > opinions on IT & business topics through brief surveysand earn cash > http://www.techsay.com/default.php?page=3Djoin.php&p=3Dsourceforge&CID=3D= DEVDE V_______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
From: <tim@ce...>  20070220 14:02:59

Dear John, On Tue, 20 Feb 2007, John Peterson wrote: > If we write out only 12 digits of the solution (which I think is true > for xda) could this explain the discrepancy? Yes, but I'm using xdr, not xda (see source code). Best Regards, Tim 
From: Roy Stogner <roystgnr@ic...>  20070220 14:01:39

On Tue, 20 Feb 2007, Tim Kröger wrote: > Hmm. Just to confuse you, I wrote the attached test program. It creates a > grid, solves a simple PDE on it, solves it twice again, then writes > everything to xdr files, reads them in and solves again. Obviously (from the > attached output, produced with ksp_monitor), reading the files in does not > exactly restore the original state. Now that's a tough one. On the one hand, the original state should be restored exactly in your example... I suppose there could be come data that's sitting in cache at 80 bits but gets reread at 64 bits, but I doubt it. On the other hand, I'm finding it hard to get too excited about a state change that's smaller than most linear solver tolerances. The gripping hand is it's not a .xda/.xdr problem at all. Uncommenting equation_systems2.reinit() causes the state to remain unchanged, whereas a call to equation_systems.reinit() causes the same sort of slight state change before any files are even written.  Roy 
From: Roy Stogner <roystgnr@ic...>  20070220 14:26:57

On Tue, 20 Feb 2007, Roy Stogner wrote: > The gripping hand is it's not a .xda/.xdr problem at all. > Uncommenting equation_systems2.reinit() causes the state to remain > unchanged, whereas a call to equation_systems.reinit() causes the > same sort of slight state change before any files are even written. The change appears to be occuring in System::project_vector  but with no elements flagged for refinement or coarsening, that function should just be copying degrees of freedom, not doing anything to change them. I can't figure out yet exactly where anything's being changed  DofMap::enforce_constraints_exactly() gets called at the end of the projection now, but that shouldn't have any effect on a uniform mesh.  Roy 