<div dir="ltr"><div><div><div><div><div><div><div><div><div>Matt,<br><br></div>I have attempted this but I am getting an error when trying to output the vector into HDF5. This is my implementation:<br><br></div><div>/* Initialize */<br></div><div>PetscDS prob;<br></div><div>Vec vlocal, ulocal, vglobal;<br></div>void (*vX[1])(....) = {velx};<br><br></div><div>/* Create DM for velocity */<br></div><div>ierr = DMClone(dm,&user.dmVel);CHKERRQ(ierr);<br>ierr = DMPlexCopyCoordinates(dm, user.dmVel);CHKERRQ(ierr);<br>ierr = PetscFECreateDefault(dm, spatialDim, 1, PETSC_TRUE, "solution_", -1, &user.feVel);CHKERRQ(ierr);<br>ierr = DMGetDS(user.dmVel, &prob);CHKERRQ(ierr);<br>ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) user.feVel);CHKERRQ(ierr);<br>ierr = DMCreateGlobalVector(user.dmVel,&vglobal);CHKERRQ(ierr);<br>ierr = PetscObjectSetName((PetscObject) vglobal, "velocity");CHKERRQ(ierr);<br></div><br></div>/* Obtain velocity */<br>ierr = DMGetLocalVector(dm,&ulocal);CHKERRQ(ierr);<br>ierr = VecDuplicate(ulocal,&vlocal);CHKERRQ(ierr);<br>ierr = DMGlobalToLocalBegin(dm,user.u,INSERT_ALL_VALUES,ulocal);CHKERRQ(ierr);<br>ierr = DMGlobalToLocalEnd(dm,user.u,INSERT_ALL_VALUES,ulocal);CHKERRQ(ierr);<br>ierr = DMPlexProjectFieldLocal(dm,ulocal,vX,INSERT_ALL_VALUES,vlocal);CHKERRQ(ierr);<br>ierr = DMLocalToGlobalBegin(dm,vlocal,INSERT_ALL_VALUES,vglobal);CHKERRQ(ierr);<br>ierr = DMLocalToGlobalEnd(dm,vlocal,INSERT_ALL_VALUES,vglobal);CHKERRQ(ierr);<br>ierr = VecViewFromOptions(vglobal,NULL,"-vec_view_x");CHKERRQ(ierr);<br><br></div>And as soon as I get to the last line, it gives me this error:<br><br><br>[0]PETSC ERROR: PetscMallocValidate: error detected at PetscViewerFileSetName_HDF5() line 251 in /home/justin/petsc-dev/src/sys/classes/viewer/impls/hdf5/hdf5v.c<br>[0]PETSC ERROR: Memory at address 0x3077ac0 is corrupted<br>[0]PETSC ERROR: Probably write past beginning or end of array<br>[0]PETSC ERROR: Last intact block allocated in VecCreate_Seq() line 38 in /home/justin/petsc-dev/src/vec/vec/impls/seq/bvec3.c<br>[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>[0]PETSC ERROR: Memory corruption: <a href="http://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind">http://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind</a><br>[0]PETSC ERROR: <br>[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>[0]PETSC ERROR: Petsc Development GIT revision: v3.5.3-2766-g0d12e5d GIT Date: 2015-04-21 21:32:23 -0500<br>[0]PETSC ERROR: ./main on a arch-linux2-c-debug named pacotaco by justin Thu Apr 23 22:14:21 2015<br>[0]PETSC ERROR: Configure options --download-chaco --download-exodusii --download-fblaslapack --download-hdf5 --download-metis --download-mpich --download-mumps --download-netcdf --download-parmetis --download-scalapack --download-tetgen --download-triangle --with-cc=gcc --with-clanguage=cxx --with-cmake=cmake --with-ctetgen --with-cxx=g++ --with-debugging=1 --with-fc=gfortran --with-valgrind=1 PETSC_ARCH=arch-linux2-c-debug<br>[0]PETSC ERROR: #1 PetscMallocValidate() line 135 in /home/justin/petsc-dev/src/sys/memory/mtr.c<br>[0]PETSC ERROR: #2 PetscViewerFileSetName_HDF5() line 251 in /home/justin/petsc-dev/src/sys/classes/viewer/impls/hdf5/hdf5v.c<br>[0]PETSC ERROR: #3 PetscViewerFileSetName() line 623 in /home/justin/petsc-dev/src/sys/classes/viewer/impls/ascii/filev.c<br>[0]PETSC ERROR: #4 PetscOptionsGetViewer() line 130 in /home/justin/petsc-dev/src/sys/classes/viewer/interface/viewreg.c<br>[0]PETSC ERROR: #5 PetscObjectViewFromOptions() line 2622 in /home/justin/petsc-dev/src/sys/objects/options.c<br>[0]PETSC ERROR: #6 main() line 1237 in /home/justin/Dropbox/dmplex-nonneg/main.c<br>[0]PETSC ERROR: PETSc Option Table entries:<br>[0]PETSC ERROR: -al 100<br>[0]PETSC ERROR: -at 0.01<br>[0]PETSC ERROR: -bcloc 496000,496000,536000,542000,0,100,503000,503000,536000,542000,0,100,496000,503000,536000,536000,0,100,496000,503000,542000,542000,0,100<br>[0]PETSC ERROR: -bcnum 4<br>[0]PETSC ERROR: -bcval 1,1,1,1<br>[0]PETSC ERROR: -dim 3<br>[0]PETSC ERROR: -dm_refine 0<br>[0]PETSC ERROR: -dt 0.01<br>[0]PETSC ERROR: -edges 3,3,3<br>[0]PETSC ERROR: -floc 499550,499600,538950,539000,0,100<br>[0]PETSC ERROR: -flow datafiles/chrom_flow.dat<br>[0]PETSC ERROR: -fnum 1<br>[0]PETSC ERROR: -ftime 0,99<br>[0]PETSC ERROR: -fval -0.1<br>[0]PETSC ERROR: -ksp_rtol 1.0e-9<br>[0]PETSC ERROR: -ksp_type cg<br>[0]PETSC ERROR: -lower 0,0,0<br>[0]PETSC ERROR: -mat_petscspace_order 0<br>[0]PETSC ERROR: -mesh datafiles/chrom_mesh.dat<br>[0]PETSC ERROR: -mu 3.95e-5<br>[0]PETSC ERROR: -nonneg 0<br>[0]PETSC ERROR: -numsteps 0<br>[0]PETSC ERROR: -options_left 0<br>[0]PETSC ERROR: -pc_type jacobi<br>[0]PETSC ERROR: -progress 1<br>[0]PETSC ERROR: -solution_petscspace_order 1<br>[0]PETSC ERROR: -tao_type tron<br>[0]PETSC ERROR: -upper 1,1,1<br>[0]PETSC ERROR: -vec_view_x hdf5:solx.h5<br>[0]PETSC ERROR: -vec_view_y hdf5:soly.h5<br>[0]PETSC ERROR: -vec_view_z hdf5:solz.h5<br>[0]PETSC ERROR: -vtuname pressure<br>[0]PETSC ERROR: -vtuprintat 1,20,80,160,320,640,999<br>[0]PETSC ERROR: -vtusteps 7<br>[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<br><br></div>Any idea what's going on here? Also, if I replace the last line with VecView(vglobal,PETSC_VIEWER_STDOUT_WORLD) I get a completely wacko error:<br><br>Vec Object:@����r?������b?�jүr?�.<br> ~E�b?�\2p}^r?���"�.b?'xDs��q?�|n��a?`K�Lwq?�n��a?�Ψ���p?�g◻u`?!�<br> |�<br> p?;�:2_?����5n?�6/w $]?��I!}�k?!c/�Z?�s#li?�=Y��X?I%�[�f?^�u=��T?�Sx��Yc?�F�u�Q?6<br> e��_?;6��a3L?z_�RX?���l�D?X�(�(�P?=�/k�9?;z�IB?<br> R�ٕp$?�j�-�?��7��-��H>�h�:����6 5�{��2@�L��]P�I#B�<U��U�O�fV`oI�7trdj�\��\��w0P��)e���a�Kf�s�uS�K���ze�����V�ΏY�g��3�2WY��Xëj�)�&M�[�s�!Q&+m��*�l�V^�Ӈй�uo��ڄD`�z 1 MPI processes<br> type: seq<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range<br>[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger<br>[0]PETSC ERROR: or see <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind</a><br>[0]PETSC ERROR: or try <a href="http://valgrind.org">http://valgrind.org</a> on GNU/linux and Apple Mac OS X to find memory corruption errors<br>[0]PETSC ERROR: PetscMallocValidate: error detected at PetscSignalHandlerDefault() line 149 in /home/justin/petsc-dev/src/sys/error/signal.c<br>[0]PETSC ERROR: Memory at address 0x2575ac0 is corrupted<br>[0]PETSC ERROR: Probably write past beginning or end of array<br>[0]PETSC ERROR: Last intact block allocated in VecCreate_Seq() line 38 in /home/justin/petsc-dev/src/vec/vec/impls/seq/bvec3.c<br>[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>[0]PETSC ERROR: Memory corruption: <a href="http://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind">http://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind</a><br>[0]PETSC ERROR: <br>[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>[0]PETSC ERROR: Petsc Development GIT revision: v3.5.3-2766-g0d12e5d GIT Date: 2015-04-21 21:32:23 -0500<br>[0]PETSC ERROR: ./main on a arch-linux2-c-debug named pacotaco by justin Thu Apr 23 22:17:26 2015<br>[0]PETSC ERROR: Configure options --download-chaco --download-exodusii --download-fblaslapack --download-hdf5 --download-metis --download-mpich --download-mumps --download-netcdf --download-parmetis --download-scalapack --download-tetgen --download-triangle --with-cc=gcc --with-clanguage=cxx --with-cmake=cmake --with-ctetgen --with-cxx=g++ --with-debugging=1 --with-fc=gfortran --with-valgrind=1 PETSC_ARCH=arch-linux2-c-debug<br>[0]PETSC ERROR: #1 PetscMallocValidate() line 135 in /home/justin/petsc-dev/src/sys/memory/mtr.c<br>[0]PETSC ERROR: #2 PetscSignalHandlerDefault() line 149 in /home/justin/petsc-dev/src/sys/error/signal.c<br><br><br></div>However, if I run the program with debugging off (i.e., --with-debugging=0), VecView gives me an answer that seems relatively correct.<br><br></div>Any idea what's going on here?<br><br></div>Thanks,<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Apr 23, 2015 at 4:00 AM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><div class="h5">On Wed, Apr 22, 2015 at 9:57 PM, Justin Chang <span dir="ltr"><<a href="mailto:jchang27@uh.edu" target="_blank">jchang27@uh.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div>Hello,<br><br></div>I am using DMPlexProjectField(...) to postprocess my final solution u into another vector sol. I have the following pointwise functions:<br><br></div><div>/* x-component of gradient */<br></div><div>void velx(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar v[])<br>{<br> v[0] = -a[0]*mu*u_x[0];<br>}<br><br>/* y-component of gradient */<br>void vely(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar v[])<br>{<br> v[0] = -a[0]*mu*u_x[1];<br>}<br><br>/* z-component of gradient */<br>void velz(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar v[])<br>{<br> v[0] = -a[0]*mu*u_x[2];<br>}<br><br></div>Where a[0] is the cell-wise permeability and mu is the (inverse of) viscosity. I am solving the diffusion equation (P1 elements). At the end of my simulation I have these lines:<br><br>ierr = VecDuplicate(user.u,&user.sol);CHKERRQ(ierr);<br>ierr = PetscObjectSetName((PetscObject) user.sol, "velocity");CHKERRQ(ierr);<br>void (*vX[1])(const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscReal *, PetscScalar *) = {velx};<br>ierr = DMPlexProjectField(dm,user.u,vX,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);<br>ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_x");CHKERRQ(ierr);<br><br>void (*vY[1])(const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscReal *, PetscScalar *) = {vely};<br>ierr = DMPlexProjectField(dm,user.u,vY,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);<br>ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_y");CHKERRQ(ierr);<br><br>if (spatialDim == 3) {<br> void (*vZ[1])(const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscReal *, PetscScalar *) = {velz};<br> ierr = DMPlexProjectField(dm,user.u,vZ,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);<br> ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_z");CHKERRQ(ierr);<br>}<br>ierr = VecDestroy(&user.sol);CHKERRQ(ierr);<br><br></div>But when i look at the results, the dirichlet constrained values of u were projected into the respective locations in sol. I thought putting INSERT_ALL_VALUES would insert the pointwise functions into all locations including the constrained ones. How can I project the above pointwise functions into the constrained fields? Or what am I potentially missing here?<br></div></div></blockquote><div><br></div></div></div><div>There are no Dirichlet constrained values in the global vector user.sol. When you call VecView() it gets projected to</div><div>an enlarged vector including the BC and output. What you want is to project that vector into a space with no constraints,</div><div>so you really need another DM/Section combination.</div><div><br></div><div>I think this will work:</div><div><br></div><div> 1) Get a local vector from dm</div><div><br></div><div> 2) Use DMPlexProjectFieldLocal(INSERT_ALL_VALUES) to fill it up</div><div><br></div><div> 3) DMClone(dm) and give it a Section with no constraints</div><div><br></div><div> 4) Use DMLocalToGlobal() into the new DM</div><div><br></div><div> 5) VecView() that new global vector</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div></div>Thanks,<span><font color="#888888"><br><br clear="all"><br>-- <br><div><div dir="ltr"><div><div><div>Justin Chang<br></div>PhD Candidate, Civil Engineering - Computational Sciences<br></div>University of Houston, Department of Civil and Environmental Engineering<br></div>Houston, TX 77004<br><a href="tel:%28512%29%20963-3262" value="+15129633262" target="_blank">(512) 963-3262</a><br></div></div></font></span></div>
</blockquote></span></div><span class="HOEnZb"><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</font></span></div></div>
</blockquote></div><br></div><br clear="all"><br>-- <br><div class="gmail_signature"><div dir="ltr"><div><div><div>Justin Chang<br></div>PhD Candidate, Civil Engineering - Computational Sciences<br></div>University of Houston, Department of Civil and Environmental Engineering<br></div>Houston, TX 77004<br>(512) 963-3262<br></div></div>