<div class="gmail_quote">On Sat, Aug 13, 2011 at 19:02, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<br>
Jed,<br>
<br>
You will need to explain to me all this code or I am going to have to rip it all out since Matt wants a release.<br>
<br>
That is, all the DMDA interpolation stuff that was carefully put in under your supervision will be removed unless I understand why it is there. I asked you once before and you didn't really answer.<br>
<br>
For example in the 3d interpolation it has a<br>
<br>
if (!vcoors) {<br>
<br>
then two blocks of code NEITHER of which actually use the coordinate information AT ALL. So why the two blocks of code? Is one wrong and one right and if so why not just get rid of the wrong one.</blockquote><div><br>
</div><div>Dave showed some cases where the old code was producing incorrect results, and then wrote a bunch of new code to get second-order accurate interpolants. I thought the old code was also messy, so I didn't carefully audit the new code that came in passing the existing tests and adding new ones. I asked Dave about the code duplication and he wrote the following just before I fell off the internet for a couple days and then forgot to follow up. :-/</div>
<div><br></div><div><div><div style="font-style: italic; ">1) We never wanted to support non-uniform grids. This would have</div><div style="font-style: italic; ">mandated we also introduce</div><div style="font-style: italic; ">
a point location routine.</div><div style="font-style: italic; ">These routines were not meant to be completely general nor work with</div><div style="font-style: italic; ">any old set of DA's.</div><div style="font-style: italic; ">
We decided (and thought every one agreed around the time of last AGU)</div><div style="font-style: italic; ">that it was fine to</div><div style="font-style: italic; ">provide an interpolation routine which would work only if the DA's</div>
<div style="font-style: italic; ">were created using DM{Refine,Coarsen}. Okay... the periodic case</div><div style="font-style: italic; ">wasn't supported....</div><div style="font-style: italic; "><br></div><div>This is true, DMRefine always does regular refinement so we don't really have a source for refinement that is not isotropic on the reference element. If we assume an isoparametrically mapped basis, then we don't need to look at coordinates at all to do the interpolation (Dave also makes this point below).</div>
<div style="font-style: italic; "><br></div><div style="font-style: italic; ">2) I was always worried about the periodic case, but there were not</div><div style="font-style: italic; ">any tests availble</div><div style="font-style: italic; ">
to verify our implementation worked/didn't work with periodic (and</div><div style="font-style: italic; ">i've never used periodic</div><div style="font-style: italic; ">DMDA's and so I didn't feel super comfortable writing a new test)</div>
</div><div style="font-style: italic; "><br></div><div>Whatever code is in there needs to work for periodic, even if it's only for Cartesian meshes. Without making risky assumptions, we do not currently have enough information to handle non-Cartesian periodic meshes because the thickness of the "wrap" element cannot be determined from the global coordinates. We should at least agree on some standard assumption here so that the code doesn't just crash if you set coordinates.</div>
<div style="font-style: italic; "><br></div><div><i><br></i></div><div><i>3) The if (!vcoors) {} block could be completely removed and we just</i></div><div><i>leave the other block, however I left it in there so that the function</i></div>
<div><i>would behave as the old version if no coordinates were defined.</i></div></div><div><i><br></i></div><div>This just makes things confusing, especially considering that the branch taken when vcoords is set doesn't even use them. The (!vcoords) branch is currently the only one that works for periodic, and it's the only one being taken (!) since</div>
<div><br></div><div><div>changeset: 19720:f6d8213d8f85</div><div>user: Jungho Lee <a href="mailto:julee@mcs.anl.gov">julee@mcs.anl.gov</a></div><div>date: Wed Jul 27 14:49:35 2011 -0500</div><div>summary: Putting in a temporary fix which seems to handle periodic boundary conditions correctly.</div>
<div><br></div><div>diff --git a/src/dm/impls/da/dainterp.c b/src/dm/impls/da/dainterp.c</div><div>--- a/src/dm/impls/da/dainterp.c</div><div>+++ b/src/dm/impls/da/dainterp.c</div><div>@@ -364,7 +364,8 @@</div><div> }</div>
<div> </div><div> /* loop over local fine grid nodes setting interpolation for those*/</div><div>- if (!vcoors) {</div><div>+ // if (!vcoors) {</div><div>+ if (1) {</div><div> </div><div> for (j=j_start; j<j_start+n_f; j++) {</div>
<div> for (i=i_start; i<i_start+m_f; i++) {</div></div><div><br></div><div><br></div><div>As far as I can tell, the ex36 tests are passing this way. Hopefully Dave can dig out a test case where it was doing something incorrect, because I doubt he wrote all that code to fix something that wasn't broken.</div>
<div><i><br></i></div><div><i><div>4) Even through the global coordinates are not directly used, the</div><div>interpolant is still valid.</div><div>We assume coordinates are interpolated across an element via</div><div> x = \sum N_i ( \xi ) x*_i</div>
<div>where N_i is the usual linear, bilinear, trilinear basis, \xi is the</div><div>local coord and x*_i is the nodal value of x</div><div>The restriction operator is then</div><div> R_ij = N_i^C( \xi_j^F )</div><div>where N_i^C is the basis on the coarse grid and \xi_j^F is the local</div>
<div>coordinate on the fine grid.</div><div>In DMDAGetInterpolation() we obviously assemble the transpose of this.</div></i></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
BTW the second version does not work for periodic conditions which is bad!</blockquote><div><br></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Is your intention to use the coordinate information on each level to generate an interpolation matrix that depends on the coordinates or not? Why is the vcoors flag used as a check here? Are you assuming that the vcoors existing has some special meaning? I understand that the code<br>
<br>
Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);<br>
Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);<br>
Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);<br>
Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);<br>
<br>
Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);<br>
Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);<br>
Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);<br>
Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);<br>
<br>
comes from the trilinear finite element basic functions to define an interpolation but since you are evaluating them on a uniform mesh<br>
<br>
for (li=0; li<nxi; li++) {<br>
xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));<br>
}<br>
for (lj=0; lj<neta; lj++) {<br>
eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));<br>
}<br>
for (lk=0; lk<nzeta; lk++) {<br>
zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));<br>
}<br>
<br>
why bother? Isn't this always going to give you a uniform grid interpolation? Also src/dm/examples/tests/ex36.c seems to be suppose to test the interpolation for nonuniform grids but it doesn't seem to work<br>
<br>
arry-smiths-macbook-pro:tests barrysmith$ ./ex36 -dim 2 -nl 1 -cmap 1<br>
[2 x 2]=>[4 x 4], interp err = 1.4873e+00<br>
<br>
and doesn't have any "nightly" output files so isn't a test at all. Is there a test that is actually run in the nightly for this?<br></blockquote><div><br></div><div>Did you try make runex36_*? They should show second order convergence. If you look at the makefile, these have been running in the nightlies since Dec 1, 2010, so I think they are working.</div>
<div><br></div><div>We need more tests, something for periodic and something showing why the (!vcoords) code is worse than the other version.</div></div>