Different behavior with single and multiple grids at the same refinement level

[carlosmunozmoncayo]

Hello, I am working with a shallow water-type system on GeoClaw. The solver works fine for most tests, but there are some scenarios where I encounter some unexpected behavior when using multiple grids, even if they all have the same refinement level. For instance, if I consider this problem (all the variables initialized as 0 except for the depth),
I get this with a single grid:

and this with multiple grids:

Since I am using a single refinement level in both cases, shouldn't the solutions be exactly the same? I get this behavior with qad.f turned on and off. I suspect this is related to the way the aux array values are passed across patch boundaries (maybe it is related to #653 ?).
My Riemann problems depend on the aux array which, for this scenario, does not depend on time:

aux array single grid	aux array multiple grids

The system I am solving has 5 conserved quantities (the first three ones are the same as the SWEs), and has larger eigenvalues, which I am accounting for in the Riemann solver speeds and getmaxspeed.f90. The routines I have more extensively modified are:

rpn2.f90
rpt2.f90
setaux.f90
src2.f90

I have also made some small modifications to:

amr2.f90
filval.f90
tick.f
valout.f90

But, to my understanding, nothing that would interfere with the way aux arrays are handled at patch boundaries. Could this be a GeoClaw or AMRClaw bug?

[mandli]

It definitely looks like something wrong in the patch boundaries. I doubt it's related to #653 as that only happens if the aux arrays are time dependent. I have seen problems like this before, but I think it's just a common symptom rather than indicative of a deeper cause. You mentioned though that you were using a different set of wave speeds to determine time steps, could it be that this is not working in the ghost cells?

First thing that would be useful to see though is the velocities (momentum) that you are getting as these types of problems tend to start in the momentum terms first.

[mjberger]

