[petsc-dev] Jed will not be happy

[Jed Brown]

On Sat, Aug 13, 2011 at 19:02, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>   Jed,
>
>    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.
>
>    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.
>
>    For example in the 3d interpolation it has a
>
>  if (!vcoors) {
>
>    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.


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. :-/

1) We never wanted to support non-uniform grids. This would have
mandated we also introduce
a point location routine.
These routines were not meant to be completely general nor work with
any old set of DA's.
We decided (and thought every one agreed around the time of last AGU)
that it was fine to
provide an interpolation routine which would work only if the DA's
were created using DM{Refine,Coarsen}. Okay... the periodic case
wasn't supported....

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).

2) I was always worried about the periodic case, but there were not
any tests availble
to verify our implementation worked/didn't work with periodic (and
i've never used periodic
DMDA's and so I didn't feel super comfortable writing a new test)

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.

*
*
*3) The if (!vcoors) {} block could be completely removed and we just*
*leave the other block, however I left it in there so that the function*
*would behave as the old version if no coordinates were defined.*
*
*
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

changeset:   19720:f6d8213d8f85
user:        Jungho Lee julee at mcs.anl.gov
date:        Wed Jul 27 14:49:35 2011 -0500
summary:     Putting in a temporary fix which seems to handle periodic
boundary conditions correctly.

diff --git a/src/dm/impls/da/dainterp.c b/src/dm/impls/da/dainterp.c
--- a/src/dm/impls/da/dainterp.c
+++ b/src/dm/impls/da/dainterp.c
@@ -364,7 +364,8 @@
   }

   /* loop over local fine grid nodes setting interpolation for those*/
-  if (!vcoors) {
+  // if (!vcoors) {
+  if (1) {

     for (j=j_start; j<j_start+n_f; j++) {
       for (i=i_start; i<i_start+m_f; i++) {


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.
*
*
*
4) Even through the global coordinates are not directly used, the
interpolant is still valid.
We assume coordinates are interpolated across an element via
 x = \sum N_i ( \xi ) x*_i
where N_i is the usual linear, bilinear, trilinear basis, \xi is the
local coord and x*_i is the nodal value of x
The restriction operator is then
 R_ij = N_i^C( \xi_j^F )
where N_i^C is the basis on the coarse grid and \xi_j^F is the local
coordinate on the fine grid.
In DMDAGetInterpolation() we obviously assemble the transpose of this.
*



> BTW the second version does not work for periodic conditions which is bad!





>  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
>
>  Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
>          Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
>          Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
>          Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
>
>          Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
>          Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
>          Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
>          Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
>
> comes from the trilinear finite element basic functions to define an
> interpolation but since you are evaluating them on a uniform mesh
>
> for (li=0; li<nxi; li++) {
>      xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
>    }
>    for (lj=0; lj<neta; lj++) {
>      eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
>    }
>    for (lk=0; lk<nzeta; lk++) {
>      zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
>    }
>
> 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
>
> arry-smiths-macbook-pro:tests barrysmith$ ./ex36 -dim 2 -nl 1 -cmap 1
> [2 x 2]=>[4 x 4], interp err = 1.4873e+00
>
> 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?
>

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.

We need more tests, something for periodic and something showing why the
(!vcoords) code is worse than the other version.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20110813/503ddaa8/attachment.html>