clawpack / pyclaw

PyClaw is a Python-based interface to the algorithms of Clawpack and SharpClaw. It also contains the PetClaw package, which adds parallelism through PETSc.
http://www.clawpack.com/pyclaw
BSD 3-Clause "New" or "Revised" License
157 stars 98 forks source link

Use fortran ordering for q array when reading ascii files. #697

Closed ketch closed 1 year ago

ketch commented 1 year ago

This fixes an error reported here: https://groups.google.com/g/claw-users/c/lcDROYpRMxY/m/sDW_gINPCAAJ?utm_medium=email&utm_source=footer

The issue is that when reading ASCII files, the solution.state.q array is created with C ordering instead of Fortran ordering.

rjleveque commented 1 year ago

I don't see how this changes anything, since passing a shape to zeros returns an array of that shape, and all its elements are 0 regardless of what order you fill it with zeros, e.g.

In [36]: zeros((2,3), order='F')
Out[36]: 
array([[0., 0., 0.],
       [0., 0., 0.]])

In [37]: zeros((2,3), order='C')
Out[37]: 
array([[0., 0., 0.],
       [0., 0., 0.]])

I would think it's when filling the array with solution values that the order matters, but that's done in nested loops.

rjleveque commented 1 year ago

Another thought: This routine is used by visclaw in reading ascii output from the Fortran codes in classic, amrclaw, and geoclaw. So if you change the order that arrays are filled in this routine, I fear all of those will break. Maybe the problem is the order that PyClaw writes the solution, rather than this routine?

ketch commented 1 year ago

The order you see when accessing the array doesn't change, but the layout in memory does. If you load a solution with the current master branch, PyClaw will tell you it's not valid, because the F_contiguous flag is false. With this patch, that flag is true.

ketch commented 1 year ago

The error reported is the result of this check failing: https://github.com/clawpack/pyclaw/blob/07ebb8019c9a85d6a142c9215717203128f970a1/src/pyclaw/state.py#L186

ketch commented 1 year ago

This is what it does:

Screen Shot 2022-12-10 at 1 13 29 PM

It sounds like you are worried that this change will transpose the arrays, but that's not the case (even when they are passed to C or Fortran routines). It will make things a little more efficient when arrays are passed to Fortran, since they won't need to be copied.

mandli commented 1 year ago

Makes sense to me and looks good.

rjleveque commented 1 year ago

Thanks for the clarification. I guess I don't see how this fixes a bug in pyclaw, but I'm also not understanding what the issue was from https://groups.google.com/g/claw-users/c/lcDROYpRMxY/m/sDW_gINPCAAJ?utm_medium=email&utm_source=footer, since the problem reported there seems to be entirely related to num_ghost missing.

But if it makes things more efficient in pyclaw, it's fine to do this. It has no effect on the way we use this in the Fortran/visclaw versions.

I wasn't worried about this change transposing arrays, but I thought the issue was that arrays were transposed relative to what pyclaw needed, and that re-writing the loops filling the arrays with values would break something else.

ketch commented 1 year ago

Sorry, I explained very poorly. The user fixed the missing num_ghost issue himself, after which he got a message Initial solution is not valid. I was able to replicate this and traced it to the fact that the code checks for the F_CONTIGUOUS flag and finds it is not True.

Indeed, if you don't ever pass the resulting array to a Fortran routine, then this won't have any effect (even in terms of efficiency).