problem with two-way nesting

[xupeng66]

Hi All,
Here is a problem that I have been struggling with for quite a while. I would like to set up two-way nested simulations to study the small-scale dynamics that are influenced by topographic features such as ridges and seamounts. I have done a few experiments so far, but there is one problem that persists in all of them.

I have attached a few plots from a recent experiment that shows the simulated flow field near a ridge segment at 2100 m depth over the refined grid that is embedded in a coarse grid. There is a map that shows the bathymetry with the full domain and the outline of the nested fine grid. The refinement ratio is 5:1. The two plots that show the flow field are from two parallel runs wherein one-way and two-way nesting are used. Comparing the two plots, the results from the two-way run seems problematic as the currents appear to be disrupted at the boundaries whereas the one-way run shows no such disruption. From the answers I got from the forum on my previous post, I suspect the problem lies in the mismatch between the refined and coarse bathymetries. However, in this current experiment, I have relaxed the refined bathymetry to match the coarse bathymetry near the nesting boundaries. Also, if the problem is caused by bathymetry mismatch, it is interesting that two-way nesting is affected to a much greater degree than one-way nesting.

My understanding of the two-way nesting implemented in ROMS is superficial. But intuition makes me think that maybe the problem lies in the feedback of information from the refined grid to the coarse grid, which is the extra step that is unique to two-way nesting. I am interested to know if anyone has had similar problems before or have any suggestions.

Many thanks!
Guangyu

[wilkin]

I have thought about this issue a lot in the context of our own multiply-nested Mid-Atlantic Bight OOI Pioneer grids.

What I suspect might be happening is deep in the 2-way fine2coarse step. After the fine grid completes its 3 (or 5 in your case) time steps, the fine solution is regridded with a (simple 3x3 or 5x5) smoothing operation performed level by level in s-coordinates. This overwrites the coarse grid within the fine grid domain, and the coarse grid takes a time-step.

If the fine and coarse grids have different bathymetry - anywhere not just near the perimeter - then the smoothed fine solution is inserted at a different vertical depth in the coarse. This could represent a significant change in the potential energy of the water column that triggers an internal gravity wave adjustment in the coarse grid.

If I'm right, we need to do some thinking about applying a vertical regridding step also.

[xupeng66]

Thanks for your reply, John. Your point makes sense to me. I initially thought the problem was limited to the bathymetry difference near the nesting interface. But my attempts to relax the fine-grid bathymetry towards the coarse-grid bathymetry at the interface did not make any difference. So the problem likely occurs when interpolating the results across the entire fine grid onto the coarse grid. Indeed, when refinement grids are used, the entire find-grid becomes the contact zone used to feed information back to the coarse grid so that interpolation is done globally. On the other hand, this problem is gone in one of my experiments wherein the coarse bathymetry is used for both grids so that both grids share the same sigma levels globally.

I agree that adding a vertical regridding step may solve the problem. Yet I cringe at the idea of diving into the source code and implement that change myself. Instead of doing that, do you think adding a sponge layer in the coarse grid near the perimeter of the fine grid to absorb those artificial internal waves originating from the interior of the fine grid can serve as a band-aid solution? I am thinking if the sponge layer is narrow enough, it may preferentially absorb small-scale internal waves and still allow the large-scale barotropic tidal waves to pass through without much distortion?

On the upside, I am pleased to report that nesting (both one-way and two-way) appears to be significantly more efficient in the latest revision of the source code compared with the one released in 2019 .

[stef]

[Edit: just saw you replied while I was posting...the following is intended as reply to wilkin's comment..]

To test this, what do you think about the following:

Mask all coarse grid points in the interior of the fine2coarse contact region, except the ones adjacent to the nest boundary, i.e. the region where the bathymetry of the child was set equal to the parent. This confines the feedback from the fine to the coarse grid to the boundary.

The strategy would be to rely on the fine grid to transport all signals from the interior to the boundary, and only there they are interpolated onto the coarse grid. Effectively it would create an "island" in the coarse grid, except the island is surrounded by a ring-like open nesting boundary condition. The assumption is that what's happening in the interior of the fine grid can be entirely opaque to the coarse grid, as long as the information gets to the boundary.

Not sure if this makes sense physically?

I guess the width of the fine2coarse region would depend on the stencil one uses for the finite difference operators in the advection operator?.

May be worth trying, especially because it may be trivial to implement. Just setting the coarse grid to nan, where the fine grid is (of course leave room for coupling).

I did this in a very simple 2-dimensional setup and it seemed to work. but it was simple, e.g. not much nonlinear advection.

[xupeng66]

I think your suggested test is a good idea. By limiting the fine2coarse contact zone to a narrow band around the nesting boundary, the two-way nesting code may also become more efficient since the costly interpolation from the fine to coarse grid is now done over a much reduced number of grid points. I will try to set up a test with my realistic settings and report the results back to this thread.

[wilkin]

This indeed would be an interesting test. I'll be interested to see what you find out, Guangyu.

I have noted in other posts that in principle it makes no difference what goes on in the portion of the parent grid that is supplanted by the child beyond the coarse2fine contact region. I've mentioned this in the context of river forcing, where rivers in the parent within the child domain have no effect.

But if this is so, then does my speculation about the vertical regridding have any merit? If masking the inner region will do no damage, then how can the vertical regridding?

Perhaps the question is whether the width of the coarse2fine region is enough to span the distance any waves would propagate in the parent? If the parent grid discrete stencil for advection and/or diffusion spans, say, 3 parent grid cells then do we need to match bathymetry precisely over this range, which would be refine_factor times the 3 parent cell stencil width.

But there won't be any economy in this ...

By limiting the fine2coarse contact zone to a narrow band around the nesting boundary, the two-way nesting code may also become more efficient since the costly interpolation from the fine to coarse grid is now done over a much reduced number of grid points

because the fine2coarse gather is done for all points in the child domain, with masking applied afterward.

[stef]

If the parent grid discrete stencil for advection and/or diffusion spans, say, 3 parent grid cells then do we need to match bathymetry precisely over this range, which would be refine_factor times the 3 parent cell stencil width.

Yes, any insight into this would be great!

If I recall correctly, my initial motivation to mask out fine-grid regions in the coarse grid was a drying-wetting problem. I think there was a coarse grid point that suddenly turned wet within the fine-grid region, and caused a blow-up. ROMS' diagnostics don't distinguish between nested cell regions, including the parts that raise the 'exit_flag' due to high velocities etc.

[jcwarner]

I have been working on the nesting for a while now. There are some simple steps that i can take to check on a few things. Is your setup easy to share out?
what ts_advection are you using? time steps?
-j

[xupeng66]

I have been working on the nesting for a while now. There are some simple steps that i can take to check on a few things. Is your setup easy to share out?
what ts_advection are you using? time steps?
-j

The simulation I am doing has a lot of input nc files that contain information on surface forcing, boundary conditions, and bottom heat flux (e.g., bhflux) that are quite large to share. I can send you the other files that you can use to set up a test run with just tidal forcing at the boundary of the coarse grid. I am keen to hear your findings from testing my case and thank you for offering to help. Could you tell me your email address?

[xupeng66]

If the parent grid discrete stencil for advection and/or diffusion spans, say, 3 parent grid cells then do we need to match bathymetry precisely over this range, which would be refine_factor times the 3 parent cell stencil width.

Yes, any insight into this would be great!

If I recall correctly, my initial motivation to mask out fine-grid regions in the coarse grid was a drying-wetting problem. I think there was a coarse grid point that suddenly turned wet within the fine-grid region, and caused a blow-up. ROMS' diagnostics don't distinguish between nested cell regions, including the parts that raise the 'exit_flag' due to high velocities etc.

Stefan, John Wilkin's reply makes me realize that I misunderstood the changes you made in your 2D test case. I was thinking that the masking could be done in the preprocessing step that generates the contact zone input file (e.g., 'roms_ngc.nc') so that the model can be tricked to think that there is an island where the coarse grid is supplanted by the fine grid so that the fine2coarse interpolation ignores that region. Maybe this is not feasible for refinement two-way nesting, where the land mask in the ngc.nc file only blocks out information input from coarse to fine grid but not the other way around. Anyways, I would appreciate if you can send me your revised fine2coarse code so that I can adapt it for a test run with my realistic settings. My email address is guangyux@uw.edu. Thanks in advance.

[jivica]

Gents,
you are running the latest ROMS version with one/two-way nesting constelations so possibly have the answer from the head;
I read long time ago about speedup of nesting algorith in ROMS by excluding some vertical interpolation, masking etc.
My question is what are the steps and how to use those options to speedup calcs?
Can you point me to the thread where it was reported speedup of new vs. old (tried to search but my google is not that good).

Thanks in advance,
I

[stef]

I was thinking that the masking could be done in the preprocessing step that generates the contact zone input file (e.g., 'roms_ngc.nc') so that the model can be tricked to think that there is an island where the coarse grid is supplanted by the fine grid so that the fine2coarse interpolation ignores that region. Maybe this is not feasible for refinement two-way nesting, where the land mask in the ngc.nc file only blocks out information input from coarse to fine grid but not the other way around.

I think the masking could be done even before generating the contact file, it's just a normal land mask.
[EDIT: I just realized there may be a confusion regarding "fine2coarse". In my posts, including this one, I mean the FORTRAN function. Now I remembered the Matlab preprocessing scripts, and saw in [2] that there is a matlab function "fine2coarse.m". Do you mean the Matlab script?]

Hmm, considering John Wilkin's post, maybe my earlier statement:

Mask all coarse grid points in the interior of the fine2coarse contact region, except the ones adjacent to the nest boundary, i.e. the region where the bathymetry of the child was set equal to the parent. This confines the feedback from the fine to the coarse grid to the boundary.

was misleading, in that it suggests that the mask itself does anything. In most cases it may just be cosmetic, in that it blacks out the regions governed by the fine grid in the history files of the coarse grid (exceptions are where these regions cause a blow-up, which I think is not relevant in your case). I guess the main emphasis should have been to make some hypotheses and propose tests for them, i.e.

1) John's (wilkin's) hypothesis could be tested by eliminating the differences in bathymetry (I'm not even sure about that, e.g. what about surface fluctuations, they change the coordinate system as well, don't they? hmm, not sure)
2) the fine2coarse feedback that's actually relevant for coupling the dynamics, happens close to the boundary.
3) To test hypothesis 1 without knowing that hypothesis 2 is valid, one could test hypothesis 1 by eliminating bathy differences in the entire fine domain.
4) To test hypothesis 1 while enforcing coupling close to the boundary, one could use the mask. This way one could leave the fine bathymetry untouched in the majority of the domain, and adjusting it only adjacent to the boundary, as John (wilkin) put it:

