[petsc-dev] Jed will not be happy

[Dave May]

Hey Barry and Jed,
  Apologies for the late response, I too have been offline for some time.

Sorry about the confusion caused by this statement Barry,
    > 1) We never wanted to support non-uniform grids.
So Jed noted, I actually meant non-uniform grid refinement.

The problem I always encountered related to the case when I used
non-uniform grids.
In such cases, the if (PetscAbsScalar(x) != 0.0)  statement would
typically be entered and an
illegal memory access would occur.

As far as I can tell using my old test code (which was turned into
ex36), by removing the statement

      if (coors) {
        /* only access the next coors point if we know there is one */
        /* note this is dangerous because x may not exactly equal ZERO */
        x = (coors[j][i].x - ccoors[j_c][i_c].x);
        if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[j_c][i_c+1].x -
ccoors[j_c][i_c].x);
        y = (coors[j][i].y - ccoors[j_c][i_c].y);
        if (PetscAbsScalar(y) != 0.0) y = y/(ccoors[j_c+1][i_c].y -
ccoors[j_c][i_c].y);
      }

from what you call the (!NEWVERSION), both implementations currently
in petsc-dev produce identical
errors for the non-periodic case in 2D.

I didn't manage to break the code using
#define NEWVERSION 0
So from my point of view, since both versions pass all the ex36_
tests, they are both good.

I also noticed that the 3D version (DMGetInterpolation_DA_3D_Q1, line
818) appears to be missing the ! symbol.


Cheers,
  Dave



On 14 August 2011 22:23, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>   Ok, I have removed the horrible branch on if (coords) from the interpolation code but left both the old and new versions there. You control the old or new version by defining NEWVERSION at the top of the file. Both versions pass all the ex36_ tests (does that mean that they are both good?).  Once we have some tests where old version breaks I can continue work on them. I have left by default to use the old version since that handles periodic  properly.
>
>  Barry
>
>
>
> On Aug 13, 2011, at 11:23 PM, Jed Brown wrote:
>
>> On Sat, Aug 13, 2011 at 23:03, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> > 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?
>>
>> Dave?
>>
>>
>> > 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 :-)).
>>
>> Some people (us included) solve problems where the coordinates are part of the solution, so having them live in bonified function spaces is useful. The code runs on non-uniform grids (it would be much simpler otherwise), but it assumes uniform refinement. You could also imagine a deformed grid that was not precisely nested within its parent. In that case, you would need point location and a bit more code complexity to handle ghost layer thickness. DA refinement and coarsening is currently all regular, even though the starting meshes may not be. Actually, an exception to this may be if you start with a coarse mesh that is graded for a boundary layer. This case could still be nicely nested, but must need to actually use coordinates, in which case the current code must be incorrect.
>>
>> Dave, was this something you originally intended to support?
>>
>>    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.
>>
>> Are you going to maintain strict nesting such that corner nodes don't move at all and edge/face nodes stay on the original coarse face/edge? I think going the other direction (coarsening from a graded mesh) is more meaningful, but both produce the same situation (described above) where we really should be using 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
>>
>> Haha, good point. So if we fix the interpolants to actually use coordinates in the interior elements, do we silently assume the wrap element is regular or somehow check if the mesh is regular nearby and error if it's not?
>>
>>   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.
>>
>> Do you not have bash-completion?
>
>