<div dir="ltr"><div><div>For your pressue Poisson solve, you can use DMDAGetReducedDMDA to generate a DMDA with only 1 dof. This reduced DMDA will have the same parallel layout as your DMDA with 4 dofs.<br><br>See <br><a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMDAGetReducedDMDA.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMDAGetReducedDMDA.html</a><br><br></div>Cheers,<br></div>  Dave<br></div><div class="gmail_extra"><br><div class="gmail_quote">On 15 October 2014 15:41, Håkon Strandenes <span dir="ltr"><<a href="mailto:haakon@hakostra.net" target="_blank">haakon@hakostra.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Thanks for the reply, both Jed and Matthew. However, I am still a bit unsure on a small detail:<br>
<br>
I have created a 3D DMDA with ndof = 4, one for each velocity component (u,v,w) and one for the pressure (p). I only solve a Poisson equation to correct p, the velocities are solved explicitly.<br>
<br>
My question is this time:<br>
When I create a matrix with DMCreateMatrix(), I get four times as many rows as I really need, i.e. it assumes that I want to solve equations for all my DOF's. The same happens when I create a vector with DMCreateGlobalVector(), it gets one element for each dof. This is not desirable, since I only solve equations for one unknown per gridpoint. If my domain (for debugginf purposes) have 250 points, I want the local portion of my matrix to have 250 rows. Instead it gets 1000.<br>
<br>
Was I wrong when I created my DMDA with ndof = 4? Should I have created one DMDA for the velocity with ndof = 3 and one with ndof = 1 to use as a basis for the pressure solution? Or have I again missed something fundamentally?<br>
<br>
Regards,<br>
Håkon<div class="HOEnZb"><div class="h5"><br>
<br>
<br>
On 08. okt. 2014 22:07, Jed Brown wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Håkon Strandenes <<a href="mailto:haakon@hakostra.net" target="_blank">haakon@hakostra.net</a>> writes:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi,<br>
<br>
I am making a new, simple incompressible CFD code. I have created a grid<br>
as a 3D DMDA, with ndof=4, one dof for each velocity component (u,v,w)<br>
and one for the pressure (p). Currently I am creating a forward Euler<br>
step, where the velocities are updated based on the velocities and<br>
pressure from the previous time step. That works great. Now I need to<br>
solve a Poisson equation for the pressure correction term, enforcing<br>
continuity onto my solution, and this is where I have a few questions<br>
for you:<br>
<br>
1)<br>
How do I set up my equation system in the most efficient way? The<br>
vectors from the DMDA is not contiguous in a global sense, i.e. if I map<br>
the local (i,j,k) indices to a global index with<br>
<br>
globId = (<a href="http://info.mx" target="_blank">info.mx</a>*info.my)*k + (info-.mx)*j + i;<br>
</blockquote>
<br>
This is a "natural index", used in some cases for IO, but not for<br>
computation.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
globId is not gontignous on any process. How do I create my coefficient<br>
matrix, r.h.s. vector and solution vector?<br>
</blockquote>
<br>
DMCreateGlobalVector()<br>
<br>
it will use what is described in the Users Manual as "PETSc ordering".<br>
Please read this section of the manual; I think it will answer your<br>
questions.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Is there some way I can make the local parts of the matrix and vectors<br>
mimic the decomposition from the DMDA? Or should I push all values<br>
into a different decomposition scheme, where each process holds one<br>
contigous slab of the equation system using the global ID's calculated<br>
above? Have I missed some really basic stuff here?<br>
<br>
2)<br>
I know that my coefficient matrix from the Poisson equation is<br>
symmetric. Currently I am still creating and inserting values on both<br>
the upper and lower diagonal. Is this really necessary? Or is there a<br>
"smart" way I can save half the storage by only creating the upper or<br>
lower part?<br>
</blockquote>
<br>
-dm_mat_type sbaij<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
3)<br>
I solve a transient problem, but my coefficient matrix is the same for<br>
all time steps and I only recompute the r.h.s. vector for each time<br>
step. I have got the impression that KSPsolve "finds out" if the A<br>
matrix have changed since last invocation or not, recalculating the<br>
preconditioner as necessary.<br>
</blockquote>
<br>
If you just call KSPSolve, it will reuse everything.  Call<br>
KSPSetOperators (assuming you are manually doing the time<br>
discretization) if the matrix changes.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
But still I feel that I should utilize this property more. Do you have<br>
any tips in how I can use this in the most efficient way possible? For<br>
example really, really efficient, but expensive, preconditioners?<br>
<br>
Thanks in advance.<br>
<br>
Regards,<br>
Håkon Strandenes<br>
</blockquote></blockquote>
</div></div></blockquote></div><br></div>