If the parent grid discrete stencil for advection and/or diffusion spans, say, 3 parent grid cells then do we need to match bathymetry precisely over this range, which would be refine_factor times the 3 parent cell stencil width.

Obviously 4) would be advantagous to 3) in that the fine model still resolves the fine bathymetry in the majority of the domain. But maybe 3 would be a great test too!


I attached a patch for my 2d toy setup in a previous thread [1], including the changes to nesting.F., but I can't recommend using those in your application, because I think they would break a 3d application (although I have not tried).

[1] viewtopic.php?p=21278#p21278

[2] https://www.myroms.org/wiki/Grid_Proces ... ine2coarse

[xupeng66]

Gents,
you are running the latest ROMS version with one/two-way nesting constelations so possibly have the answer from the head;
I read long time ago about speedup of nesting algorith in ROMS by excluding some vertical interpolation, masking etc.
My question is what are the steps and how to use those options to speedup calcs?
Can you point me to the thread where it was reported speedup of new vs. old (tried to search but my google is not that good).

Thanks in advance,
I

Hi Jivica, through my experiments, both two-way and one-way nesting simulations have become much faster after I updated the source code to the latest revision that was released late last year. I have not done any systematic tests to quantify the speedup and do not know what changes have been made that lead to the speedup, though. Also, I did not use any special CPP options and am not sure if there are any that can be turned on to further speedup nesting simulations. It would be great if some ROMS developer can step in to explain.

[wilkin]

The removal of unnecessary repeated operations in nesting was reported in this ticket:
https://www.myroms.org/projects/src/ticket/861

You can find notes like this by going to the search box in trac https://www.myroms.org/projects/src/search
I simply entered "nesting" and the ticket was the 3rd one there.

[Aljaz]

In my case, when I had a similar boundary problem. The problem was that I didn't have the time step ratios exactly like the grid refinement ratios.

Hope this helps

Aljaz

[xupeng66]

The removal of unnecessary repeated operations in nesting was reported in this ticket:
https://www.myroms.org/projects/src/ticket/861

You can find notes like this by going to the search box in trac https://www.myroms.org/projects/src/search
I simply entered "nesting" and the ticket was the 3rd one there.

John, thanks for bringing up this ticket system. It is a good source of information on revisions to the source code. That ticket on the most recent changes in nesting is also very informative. In particular, it mentions

Currently, the vertical interpolation weights are used in composite grids because their grids are usually not coincident. They are not needed in refinement grids because the donor and receiver grids have the same number of vertical levels and matching bathymetry.

The second sentence seems to suggest that the vertical regridding is not done for refinement grids because those grids are supposed to have matching bathymetry. Having said that, it would be better to add the vertical regridding in refinement nesting since it seems natural to use high-resolution bathymetries for fine grids. I am in the process of setting up an experiment to test the changes suggested by Stefan. Regardless of the test results, it is probably better to add an additional CPP option in ROMS and change the nesting source code so that the vertical regridding can be turned on and off in refinement nesting as needed.

[wilkin]

...it is probably better to add an additional CPP option in ROMS and change the nesting source code so that the vertical regridding can be turned on and off in refinement nesting as needed.

I don't know if anyone has much experience with exercising that code. The nesting has many features that were coded but have yet to be tested in earnest. Figure out how to hack the code to activate the vertical regridding (logical get_Vweight = .TRUE.) with refinement nesting and see what happens. Report back .

[jivica]

