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