[petsc-dev] Jed will not be happy

[Barry Smith]

On Aug 13, 2011, at 8:30 PM, Jed Brown wrote:

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

    It is possible the old code was producing incorrect results; but it is actually suppose to be doing pretty much the same as the new code in its design and logic (though not as clear) so I am not sure why it would generate 
different answers except for bugs. Any examples showing a difference?

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

    I'm totally confused. So it is only suppose to do uniform grids with uniform grid refinement? Ok but that is what the old code did. So is the new code just a "fix" for the old code (plus broken periodic) and should it have just replaced the old code? Why does the test code create a nonuniform coordinate vector then if it isn't used. Are all the tests done as tests on interpolating coordinates? Not on interpolating some function on the grid (yes I know coordinates are a particular case of functions on the grid, but no one I care about :-)).


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

    That wasn't my plan. If I do a DARefine (say by 2) but then set in the coordinates on the finer mesh some scale coordinates is the DAGetInterpolation suppose to ignore the scaled coordinates?  Looks like it.

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

   I don't understand this AT ALL. Currently the code completely ignores all coordinate information for non-periodic boundary conditions (assuming I guess uniform everything?) but you get all hot and bothered by the fact that you don't have coordinate values to do the interpolation for periodic case so you have to skip it. This seems completely contradictory to me?  I can use the periodic interpolation generated (with the old code) to do multigrid on a periodic domain just fine.  The "new" code could likely be very easily fixed to do the periodic case by using the same business with the unit element on the "extra" elements along the edge (like is done in the old code). Is all your rejection of the "periodic case" because you cannot interpolate "coordinates" in the periodic region? Big hairy deal

   
> 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_*?

   No I tried make runex36 per the standard naming convention and it said no such test. I didn't bother looking further than that figuring no one would break the standard naming convention.


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

   For me the DMGetInterpolation() is for getting an interpolation that maps between functions that live on two grids (for example for multigrid); you guys seem obsessed with using it to interpolate coordinates? Is that the root of the confusion? 

  Sounds like I should go through and "merge" the old code and the new code; currently I think they should actually generate the same answers for "uniform meshes" which seem to be the only thing anyone is supporting so I wonder why they do not.  Any examples that show a difference between the old and new so I can debug the old?

   Barry


> 
> We need more tests, something for periodic and something showing why the (!vcoords) code is worse than the other version.