Possibly you don't have to change the code, only your contact netcdf file and its flags coincident and mosaic 0/1, as they should be used to IF branch the nesting code.

Ivica

[xupeng66]

Hi All,

I have done a few more tests since my last post on this issue. The following are some interesting new findings.

First of all, as suggested by Stefan, I did a test wherein the fine and coarse girds share the same bathymetry. Surprisingly, the problem with the nesting interface is still there. I have attached a plot showing the coarse-grid flow velocity where there is apparent discontinuity along the boundaries of the fine grid. The refinement ratio is 3:1 in this experiment. This suggests the current problem is not due to the mismatch in bathymetry within the interior of the fine grid, although I believe the lack of vertical regridding in fine2coarse can still cause trouble.

Second, I have compared the fine and coarse grid bathymetries in the corresponding history files to make sure those bathymetries are indeed consistent. Surprisingly, the bathymetries match up everywhere except for the the outermost layer of grid cells, which are presumably the 'ghost points'. What is more puzzling is that the bathymetry in the fine-grid history file differs from the bathymetry in the fine-grid grid file at those 'ghost points'. So it appears that the fine-grid bathymetry has been modified at the ghost points somewhere inside the nesting code, which leads to the mismatch with the coarse-grid bathymetry.

Based on the findings above, I am now suspecting the problem in my two-way nesting experiment may be caused by the bathymetry mismatch at the ghost points on the fine grid. I am curious to know if this is plausible. If so, could anyone tell me why the fine-grid bathymetry is modified from what is in the grid file at those ghost points?

Thanks!

[jcwarner]

the exterior points of the child grid are replaced with information from the contact file during initialization, in get_grid. At those points, data is interpolated from the parent and the nesting needs characteristics of the parent at those locations. I have found some issues with this for the masking, in very shallow water. But we also have a nest running now for the Gulf Stream along the OBX, let me look to see how the currents are at the boundary.

[xupeng66]

the exterior points of the child grid are replaced with information from the contact file during initialization, in get_grid. At those points, data is interpolated from the parent and the nesting needs characteristics of the parent at those locations. I have found some issues with this for the masking, in very shallow water. But we also have a nest running now for the Gulf Stream along the OBX, let me look to see how the currents are at the boundary.

Thanks a lot, John. One thing to add is the issue only occurs in two-way nesting. My one-way results with the same girds and settings look okay. I will look into the bathymetry in my ngc file to see if there is anything suspicious. If so, the problem may be in the Matlab program (contact.m) I used to create the contact file.

[wilkin]

I never saw anything like you are showing for 2-way nesting (3 grids) in our Mid-Atlantic Bight model. Maybe we need a comparison of what our respective codes generate for the contacts point file.

[jcwarner]

we have a 3 nest of USeast - Carolinas - Outer Banks and i am not seeing this behavior of the currents along the boundaries. can you put the grids and the contact files somewhere accessible?

[xupeng66]

I never saw anything like you are showing for 2-way nesting (3 grids) in our Mid-Atlantic Bight model. Maybe we need a comparison of what our respective codes generate for the contacts point file.

we have a 3 nest of USeast - Carolinas - Outer Banks and i am not seeing this behavior of the currents along the boundaries. can you put the grids and the contact files somewhere accessible?

Thanks for checking. I think I have found the source of the error. It appears that the mismatch in bathymetry is due to the different methods used to interpolate coarse-grid bathymetry onto the fine grid. In my fine-grid grid file, the bathymetry is interpolated linearly from the coarse grid. I just saw that in the Matlab program 'contact.m' for creating the contact file, the bathymetry in the contact zone surrounding the fine grid is interpolated from the coarse-grid using the 'spline' method. This is why there is a mismatch in bathymetry at the 'ghost points' of the fine grid. I have reproduced the bathymetry in the fine-grid grid file using spline interpolation. Hopefully, the results will look okay this time. I will post my grid and contact files if the issue persists.

[jcwarner]

either way please let us know. would like to figure this one out.

[jivica]

Really good detection of any problems with nesting is to compute free surface gradient, for both parent and child and then plot them together.
I did that for my 2-way case using simple python script, here is snip of the most important part:

Code: Select all

....
def roms_fs_grad(zeta, pm, pn):
    dphi_y, dphi_x = np.gradient(zeta, axis=0)
    return dphi_x*pm, dphi_y*pn

def make_plot_grad_two(rp, pvar,rc, cvar, title, fileout):
    fig, ax = make_map(projection = ccrs.PlateCarree())
    ax.set_extent(bbox, crs = ccrs.PlateCarree())
    ax.add_geometries(shp.geometries(), ccrs.PlateCarree(), facecolor='gray')
    cf = ax.contourf(rp.lon_rho, rp.lat_rho, pvar, levels = clevs, cmap = cmap, 
                     transform = ccrs.PlateCarree(), extend = 'max')
    c2 = ax.contourf(rc.lon_rho, rc.lat_rho, cvar, levels = clevs, cmap = cmap, 
                     transform = ccrs.PlateCarree(), extend = 'max')
    cb = plt.colorbar(cf, orientation = 'vertical', ticks = zrange, extend = 'max', 
                      extendrect = False, shrink=0.6, pad = 0.02, label = 'cm/km')
    plt.title(title)
    if fileout != None:
        plt.savefig(fileout, dpi = 144, transpearent = True)
    plt.close()

p = xarray.open_mfdataset('archive_qck_000*.nc', concat_dim='ocean_time')
c = xarray.open_mfdataset('archive_qck_nest_000*.nc', concat_dim='ocean_time')

px,py = roms_fs_grad(p.zeta[160,:,:],p.pm, p.pn)
cx,cy = roms_fs_grad(c.zeta[160,:,:],c.pm, c.pn)
grad_p = np.abs(px + 1j*py)*1e5
grad_c = np.abs(cx + 1j*cy)*1e5

# plot them for time index 160
make_plot_grad_two(p,grad_p[160,:,:],c,grad_c[160,:,:],
                   'gradient of free surface', 'gradient_zeta_parent_child.png')

My child grid (500m) doesn't match parent bathymetry (2500m) inside its whole domain (it is better-high res in the middle part),
but is exactly the same within 10 grid cells around its bry, and then further inside I am linearly tappering to fine-child high resoution bathymetry (within next 40 cells).

Recently, I saw some strange behaviour at the shallow boundary (southern) using ultimate test of free surface gradient.
I am expecting to have sharper solution for internal tides within nested domain, but think there are some unphysical parallel stripes close to bry.

Will check about interpolation in contact.m and true values of bathy in the parent and child history output.

Regards,
Ivica

[xupeng66]

either way please let us know. would like to figure this one out.

Here are my new test results. The first three files shows the flow velocity at 2100 m depth and the last three plots show gradients of sea-surface elevation as suggested by Ivica. I have also attached my grid files and contact file.

