<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 26, 2014 at 10:13 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><span class="">On Tue, Nov 25, 2014 at 6:40 PM, 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:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div><div><div class="gmail_quote">On Tue, Nov 25, 2014 at 11:36 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Dear all,<br>I'm solving a system using petsc and I get a global Vec, say X . It uses DMDA with 4 dofs (3 velocity + 1 pressure). X contains velocity at the cell faces since I solve using staggered grid.<br>Now I'd like to create one array for velocity with values at cell centers and another array for pressure (not the Vec so that I can send the pointer to the array to other part of the code that does not use Petsc).<br><br></div><div>Currently what I do is :<br>------ Vec X, X_local;<br>------ PetscScalar *X_array;<br><br>// Scatter X to X_local and then use:<br><br>------- VecGetArray(X_local, &X_array) <br></div><div><br>And then have a function, say <br>getVelocityAt(x, y, z, component) {<br></div><div>// interpolate velocity at (x,y,z) cell center using X_array<br></div><div>}<br><br>The function getVelocityAt() gets called from outside petsc in a loop over all (x,y,z) positions.<br></div><div>This is not done in parallel.<br></div><div><br>Now, how do I use Petsc to instead interpolate the cell center velocities in parallel and store it<br>in an array say<br>PetscScalar *X_array_cellCenter;<br>?<br></div><div>This would need to have size one less along each axis compared to the orginal DMDA size.<br></div><div>This way I intend to return X_array_cellCenter to the code outside Petsc.<br></div></div>
</blockquote></div><br></div></div>SNES ex30 is an example of a staggered grid code using DMDA. It does this kind of interpolation,</div><div class="gmail_extra">and puts the result in a Vec.</div></div></blockquote></span><div><br>After looking at the example, I tried the following but I have a problem:<br><br>The idea I take from the e.g is to use DMDALocalInfo *info to get local<br>grid coordinates and use for loop to iterate through local grid values that are<br>available in the Field **x array to do the interpolation.<br><br>Here is what I tried in my case where I get the solution from a ksp attached to a dmda:<br><br>// Get the solution to Vec x from ksp. <br>// ksp is attached to DMDA da which has 4 dofs: vx,vy,vz,p;<br>KSPGetSoution(ksp, &x);<br>//x has velocity solution in faces of the cell.<br><br>//Get the array from x:<br>Field ***xArray; //Field is a struct with 4 members: vx,vy,vz and p.<br>DMDAGetVectorArray(da,xVec,&xArray); <br><br>//Now I want to have a new array to store the cell center velocity values<br>//by interpolating from xArray.<br><br>DMCreateGlobalVector(daV, &xC); //daV has same size as da but with dof=3<br>FieldVelocity ***xCArray; //FieldVelocity is a struct with 3 members: vx,vy and vz.<br>DMDAVecGetArray(daV,xC,&xCArray);<br><br>//Do the interpolation<br>DMDAGetLocalInfo(da,&info);<br>for (PetscInt k=info.zs; k<info.zs+<a href="http://info.zm" target="_blank">info.zm</a>; ++k) {<br> for (PetscInt j=info.ys; j<info.ys+info.ym; ++j) {<br> for (PetscInt i=info.xs; i<info.xs+info.xm; ++i) {<br> //do interpolation which requires access such as:<br> if(i==0 || j==0 || k==0) {<br> xCArray[k][j][i].vx = 0; //boundary condition;<br> xCArray[k][j][i].vy = 0;<br> xCArray[k][j][i].vz = 0<br> } else {<br> xCArray[k][j][i].vx = 0.5*(xArray[k][j][i].vx + xArray[k][j][i-1].vx);<br> xCArray[k][j][i].vy = 0.5*(xArray[k][j][i].vy + xArray[k][j][i-1].vy);<br> xCArray[k][j][i].vz = 0.5*(xArray[k][j][i].vz + xArray[k][j][i-1].vz);<br> }<br><br>I get runtime error within this loop. <br>What am I doing wrong here ?<br><br></div></div></div></div></blockquote><div><br></div><div>Ok, I now see that KSPGetSolution(ksp, &x) gives a GLOBAL vector in x and not a local vector; hence the runtime memory error above since I can't access ghost values in indices such as [i-1] wth the global vec.<br></div><div>Was there any reason not to make this information explicit in the docs to prevent the confusion I had ? Or should it have been obvious to me and that I'm missing something here ? <br><a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPGetSolution.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPGetSolution.html</a> <br></div><div>Thanks <br></div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div> <br></div><span class=""><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><br></div><div class="gmail_extra"> Thanks,</div><div class="gmail_extra"><br></div><div class="gmail_extra"> Matt<span><font color="#888888"><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></span></div><br></div></div>
</blockquote></div><br></div></div>