Filpatch handles the 4th and 5th equations differently than the first 3, and assumes they are the Boussinesq variables. (For example it doesn't limit them and check for positivity, etc). There might be other places like this. I assume this doesn't happen with the regular shallow water with 3 variables, correct? So this gives a place to look. I can help if you need it.

…

On Jul 7, 2025, at 7:30 AM, Carlos Muñoz ***@***.***> wrote: carlosmunozmoncayo created an issue (clawpack/geoclaw#659) <#659> Hello, I am working with a shallow water-type system on GeoClaw. The solver works fine for most tests, but there are some scenarios where I encounter some unexpected behavior when using multiple grids, even if they all have the same refinement level. For instance, if I consider this problem (all the variables initialized as 0 except for the depth), I get this with a single grid: image.png (view on web) <https://github.com/user-attachments/assets/4429e6c3-6da8-414f-af8d-69d6d6392134> image.png (view on web) <https://github.com/user-attachments/assets/c7e0bb8f-f124-4294-935a-5f7f672b1c6e> and this with multiple grids: image.png (view on web) <https://github.com/user-attachments/assets/c86ab40e-13ae-4398-be67-86d0b18212da> image.png (view on web) <https://github.com/user-attachments/assets/89e90239-8226-4ae5-9103-2e16c46f580c> image.png (view on web) <https://github.com/user-attachments/assets/5a489efa-f139-4833-b8ba-6296881f2168> Since I am using a single refinement level in both cases, shouldn't the solutions be exactly the same? I get this behavior with qad.f turned on and off. I suspect this is related to the way the aux array values are passed across patch boundaries (maybe it is related to #653 <#653> ?). My Riemann problems depend on the aux array which, for this scenario, does not depend on time: aux array single grid aux array multiple grids image.png (view on web) <https://github.com/user-attachments/assets/4ad17446-38ea-4f94-89cb-c6189a39291e> image.png (view on web) <https://github.com/user-attachments/assets/a276994c-0a7f-4604-96b6-7a0f7fa31735> The system I am solving has 5 conserved quantities (the first three ones are the same as the SWEs), and has larger eigenvalues, which I am accounting for in the Riemann solver speeds and getmaxspeed.f90. The routines I have more extensively modified are: rpn2.f90 rpt2.f90 setaux.f90 src2.f90 I have also made some small modifications to: amr2.f90 filval.f90 tick.f valout.f90 But, to my understanding, nothing that would interfere with the way aux arrays are handled at patch boundaries. Could this be a GeoClaw or AMRClaw bug? — Reply to this email directly, view it on GitHub <#659>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAGUGC3R6ZWLSCBRLUNZBTD3HJK5ZAVCNFSM6AAAAACA57JCK2VHI2DSMVQWIX3LMV43ASLTON2WKOZTGIYDQNRTHA2TENA>. You are receiving this because you are subscribed to this thread.

[mjberger]

to debug it I would use the coarsest grid you can, and check after 1 time step that the aux arrays are exactly the same in both cases.

…

On Jul 7, 2025, at 7:30 AM, Carlos Muñoz ***@***.***> wrote: carlosmunozmoncayo created an issue (clawpack/geoclaw#659) <#659> Hello, I am working with a shallow water-type system on GeoClaw. The solver works fine for most tests, but there are some scenarios where I encounter some unexpected behavior when using multiple grids, even if they all have the same refinement level. For instance, if I consider this problem (all the variables initialized as 0 except for the depth), I get this with a single grid: image.png (view on web) <https://github.com/user-attachments/assets/4429e6c3-6da8-414f-af8d-69d6d6392134> image.png (view on web) <https://github.com/user-attachments/assets/c7e0bb8f-f124-4294-935a-5f7f672b1c6e> and this with multiple grids: image.png (view on web) <https://github.com/user-attachments/assets/c86ab40e-13ae-4398-be67-86d0b18212da> image.png (view on web) <https://github.com/user-attachments/assets/89e90239-8226-4ae5-9103-2e16c46f580c> image.png (view on web) <https://github.com/user-attachments/assets/5a489efa-f139-4833-b8ba-6296881f2168> Since I am using a single refinement level in both cases, shouldn't the solutions be exactly the same? I get this behavior with qad.f turned on and off. I suspect this is related to the way the aux array values are passed across patch boundaries (maybe it is related to #653 <#653> ?). My Riemann problems depend on the aux array which, for this scenario, does not depend on time: aux array single grid aux array multiple grids image.png (view on web) <https://github.com/user-attachments/assets/4ad17446-38ea-4f94-89cb-c6189a39291e> image.png (view on web) <https://github.com/user-attachments/assets/a276994c-0a7f-4604-96b6-7a0f7fa31735> The system I am solving has 5 conserved quantities (the first three ones are the same as the SWEs), and has larger eigenvalues, which I am accounting for in the Riemann solver speeds and getmaxspeed.f90. The routines I have more extensively modified are: rpn2.f90 rpt2.f90 setaux.f90 src2.f90 I have also made some small modifications to: amr2.f90 filval.f90 tick.f valout.f90 But, to my understanding, nothing that would interfere with the way aux arrays are handled at patch boundaries. Could this be a GeoClaw or AMRClaw bug? — Reply to this email directly, view it on GitHub <#659>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAGUGC3R6ZWLSCBRLUNZBTD3HJK5ZAVCNFSM6AAAAACA57JCK2VHI2DSMVQWIX3LMV43ASLTON2WKOZTGIYDQNRTHA2TENA>. You are receiving this because you are subscribed to this thread.

[carlosmunozmoncayo]

Thank you for your replies @mandli and @mjberger!
@mandli This is what I get for the velocities:

It seems like the problem starts once the waves interact with patch boundaries close to sharp changes in the aux array. I don't think the main cause is the Riemann solution speeds because I get the same behavior with a fixed time step size and no second-order corrections (the solver uses f-waves, so the speeds would only affect the direction of propagation).

@mjberger I am using filpatch from the shallow water equations. I understand that, in that one, all the variables but the depth are treated equally. In this scenario, only intfil, trimbd, and bc2amr are called from filpatch, since all the patches are at level 1. I will try with a coarser setup. The idea would be to access the ghost cell values of a grid from, e.g., b4step and compare them with the boundary values of the neighboring grid, right?

[mandli]

To look at the bathymetry on the entire patch, b4step is a good place to do that. I have a sneaking suspicion that the representation of the bathymetry is probably not different, but could to check. The fact that this is happening with only one level removes a lot of possibilities though, which is helpful.

A couple of questions for you:

• What sort of bathy jumps are required for this to occur? Are they required? Do large jumps need to be aligned with patch boundaries?
• Does this occur when in lat-long and cartesian coordinates?

There are a few other things you could try so that we can trace down the problem:

• Try different options for the transverse waves, from turning off limiting to turning off the transverse solver all together.
• Try turning off the qad routine (if you look at it there is a return right at the top that can be uncommented). At the very least this will remove some complexity from the problem.
• You could set max1d to something larger and see how that impacts the simulation (it is a parameter in the rundata.amrdata object).
• Try simplifying the example as much as possible, i.e. only one jump in bathymetry and a plane wave aligned with the jump. This can try and remove and 2D effects.

[mjberger]

i suggest restarting from just before the first problem. Print out all state and aux variables in the patches and one full grid (on the coarsest level you can), including the locations. Then you sort the location since the integer indices won't be the same, but the x and y values will. Then you can diff then.

…

On Jul 13, 2025, at 6:28 PM, Carlos Muñoz ***@***.***> wrote: carlosmunozmoncayo left a comment (clawpack/geoclaw#659) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3067342436&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=Re_ywEC2SdLGQDdwWPg7ZAQTJuD5HHs4CVEBtEZsI_0&e=> Thank you for your answers @mandli <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_mandli&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=_9DS9qirXJeBFEUDambmtcMMb5L711PbFcPe-ub-uAM&e=> and @mjberger <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_mjberger&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=-9keVvYki2RYUSQnQPkE3sGfDdSySeVIXl5m4B-6WpU&e=>! @mandli <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_mandli&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=_9DS9qirXJeBFEUDambmtcMMb5L711PbFcPe-ub-uAM&e=> This is what I get for the velocities: image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_d4cc7d2e-2Ddb03-2D4a80-2D8982-2Dcbab361783f9&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=ebA3aVJHnj-sqVYdxYHWsd2qni0Fh0jOLbdasaLrj4E&e=> image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_19d8fbad-2D7514-2D4f1e-2Da2c1-2D93c1b3f15ca0&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=FILROl9zQu3SvuWlPHQ7ksbC7DOfT0manDUW5SOcqXQ&e=> image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_7a9dcbdb-2D47d9-2D4cf0-2D812c-2D1f3493736a85&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=uVR3_k7jAPfh5i06TKpgtdnYcWoWg2wSmL2JR7wXeZs&e=> image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_27cbb635-2D91c1-2D45e9-2D827e-2D0609d2428e39&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=mPzPCN9kEjcQY7ph4bVv36mkfaQWR4MAK4fjWcRjmok&e=> image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_36193958-2D2121-2D49b7-2Db5e0-2D303530b2c864&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=rosits6oUq3CEi6V4ghsVKu3Goj8TTOIj6h6R8wDnt0&e=> image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_83cce5a1-2D2b65-2D41a0-2D942d-2D343d8dfb4ab6&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=pYGefiC3cXDNZHtkegK9ZmdsaZppD6LAiXBMlejWNR8&e=> It seems like the problem starts once the waves interact with patch boundaries close to sharp changes in the aux array. I don't think the main cause is the Riemann solution speeds because I get the same behavior with a fixed time step size and no second-order corrections (the solver uses f-waves, so the speeds would only affect the direction of propagation). @mjberger <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_mjberger&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=-9keVvYki2RYUSQnQPkE3sGfDdSySeVIXl5m4B-6WpU&e=> I am using filpatch from the shallow water equations. I understand that, in that one, all the variables but the depth are treated equally. In this scenario, only intfil, trimbd, and bc2amr are called from filpatch, since all the patches are at level 1. I will try with a coarser setup. The idea would be to access the ghost cell values of a grid from, e.g., b4step and compare them with the boundary values of the neighboring grid, right? — Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3067342436&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=Re_ywEC2SdLGQDdwWPg7ZAQTJuD5HHs4CVEBtEZsI_0&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AAGUGCZVMRUUXTTJ3DRBQYD3ILMRXAVCNFSM6AAAAACA57JCK2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTANRXGM2DENBTGY&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=kYPogjHZvjV7pGQhSlBOXHGFR5GvpmH8D-vScCG8GKb-39rrQfp0cq5gPQ0bt27Y&s=67bHQWOPLpwQUHcJExoKFXxfSKdaQvJC4KOsdwR1SCo&e=>. You are receiving this because you were mentioned.

[carlosmunozmoncayo]

This only occurs when I have large bathymetry jumps close to patch boundaries. However, my additional aux variable is a piecewise affine function of the bathymetry, so large bathy jumps also lead to large jumps in this variable. I suspect this aux variable is the source of the discrepancy. So far, I've only observed this in lat-long coordinates; I’ll try to reproduce it in Cartesian coordinates.

These are really insightful suggestions, I'll try them out and report back. Thank you!

[carlosmunozmoncayo]

Hello, I've tried your suggestions, and I think I found the source of the discrepancy: some components of the aux array at ghost cells for patches seem not to be properly defined.

Considering the same setup, but with only 4 patches, I am taking a single step and writing to a file the q and aux arrays close to the boundaries of each patch. The values at ghost cells match those of the interior boundaries of neighboring grids for the first components of the aux array (bathymetry, capacity function, and length ratio), but they differ for the additional component of the aux array I need for this system (I call it indicator in the figure below):

In my setaux routine, I am setting this component by looping only through the interior cells. It looks like this:

neighborhood_idx(:,1) = [-1,-1,-1,0 ,0, 1,1,1]
neighborhood_idx(:,2) = [-1, 0, 1,-1,1,-1,0,1]
do jj=1,my
    do ii=1,mx
        h_at_rest = sea_level - aux(1,ii,jj)
        !Go around the neighborhood of this cell
            do kk=1,8
                iidx = ii + neighborhood_idx(kk,1)
                jidx = jj + neighborhood_idx(kk,2)
                h_at_rest = min(h_at_rest,sea_level - aux(1,iidx,jidx))
            end do
        aux(csq_index,ii,jj) = f(h_at_rest)
    enddo
enddo

I think this issue would be fixed if I looped all the way through the ghost cells here, like it's done for the capacity function. However, I cannot do this since this aux component requires information from the neighboring cells. So that component of the aux array stays unassigned (it seems like it defaults to zero even though I am not explicitly doing it in setaux).

Even if I was able to define the aux array at ghost cells from setaux, shouldn't those values get overridden by values at the boundary of neighboring patches?

[mandli]

Which version of the code are you working with?

[carlosmunozmoncayo]

I'm working with version 5.10.0. Could this be affected by changes in recent releases?

[mandli]

Possibly, there were some changes we made due to weirdness around how ghost cells and aux were filled. I had made a comment above about that only being relevant to time dependent cases, but I am not so sure. Would be good to check if you see anything different.

[mjberger]

The new stuff only affects setting the ghosts cells for time dependent aux arrays, and does not loop over grids to interpolate ghost cells a second time. Are these components time-dependent? Filpatch should be setting ghost cells for these components using interior cells from other patches. I would verify that you are really seeing this indicator set in filpatch. You can make 2 grids instead of 4 to make it easier to print out all the ghost cells (you can do these by setting the number of cells not equal and max1d to in-between. If you can't find it I have debugger I can help you with if you can sent me something easy to install. Marsha

…

On Sep 16, 2025, at 2:35 PM, Kyle Mandli ***@***.***> wrote: mandli left a comment (clawpack/geoclaw#659) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3299906312&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=bDxilm9-gwV2uBK_VICPJwyBdEJWcgM_CUrDe_dCnXct6tNrMXdczcTIC8boZf44&s=sWRKu76CCUALh1K9jpAjGnYwC-zuh6NrJYYuLO7wZ1A&e=> Possibly, there were some changes we made due to weirdness around how ghost cells and aux were filled. I had made a comment above about that only being relevant to time dependent cases, but I am not so sure. Would be good to check if you see anything different. — Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3299906312&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=bDxilm9-gwV2uBK_VICPJwyBdEJWcgM_CUrDe_dCnXct6tNrMXdczcTIC8boZf44&s=sWRKu76CCUALh1K9jpAjGnYwC-zuh6NrJYYuLO7wZ1A&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AAGUGC4KBF6JTRM7M2JZHCL3TBKAFAVCNFSM6AAAAACA57JCK2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTEOJZHEYDMMZRGI&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=bDxilm9-gwV2uBK_VICPJwyBdEJWcgM_CUrDe_dCnXct6tNrMXdczcTIC8boZf44&s=372I0rlEmetQNXA1SrnYQP5R-NbZgCgfDSnE_4SFicY&e=>. You are receiving this because you were mentioned.

[carlosmunozmoncayo]

The additional components of the aux array are not time-dependent. In this setup, I am not modifying the aux array in any other place than in the loop above from setaux, and I am not calling setaux from b4step. In any case, I will try with the changes from the latest release to see if there's any difference.

I think this issue does not occur only for the additional aux array components, but it shows up because, for these, I am not looping over the ghost cells. For instance, if I modify the loop in which the capacity function is set in setaux from this

do jj = 1 - mbc, my + mbc
    do ii = 1 - mbc, mx + mbc
        ...
    end do
end do

to this

do jj = 1 , my
    do ii = 1 , mx
        ...
    end do
end do

I run into the same issue. Instead of this (looping over the ghost cells in setaux):

I get this (setting only interior cells in setaux):

Is this expected behavior? I am printing all of these values from the b4step function.

[mjberger]

I"m looking at your loops more closely I have them backwards in my mind. In the first one you are setting the ghost cells, which gives you the good interior ghost cells on the outside, but not the second one. The interior "ghost cells" should be set in filpatch from the other interior grids that intersect. That is where you should look.

…

On Sep 18, 2025, at 6:32 PM, Carlos Muñoz ***@***.***> wrote: carlosmunozmoncayo left a comment (clawpack/geoclaw#659) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3309904178&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=Evh8chm3aEvk7ZePXRcTYqB3adUu9w2EsOdChB99qHHqubq5skG1-ddwQ7jKa2zu&s=kjsWQ99V5y06NWRfBWRU9zXxuDew7U9xzmL2PB9hH24&e=> The additional components of the aux array are not time-dependent. In this setup, I am not modifying the aux array in any other place than in the loop above from setaux, and I am not calling setaux from b4step. In any case, I will try with the changes from the latest release to see if there's any difference. I think this issue does not occur only for the additional aux array components, but it shows up because, for these, I am not looping over the ghost cells. For instance, if I modify the loop in which the capacity function is set in setaux from this do jj = 1 - mbc, my + mbc do ii = 1 - mbc, mx + mbc ... end do end do to this do jj = 1 , my do ii = 1 , mx ... end do end do I run into the same issue. Instead of this (looping over the ghost cells in setaux): image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_e68db3f2-2D148c-2D4c0e-2Db74a-2Dd1dd9bfe6124&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=Evh8chm3aEvk7ZePXRcTYqB3adUu9w2EsOdChB99qHHqubq5skG1-ddwQ7jKa2zu&s=Q0-dA3pK4KE5MU2moltqkWEJ_5e-b-We-1vO6QGeZcg&e=> I get this (setting only interior cells in setaux): image.png (view on web) <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_user-2Dattachments_assets_f50cd0c0-2Dfae2-2D4848-2Dab34-2D85c61812973d&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=Evh8chm3aEvk7ZePXRcTYqB3adUu9w2EsOdChB99qHHqubq5skG1-ddwQ7jKa2zu&s=92g-xmEkB6z6_CMQALZtpl9fM1nDwwvZHMvgbuMaqXM&e=> Is this expected behavior? I am printing all of these values from the b4step function. — Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_clawpack_geoclaw_issues_659-23issuecomment-2D3309904178&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=Evh8chm3aEvk7ZePXRcTYqB3adUu9w2EsOdChB99qHHqubq5skG1-ddwQ7jKa2zu&s=kjsWQ99V5y06NWRfBWRU9zXxuDew7U9xzmL2PB9hH24&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AAGUGC4KYES3WZ736SFZ3433TMXGPAVCNFSM6AAAAACA57JCK2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTGMBZHEYDIMJXHA&d=DwMCaQ&c=slrrB7dE8n7gBJbeO0g-IQ&r=4PEQeJ-eTbHJv3d0pLWzLQ&m=Evh8chm3aEvk7ZePXRcTYqB3adUu9w2EsOdChB99qHHqubq5skG1-ddwQ7jKa2zu&s=RUGEV76uuZAuLWov_lzFNPp3ykUWxk-ncmhe_beRpnU&e=>. You are receiving this because you were mentioned.

[carlosmunozmoncayo]

It seems like setting

do_aux_copy = .true.

in bound.f90 fixes the issue:

Now I am getting the same behavior both single and multiple grids. I guess that do_aux_copy = .false. is the default since, when possible, setting the aux array at ghost cells from setaux is faster than copying it.
Thank you very much for your help!