First of all, the velocity appears to have improved after correcting the mismatch of bathymetry in the contact file in a sense that there is now less disruption near the nesting interface. However, the interface still stands out in the plot of coarse-grid velocity, so the propagation from fine to coarse grid is still not perfect. Again, this issue only shows up in two-way nesting. In addition, as suggested by Ivica, I computed the gradients of sea-surface. Surprisingly, the results also show horizontal stripes that look suspicious. Also, those stripes are present in both two-way and one-way results, which are very similar to each other. I am thinking maybe those stripes in sea-surface gradients are waves that bounce back from the nesting interface because they are not resolvable on the coarse grid? I would think the averaging in fine2coarse will smooth out waves on fine scales but probably the smoothing is not 100%.

I never saw anything like you are showing for 2-way nesting (3 grids) in our Mid-Atlantic Bight model. Maybe we need a comparison of what our respective codes generate for the contacts point file.

John, I have attached the Matlab code (contact.m) I used to generate the contact file. It seems that I changed the code at some point because SVN shows that it is different from the original. Not sure if my changes are the culprit of all the problems I am having.

[jivica]

I just checked effect of the different options in the contact.m for interpolation (akima, linear, cubic) and results are the same -- parallel stipe pattern still exist in grad zeta (differences are super small for different interp options). I am doing diffs after 7 days in simulation after transients are gone. My feeling is we have problem with grids time stepping mismatch and to eliminate this I will run both grids with the same time stepping. This is waste of CPU but could eliminate time stepping possiblity. I even run one-way nesting and the stipes are still there.

Another thing is puzzling me, it is the difference between bathymetry in the grid file and the roms output (qck for example).
Simple NCO command: ncbo -y sub -vh grid.nc qck_0001.nc -o diff_depth.nc gives different bathymetry along the open boundary rim which can go +- 5m in my case (250m depth).
This should be exactly the same as I am not updating depth i.e. depth evolution in time via sediments, and this is static field in ROMS.

Cheers,
Ivica

[xupeng66]

Hi Ivica,

I think the difference in bathymetry between the grid file and the output file you see is the same as what I saw in my previous test. See the quote below from one of my earlier posts.

Second, I have compared the fine and coarse grid bathymetries in the corresponding history files to make sure those bathymetries are indeed consistent. Surprisingly, the bathymetries match up everywhere except for the the outermost layer of grid cells, which are presumably the 'ghost points'. What is more puzzling is that the bathymetry in the fine-grid history file differs from the bathymetry in the fine-grid grid file at those 'ghost points'. So it appears that the fine-grid bathymetry has been modified at the ghost points somewhere inside the nesting code, which leads to the mismatch with the coarse-grid bathymetry.

To eliminate that difference, you need to make sure the bathymetry in your fine-grid grid file is interpolated from the coarse-grid bathymetry using the same method (e.g., spline) as in contact.m. Having said that, I agree that those horizontal stripes we see in our respective results are caused by an issue other than the mismatch in bathymetry. Maybe it is a problem of time stepping as you pointed out. I am curious to hear what you find out from your test with the same time stepping between coarse and fine grids.

[takashi]

I have also encountered a similar problem of the artificial u/v fluctuation at the boundaries between the coarse and refine grids.

In my understanding, the boundary values of ubar/vbar for the refine grid are stored in REFINED(cr)%DU_avg2(1,m,tnew)/REFINED(cr)%DV_avg2(1,m,tnew) through the get_persisted2d subroutine and are set through the routine in the u2dbc_im.F/u2dbc_im.F.
On the refine grid boundary, all the ubar/vbar values of the refine grid cells inside a cell of the coarser grid are set at the same value as the coarser grid ubar/vbar value through the get_persisted2d subroutine.
It means that the same velocity values' blocks (five cells' block in the case of refine_factor = 5) are aligned on the boundary, and the discontinuous velocity exists between the blocks.

I understand that this treatment is for conserving mass-flux through the refine grid boundaries, but it seems that this treatment makes step-like changes in u/v velocity on the boundaries and is sometimes induced small artificial eddies near the discontinuous area of velocity.

To resolve this problem, I modified the Lines 5002-5025 of get_persisted2d subroutine in nesting.F with new CPP flag of REFINE_LINIEAR_DUV_AVE2_BOUNDS as follows:

Code: Select all

        IF (((Istr.le.Idg).and.(Idg.le.Iend)).and.                      &
     &      ((Jstr.le.Jdg).and.(Jdg.le.Jend))) THEN
          IF (on_boundary(m+m_add).gt.0) THEN
            IF ((on_boundary(m+m_add).eq.1).or.                         &
     &          (on_boundary(m+m_add).eq.3)) THEN       ! western and
!!! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TN: add
# if defined REFINE_LINIEAR_DUV_AVE2_BOUNDS
              Ac(1,m) = contact(cr)%Lweight(1,m)*Ad(Idg,Jdg) +          &
     &                  contact(cr)%Lweight(2,m)*Ad(Ip1,Jdg) +          &
     &                  contact(cr)%Lweight(3,m)*Ad(Ip1,Jp1) +          &
     &                  contact(cr)%Lweight(4,m)*Ad(Idg,Jp1)
              Ac(2,m) = Ac(1,m)  ! not used?
              Ac(3,m) = Ac(1,m)  ! not used? 
              Ac(4,m) = Ac(1,m)  ! not used? 
# else
!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< TN: add
              j_add=INT(REAL(Jrg-1,r8)*Rscale)          ! eastern edges
              j=J_bottom(rg)+j_add
              Ac(1,m)=Ad(Idg,j)
              Ac(2,m)=Ad(Idg,j)
              Ac(3,m)=Ad(Idg,j)
              Ac(4,m)=Ad(Idg,j)
!!! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TN: add
# endif
!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< TN: add
            ELSE IF ((on_boundary(m+m_add).eq.2).or.                    &
     &               (on_boundary(m+m_add).eq.4)) THEN  ! southern and
!!! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TN: add
# if defined REFINE_LINIEAR_DUV_AVE2_BOUNDS
              Ac(1,m) = contact(cr)%Lweight(1,m)*Ad(Idg,Jdg) +          &
     &                  contact(cr)%Lweight(2,m)*Ad(Ip1,Jdg) +          &
     &                  contact(cr)%Lweight(3,m)*Ad(Ip1,Jp1) +          &
     &                  contact(cr)%Lweight(4,m)*Ad(Idg,Jp1)
              Ac(2,m) = Ac(1,m)  ! not used?
              Ac(3,m) = Ac(1,m)  ! not used? 
              Ac(4,m) = Ac(1,m)  ! not used? 
# else
!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< TN: add
              i_add=INT(REAL(Irg-1,r8)*Rscale)          ! northern edges
              i=I_left(rg)+i_add
              Ac(1,m)=Ad(i,Jdg)
              Ac(2,m)=Ad(i,Jdg)
              Ac(3,m)=Ad(i,Jdg)
              Ac(4,m)=Ad(i,Jdg)
!!! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TN: add
# endif
!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< TN: add
            END IF
          ELSE
            Ac(1,m)=Ad(Idg,Jdg)               ! contact point is not at
            Ac(2,m)=Ad(Ip1,Jdg)               ! physical boundary, just
            Ac(3,m)=Ad(Ip1,Jp1)               ! set values for spatial
            Ac(4,m)=Ad(Idg,Jp1)               ! interpolation (not used)
          END IF
        END IF

This edit seems not a smart solution and may induce a new problem related to mass-flux conservation..., but, in my case, but the problem improved after applying the above modification.
I hope this helps.

[xupeng66]

Hello,

I have done a couple more tests and found that those artificial stripes in sea-surface gradients only appear when tidal forcing is applied. Those stripes disappeared when I turned off tides. My understanding is that tidal forcing is applied at the open boundary of the coarse grid in a nested simulation with refinement grids. The tides are then introduced into the fine grid through the contact points. I do not know why adding tides in my simulation has led to the artifacts in sea-surface gradients. The following are the CPP options I used. I am hoping someone can let me know if any of those options are inappropriate.

# define NESTING
# define ONE_WAY
# define SOLVE3D
# define SPHERICAL
# define CURVGRID
# define UV_ADV
# define UV_COR
# define UV_LOGDRAG
# define UV_SMAGORINSKY
# define MIX_GEO_UV
# define DJ_GRADPS
# define TS_SMAGORINSKY
# define MIX_GEO_TS
# define SALINITY
# define NONLIN_EOS
# define BULK_FLUXES
# define ATM_PRESS
# define ANA_SSFLUX
# define ANA_BTFLUX
# define ANA_BSFLUX
# define RADIATION_2D
# define STATIONS
# define AVERAGES

# define GLS_MIXING
# if defined GLS_MIXING
# define CANUTO_A
# define N2S2_HORAVG
# define RI_SPLINES
# endif

# define SSH_TIDES
# define UV_TIDES
# define ADD_FSOBC
# define ADD_M2OBC
# define RAMP_TIDES

[wilkin]

Boundary tides are going to input a lot more barotropic energy than without. I don't think it's directly relevant to the nesting issue you are seeing, just making it more apparent because there is so much more high frequency zeta variability.

Clutching at straws a bit, but ... try without SMAGORINSKY.

And try #define MIX_S_UV instead of #define MIX_GEO_UV

Also, what is your NDTFAST? Are you resolving the split-explicit coupling well.

[xupeng66]

Hi John and Ivica,

As Ivica suspected

My feeling is we have problem with grids time stepping mismatch and to eliminate this I will run both grids with the same time stepping

,
I reduced my coarse-grid time step from 30 sec to 10 sec to match the one used for the fine grid. The refinement ratio between the coarse and fine grids is 3 in my experiment. Surprisingly, those horizontal stripes in sea-surface are gone on the fine grid as shown in the attached plots. The first one is the previous result obtained with coarse grid dt = 30 sec and fine-grid dt=10 sec, the second one is the new result obtained with coarse-grid dt = fine-grid dt=10 sec.

So Ivica, it seems that your suspicion of a time-stepping problem is on point. Could you please explain why using different time steps in refinement nesting could cause problems like this? My understanding is that the nesting in ROMS is designed to have the coarse and fine grids run at different time steps that match the grids' respective resolutions.

John, to answer your question, the 'NDTFAST' in my simulation is 30, which should satisfy the CFL criterion based on my computation. The deepest spot in my full domain is at h = 2800 m depth. So the phase speed of barotropic waves is sqrt(g*h) = 166 m/s. The grid spacing on my fine grid is 80 m. So the maximum dT in the barotropic mode is dt = 80 m /166 m/s = 0.5 sec. The actual barotropic dt used for the fine grid is 10 sec/30 = 0.33 sec. So the CFL is about 0.66. Maybe that is this still too high?

[stef]

There is approximately 5km between the stripes and the barotropic wave speed is 166m/s, so there is a period of 30s, which is in sync with the baroclinic time step? So assuming the signal is generated at the nesting boundary and propagating into the fine domain, maybe check the effect of the cpp option TIME_INTERP_FLUX. You mentioned in a previous post that waves "bounce back from the nesting interface". Is that really how they propagate? I

But note that Hernan mentions in an old post this should not be activated:

viewtopic.php?p=15221#p15221

I have not done any 3d experiments but the temporal interpolation of the fluxes also made me wonder in the 2d case. If I remember correctly, u/vbar and zeta are temporally interpolated in the coarse2fine region, but not at the actual boundary condition in u2dbc. I have not looked into it.

[xupeng66]

There is approximately 5km between the stripes and the barotropic wave speed is 166m/s, so there is a period of 30s, which is in sync with the baroclinic time step? So assuming the signal is generated at the nesting boundary and propagating into the fine domain, maybe check the effect of the cpp option TIME_INTERP_FLUX. You mentioned in a previous post that waves "bounce back from the nesting interface". Is that really how they propagate? I

But note that Hernan mentions in an old post this should not be activated:

viewtopic.php?p=15221#p15221

I have not done any 3d experiments but the temporal interpolation of the fluxes also made me wonder in the 2d case. If I remember correctly, u/vbar and zeta are temporally interpolated in the coarse2fine region, but not at the actual boundary condition in u2dbc. I have not looked into it.

Hi Stef, I think your observation is quite revealing. The baroclinic time step for the coarse grid happens to be 30 sec. The one for the fine grid is 10 sec. I presume the data on the coarse grid is fed into the fine grid through the contact points at each coarse-grid time step or 30 sec. So it seems that some perturbation is introduced into the fine grid each time it acquires information from the coarse grid. I now suspect my contact point file is till problematic. I just saw that one of my previous simulations does not have the problem of eta stripes, although that one has much lower resolutions. I am now backtracking to see which changes I have made may have caused the problem. So one step forward, two steps back

[xupeng66]

I just took another look at my old results and it turned out that the horizontal stripes in eta gradients were also present. I missed that because I did not think to check sea-surface gradients back then.

That said, I also reread Arango's post on setting up nesting grids: viewtopic.php?t=3663. One thing I have missed before is that it is preferable to build nested grids using the fine-to-coarse strategy, which involves building a large, fine intermediate grid first and extract from it the final coarse and refined grids to use. I have been doing it the opposite way, which is to build a large, coarse grid and then extract from it the fine grid using coarse2fine.m. The reason given by Arango for taking the first approach is that

This will give us the grid data needed for all contact points

. I am struggling to understand this. So what grid data would be missing if the coarse-to-fine strategy is used like what I have been doing?

Anyways, I am now back to the drawing board and will rebuild my grids using the fine-to-coarse strategy and restart from there.

[arango]

Yes, that's the way to build complicated grids. One thing that you may consider is to avoid sudden bathymetry changes like seamount and ridges across the contact areas between coarse and fine grids. I have done telescoping grids in the Middle Atlantic Bight, which include the Gulf Stream, with tidal forcing and I didn't notice the issues that you are having. I turn on the writing of 2D and 3D vorticity always to check the behavior across the nested grid's contact points. It will give you more information about any leakage, circulation, and physics in the contact area.

Building the grids from fine-to-course using an intermediate grid is the way to go! You will have the target bathymetry in the finer grid and can have smooth transitions in the bathymetry gradients when coarsening. You will have all the needed grid variables in the contact points and masked areas: h, f, pm, pn, and angle with a smooth transition. That is, you don't have to recompute such values in an inconsistent way, and have the appropriate land/sea mask for the finer grid. The coarser fields are much easier to compute. Also, it is much easy to impose volume conservation.

[xupeng66]

Thanks for your reply, Hernan. I agree that having the nesting interface go through steep topography is not a good idea. However, in my case, it is hard to avoid abrupt changes in bathymetry since my focus site is on the crest of a linear ridge so the interface is bound to cross the ridge at some point. I now see the point of taking the fine2coarse approach. I will rebuild my grids and see if that solves the problem.

[xupeng66]

Yes, that's the way to build complicated grids. One thing that you may consider is to avoid sudden bathymetry changes like seamount and ridges across the contact areas between coarse and fine grids. I have done telescoping grids in the Middle Atlantic Bight, which include the Gulf Stream, with tidal forcing and I didn't notice the issues that you are having. I turn on the writing of 2D and 3D vorticity always to check the behavior across the nested grid's contact points. It will give you more information about any leakage, circulation, and physics in the contact area.

Building the grids from fine-to-course using an intermediate grid is the way to go! You will have the target bathymetry in the finer grid and can have smooth transitions in the bathymetry gradients when coarsening. You will have all the needed grid variables in the contact points and masked areas: h, f, pm, pn, and angle with a smooth transition. That is, you don't have to recompute such values in an inconsistent way, and have the appropriate land/sea mask for the finer grid. The coarser fields are much easier to compute. Also, it is much easy to impose volume conservation.

Hi Hernan, while I am working on rebuilding my nested grids using the recommended fine-to-coarse approach, Matlab reports two errors in contact.m and fine2coarse.m respectively, both of which are caused by the use of spherical coordinates. The one in contact.m happens on line 329 when assigning the structure 'R' to S.refined(cr=2). The problem occurs because when the grids are spherical, the structure 'R' still has Cartesian coordinate variables such as 'x_rho' and 'y_rho' but those variables are missing in S.refined. So Matlab reports a problem of 'subscripted assignment between dissimilar structures'. The error in fine2coarse.m is similar, when the grids are spherical, the Cartesian coordinate variables are all empty but fine2coarse is still trying to use them as if those variables are filled. Both of the errors are easy to fix. I changed contact.m a while ago to make it work for the nested simulations posted in this thread. However, those error messages make me wonder maybe those Matlab nesting scripts have not been fully tested on spherical grids. In addition, maybe it is better to use Cartesian coordinates for refinement nesting because maybe it is easier to conserve areas between the coarse and fine grids in Cartesian coordinates. Please let me know if this is the case.

[arango]

What software did you use to generate your original GRID NetCDF? Many users really miss the fact that although the application has a spherical grid, it is highly recommended to have the associated Cartesian coordinates to compute the metrics pm and pn since ROMS computational grid is in curvilinear coordinates. Recall that users have the option of having selected geographical map projections when generating the grid. This kind of information is very useful in generating nesting grids.

I just updated the ROMS Matlab repository. Check trac ticket for more information. Check function roms_metrics.m and study it. Notice that the function grid_metrics.m is deeprecated and will be removed in the future

[xupeng66]

What software did you use to generate your original GRID NetCDF? Many users really miss the fact that although the application has a spherical grid, it is highly recommended to have the associated Cartesian coordinates to compute the metrics pm and pn since ROMS computational grid is in curvilinear coordinates. Recall that users have the option of having selected geographical map projections when generating the grid. This kind of information is very useful in generating nesting grids.

I just updated the ROMS Matlab repository. Check trac ticket for more information. Check function roms_metrics.m and study it. Notice that the function grid_metrics.m is deeprecated and will be removed in the future

Thanks again for the quick response, Hernan. I have just updated my local copy of the 'roms_tool' repository. As for the software I used to build my grids, I have been using 'GridBuilder', which I find to be quite convenient with a user friendly GUI interface. However, I now see that the lack of Cartesian coordinates in the spherical grid files created by GridBuilder could cause problems. I will add the Cartesian coordinates back in and that should solve the problems I have with contact.m and fine2coarse.m. I will keep you posted on my progress.

[xupeng66]

I just updated the ROMS Matlab repository. Check trac ticket for more information. Check function roms_metrics.m and study it. Notice that the function grid_metrics.m is deeprecated and will be removed in the future

Hi Hernan, Here is one quick update. I have just tried roms_metrics.m. The x/y coordinates calculated from lon/lat did not look right. I found a potential error on line 288:
y_rho(i,j+1) = x_rho(i,j) + dy(i,j+1); I think the 'x_rho' on the right hand side should be 'y_rho'. I have made the change and the x/y coordinates look okay now.

[arango]

Yes, thank you!

[xupeng66]

Building the grids from fine-to-course using an intermediate grid is the way to go! You will have the target bathymetry in the finer grid and can have smooth transitions in the bathymetry gradients when coarsening. You will have all the needed grid variables in the contact points and masked areas: h, f, pm, pn, and angle with a smooth transition. That is, you don't have to recompute such values in an inconsistent way, and have the appropriate land/sea mask for the finer grid. The coarser fields are much easier to compute. Also, it is much easy to impose volume conservation.

Hi Hernan, I went ahead to rebuild my nested grids following your fine-to-coarse approach. However, I did not go very far before running into further issues, so some more help will be appreciated.

Here is what I did. First, I built a large intermediate grid with a high resolution that is I would like to use for the fine grid. I then extracted a coarse grid from the intermediate grid using fine2coarse.m with Gfactor=3, which lowers the resolution three times. I did not specify the indices (e.g., imin, imax, jmin, jmax) so that fine2coarse extracted the largest coarse grid possible from the intermediate grid. I then extracted a small, fine grid from the intermediate grid at the same resolution using grid_extract.m.

The question came up in the next step, which is to plug the coarse and fine grids into contact.m to create the ngc file. Because those grids are extracted from the intermediate grid, the global attributes in those grid files such as 'refine_factor', 'parent_Imin', 'parent_Jmin', etc. are all relative to the intermediate grid, which is not used in contact.m. This seems to confuse the program and cause it to stop with an error of 'inconsistent refined grid extraction'.

I figure I must have not followed the fine-to-coarse approach correctly. So could you please explain with a little more details about the procedures for extracting a coarse grid with a fine grid nested inside from an intermediate grid and how to generate the ngc file using contact.m from there?

[arango]

Yes, I recall issues with the global attributes. In such cases, I manipulate the global attributes using my matlab/netcdf functions that are distributed in the ROMS Matlab repository:

Code: Select all

>> Imin=nc_getatt('my_file.nc', 'parent_Imin');
>> Imax=nc_getatt('my_file.nc', 'parent_Imax');
>> Jmin=nc_getatt('my_file.nc', 'parent_Jmin');
>> Jmax=nc_getatt('my_file.nc', 'parent_Jmax');
>> refine_factor=nc_getatt('my_file.nc', 'refine_factor');

>> nc_attadd('my_file.nc', 'parent_extraction', int32([Imin Imax Jmin Jmax refine_factor]));

>> nc_attdel('my_file.nc', 'refine_factor')
>> nc_attdel('my_file.nc', 'parent_Imin');
>> nc_attdel('my_file.nc', 'parent_Imax');
>> nc_attdel('my_file.nc', 'parent_Jmin');
>> nc_attdel('my_file.nc', 'parent_Jmax');

Since it is important information to have, I read old attribute values, create a new global attribute with such values, and delete the problematic attributes. It is very simple. Users of ROMS need to learn how to create and modify the NetCDF files. We provide lots of generic Matlab functions to do so. If the user is allergic to Matlab, they can use other software that is available elsewhere.

[xupeng66]

Yes, I recall issues with the global attributes. In such cases, I manipulate the global attributes using my matlab/netcdf functions that are distributed in the ROMS Matlab repository:

Code: Select all

>> Imin=nc_getatt('my_file.nc', 'parent_Imin');
>> Imax=nc_getatt('my_file.nc', 'parent_Imax');
>> Jmin=nc_getatt('my_file.nc', 'parent_Jmin');
>> Jmax=nc_getatt('my_file.nc', 'parent_Jmax');
>> refine_factor=nc_getatt('my_file.nc', 'refine_factor');

>> nc_attadd('my_file.nc', 'parent_extraction', int32([Imin Imax Jmin Jmax refine_factor]));

>> nc_attdel('my_file.nc', 'refine_factor')
>> nc_attdel('my_file.nc', 'parent_Imin');
>> nc_attdel('my_file.nc', 'parent_Imax');
>> nc_attdel('my_file.nc', 'parent_Jmin');
>> nc_attdel('my_file.nc', 'parent_Jmax');

Since it is important information to have, I read old attribute values, create a new global attribute with such values, and delete the problematic attributes. It is very simple. Users of ROMS need to learn how to create and modify the NetCDF files. We provide lots of generic Matlab functions to do so. If the user is allergic to Matlab, they can use other software that is available elsewhere.

Hi Hernan and all,

Thanks a lot for all the instructions. I managed to rebuild my grids using the fine-to-coarse approach and finished a test run using one-way nesting. The refinement ratio between the coarse and fine grids is 3. I have attached an animation showing hourly sampled sea-surface gradients over a 2-day span since the start of the simulation. Unfortunately, those much dreaded horizontal strips are still there despite the change in the domain size of the fine grid and switching to the fine-to-coarse approach. Note that those stripes seem to come and go like a surface wave originating from the southern nesting interface. As Stefan pointed out previously, the wavelength suggests a period of 30 sec, which is the coarse-grid time step. So it seems that a perturbation in sea surface is introduced at the southern interface each time the fine grid reads outputs from the coarse grid. I do not know if this is a time stepping issue as suspected by Ivica but my previous test seemed to suggest those stripes are gone after using the same time step for both fine and coarse grid.

Having said that, I still suspect there is something wrong in my grids and the ngc file created using 'contact.m'. Hernan, in order to make the grids extracted from the intermediate grid work for contact.m, I manipulated the global attributes in the coarse and fine grids as follows. I added the attribute 'refine_factor=0' to the coarse grid file and change the value of that attribute from '1' to '3' in the fine grid file to match the refinement ratio between those grids. The reason for those changes is that 'contact.m' seems to determine whether the grids are nested in the refinement configuration based on the value of 'refine_factor'. As mentioned in the instructions in the Matlab code, 'refine_factor' needs to be 0 for the coarse grid. I have attached my grid files and ngc file. I would much appreciate if you can take a look and see if there is something wrong. I am hoping a fresh pair of eyes can spot some error that has long eluded me.

[jivica]

As you know time stepping in ROMS is a elaborated and complex beast and exchanging info between grids for fast varying variables is a nighmare.
If you are modelling non-tidal, slowly evolving dynamics you are most likely OK.
However, if you are in the space of high frequency motions then, I think, there is no simple solution (except unefficient method of using the same time step).
We did long time ago simple nesting 1-way -- downscalling experiment in tidaly (internal tides as well) driven region and there was a problem of introducing spurious oscillations, similar like aliasing problem (Janekovic, I., & Powell, B. (2012). Analysis of imposing tidal dynamics to nested numerical models. Continental Shelf Research, 34, 30-40. https://doi.org/10.1016/j.csr.2011.11.017).

Not sure, but it could help to use quadratic time interpolation when exchainging info between grids (instead of linear if it is used).

Regards,
Ivica

[stef]

[edit: deleted]

[stef]

But if there is a problem with the grid topology, how come the experiment with the same timestepping does run fine. Hmm, it's surprising, isn't it?

[xupeng66]

As you know time stepping in ROMS is a elaborated and complex beast and exchanging info between grids for fast varying variables is a nighmare.
If you are modelling non-tidal, slowly evolving dynamics you are most likely OK.
However, if you are in the space of high frequency motions then, I think, there is no simple solution (except unefficient method of using the same time step).
We did long time ago simple nesting 1-way -- downscalling experiment in tidaly (internal tides as well) driven region and there was a problem of introducing spurious oscillations, similar like aliasing problem (Janekovic, I., & Powell, B. (2012). Analysis of imposing tidal dynamics to nested numerical models. Continental Shelf Research, 34, 30-40. https://doi.org/10.1016/j.csr.2011.11.017).

Not sure, but it could help to use quadratic time interpolation when exchainging info between grids (instead of linear if it is used).

Regards,
Ivica

But if there is a problem with the grid topology, how come the experiment with the same timestepping does run fine. Hmm, it's surprising, isn't it?

Hi Stef and Ivica,

After scrutinizing the results from the simulation conducted with the same dt between the coarse and fine grids, I now think those stripes are present in both cases. I have attached two snapshots of sea-surface gradients from two parallel runs respectively. It appears that those stripes are also present when the coarse-grid dt is reduced to match the fine-grid dt. The difference is that those stripes are now more closely spaced because the coarse-grid outputs are fed into the fine-grid at shorter time intervals (10 sec vs 30 sec before). So I still suspect the problem is due to the some hidden error in my contact file.

Ivia, I am curious to know whether you have solved the problem of artificial stripes near the nesting interface in your case. Do you think your problem and mine share any similarities?

[xupeng66]

Hi All,

I have done a few more tests over the past couple of weeks. In these new tests, I have removed bottom topography, open boundary, and surface forcing so that the simulations are forced with tides alone. Although this did not solve the problem of those artificial stripes in sea-surface gradients, the much simplified settings did make those artifacts stand out more clearly. It is evident that those stripes originate from the nesting interface and propagate from there into the small domain. What is more surprising is that those artifacts also appear near the open boundaries of the coarse grid. This makes me wonder if those artifacts are a general problem with open boundary conditions rather than an exclusive problem with nesting.

I also did a sanity test on the Lake Jersey test case using the original input files I downloaded from the 'test case' repository. I have attached two animations showing sea-surface elevation and its gradients over the mother grid (i.e., grid a). There is only one nested grid (grid c) in this experiment. The refinement ratio is 5:1. Surprisingly, artificial stripes in sea-surface gradients showed up again both within the daughter grid and near the open boundary of the mother grid. This finding seems to suggest the problem is likely not caused by some error in the input files I created myself. Also, since the Lake Jersey test case only has wind forcing, those artifacts are not exclusively related to tides either.

I am keen to hear your thoughts on this.

Guangyu

[xupeng66]

Hi All,

Here are some more test results from a realistic experiment. The model grids include a coarse grid and a fine grid embedded within that is centered at the Axial Seamount in the Northeast Pacific. The top of the seamount is at 1500 m depth. The grids were built using the fine-to-coarse approach recommended by Hernan with a refinement ratio of 5:1. I have done two experiments using one-way and two-way nesting respectively. The model flows are driven by tidal, and sub-tidal forcing at the coarse-grid lateral boundaries plus surface wind/pressure. The same settings were used for both experiments except for the nesting scheme (one-way vs two-way). Special care was taken to make sure the bathymetry between the two grids matches so that the grids share the same sigma levels.

The plots below show 2D velocities and the 3D velocities sampled at 1000 m depth. The two-way results are on the left and the one-way results are on the right. For those one-way results, I averaged the fine-grid solutions based on the 5:1 refinement ratio and used the averages to replace the coarse-grid solutions at collocated grid points. It is noticeable that the two-way results show apparent discontinuities near the nesting interface. The disruption in flow mostly occurs at the locations where the bathymetry changes significantly across the nesting interface. In comparison, the one-way results are much smoother, in which the nesting interface appears almost transparent, which is a sign that the interface is behaving properly. Such differences between one-way and two-way experiments are persistent in all the tests I have done that feature realistic bathymetry as shown in my previous posts.

I know that the two-way nesting differs from its one-way counterpart mainly by having an additional fine2coarse step, in which the fine-grid solutions are averaged and substituted for the coarse-grid solutions at collocated grid points. It seems to me that the averaging currently done in the fine2coarse subroutine could become problematic when the bathymetry changes significantly over the five fine-grid cells that occupies a same coarse-grid cell. In that case, the simple averaging done in fine2coarse for 2D variables, which does not take into account changes in bottom depth, may not yield results that match the coarse-grid solutions. This will result in discontinuities in mass fluxes across the nesting interface when the fine grid 2D variables are passed to the coarse grid.

If my supposition above is correct, I am wondering if that is the reason why it is recommended to place the nesting interface where the bathymetry is smooth. I am also wondering if one can modify the averaging done in fine2coarse to make it more robust in the presence of significant bathymetry across the fine-grid cells that are averaged over. Any thoughts?

[wilkin]

Take a look at ...

Debreu, L., Marchesiello, P., Penven, P. and Cambon, G., 2012. Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean Modelling, 49, pp.1-21, http://dx.doi.org/10.1016/j.ocemod.2012.03.003

Let us know if you have an suggestions on how to improve the 2-way communication.

[jcwarner]

it might be on the call to set_depth before the call to fine2coarse of the r3dvars in nesting. just a long shot but if i could get a copy of this setup, i can dig into this.
-j

[arango]

As you may have noticed, the nesting routines in nesting.F module were designed in a generic way to facilitate expansion and different schemes. The weighted area average is turned off in fine2coarse (AreaAvg=.FALSE.). One could try volume-weighted averaging (not coded). Right now, we use bilinear interpolation. We could implement some type of conservative interpolation to minimize this problem, but it is more involved and require some thinking and testing.

One of the reasons that I recommend putting the contact points away from sharp bathymetry gradients is that it will be tricky to have each coarse grid cell with the same cumulative volume as the finer subgrid cells inside for ratios 3:1, 5:1, and so on. This is controlled by your input bathymetry in all nested grids. Some bathymetry volume conservation smoothing is required. Did you try a 3:1 ratio instead? Maybe your 5:1 ratio is problematic for this application.

Also, maybe some conservative fluxes term(s) need to be coded and computed in the shallow-water kernel for each fast time step. It will complicate matters! The AGRIF implementation considers such treatment.

Anyway, why do you need nesting in this application? It will not be too expensive to run this application in a single grid at the desired target resolution. This application is not that big even for my 8 CPUs MacPro laptop. Nowadays, computers are very powerful and disk storage is very cheap. I always advise users to remove the complexity in their applications and try the easiest approaches first. Then, add the complexity.

[stef]

If you aren't using identical bathymetry in the contact regions, the following may be interesting too:

viewtopic.php?f=14&t=5438&p=21072#p21072

I haven't updated my ROMS code in a while, so I'm not sure if this is still relevant.

I thinkt the problem was that the barotropic volume flux through the coarse cell is simply divided by the refinement ratio, and then prescribed at the fine sub-cells, regardless of how deep they are.

[xupeng66]

Hi All, sorry for taking so long to write back.

To answer Stef first, I read your old post and see how that could be problematic. But I do not think that is the cause because that issue seems to affect the input of information from coarse to fine grid and thus the artifacts seen in my results would be present in both one-way and two-way nesting. Please let me know if you agree.

Arango, I am very surprised to hear that you could run such simulations on your laptop with only 8 CPUs. The grid size of the coarse grid is 1000 by 1000 at 250 m intervals with 64 sigma layers. The time step size is 30 sec. A 30-day single-grid simulation takes me about a day to finish using 72 CPUs. Increasing the resolution 5 fold will need a grid size of 5000 by 5000 and time steps 5 times smaller, which will take much longer to finish and too much storage space on the computer I am using. Maybe there is some weird setup issue that slows down the simulation at my end substantially?

Wilkin, I read the paper you recommended and saw that the 'fine2coarse' step, which is called 'update' by those authors, is quite different from what is implemented in the current version of ROMS. There are two major differences that I noticed. First, the update of coarse-grid solution only occurs in a 'feedback' interface on the coarse that lies near the boundary the fine-grid rather than the entire interior of the fine grid. Second, as shown in their Figure 7, taking the average of fine-grid solutions does not properly damp the small-scale features that are unresolved on the coarse-grid. The procedures described in that paper also includes 'flux correction' as mentioned by Hernan. I presume the nesting scheme described in that paper is used in the other implementation of ROMS that is called CROCO. I have downloaded that model and will play with it when I have more fee time.

Lastly, my temporary solution is to stick to one-way nesting to generate results for my project. But I am keen on continuing this dialogue with you guys. I think the online nesting (both one-way and two-way) capability of ROMS is a powerful feature that seems to be underused since the feature was added several years ago. It would be great if someone who is more experienced with nested simulations can write a general post summarizing the pros and cons of the different nesting schemes (online vs offline, one-way vs two-way) available.

[stef]

But I do not think that is the cause because that issue seems to affect the input of information from coarse to fine grid and thus the artifacts seen in my results would be present in both one-way and two-way nesting. Please let me know if you agree.

Yes, agreed! I have yet to try out one-way nesting..

[xupeng66]

Hi All,

I am writing this post to let everyone know that the problem with two-way nesting that prompted me to start this thread has been resolved by John Warner after revising the treatment of fine-to-coarse information exchange in nesting.F. The details of those changes are described in this ticket: https://www.myroms.org/projects/src/ticket/887. The attached is a plot of flow velocity from my latest two-way nested simulation. The artifacts near the nesting interface seen in previous results are gone after using the revised nesting.F. I would like to take this opportunity to thank John and all of you for being so supportive.

Cheers!
Guangyu