Question about automatically determining the time step

[Jiangchao3]

hi @krober10nd

I am confused about the method of automatically determining the time step.

In Example_3_ECGC, the time step was set to 0.

When I followed this setting, I found that the time step dt calculated based on the nearshore resolution is very small.

Therefore, the mesh resolution calculated according to the original code of Example 3 is very small, and it runs very slowly.

When I set dt to 2, as in Example_2_NY (Encourage mesh to be stable at a 2 s time step), It turned out to look good.

If using the automatic step size method, is such a small value of dt reasonable to use in the simulation? Or should I specify it myself?

Best wishes
Jiangchao

[krober10nd]

Hey @Jiangchao3 This seems initially at least to be a bug of some sort perhaps when the Courant limiting capabilities were improved #51. I'll take a deeper look into this.

For now, I'd rather you select the timestep yourself.

Keep in mind that after make generation is complete, you should check the maximum Courant number for your timestep.

max(CalcCFL(m,2.0))

and if you still have some problems, you can also bound it after the generation by :

m_new = BoundCr(m,2.0); 

[Jiangchao3]

ok, thanks for your suggestion. I am looking forward to the improvement.

[Jiangchao3]

Hi @krober10nd, I have another question about the low-pass filter (or the parameter of 'fl ')

In the userguide, it's said the local Rossby radius of deformation is used automatically for a low-pass filter cutoff:

And in the section of "edgefx: the mesh size function class", the 'fl' is described as

In the Example_4_PRVI, a negative value of 'fl' is used as the filter equal to Rossby radius divided by 50:

In the article of "Global Storm Tide Modeling with ADCIRC v55: Unstructured
Mesh Design and Performance", the parameter 'FL' is defined ('FL' is defined by the variable parameter alpha), and three value of 1/80, 1/20, 1/5 was applied to the experiments.

I am confused about the definition of the parameter 'fl' or 'FL'：

Is the 'fl' in the userguide and the 'FL' in the article the same parameter?

Compare the equation of FL (L) in the article with the equation of local Rossby radius in the userguide. The parameter alpha is very special.
My understanding：if I set fl=[1,0], the alpha equals to 1, the low-pass filter cutoff is the local Rossby radius. If the value of 'FL' equals to 1/80 (for example as in the article), I should set the fl = -80, this is a value equals to (local Rossby radius/80). Is that correct?

I wanted to use a system default option. However, when I only set a 'slp' parameter, don't set a 'fl' parameter, 'fl' equals to the default value 0 and it will prompt "INFO: Slope filter is off". This means the low-pass filter is off. If I want to use a system default option, should I use the local Rossby radius and set fl=[1,0]?

Please forgive me for writing so much！I am really confused about this parameter.

[krober10nd]

Hey no problem. I'm glad you're inspecting.

There's two filtering methods: a band or high or low pass filter which is activated with a positive value of 'fl'. A band pass filter is activated with an array of two options i.e low and high pass cutoffs.

Whereas the negative option instead uses a local Rossby radius of deformation to determine the low pass cutoff instead of having the user arbitrarily select the lowpass value. And in this negative fl case, the LR Radius is divided by the abs(fl).

Historically we started with just the high and band pass filtering options with the positive fl parameter and that's when we wrote the user guide. Unfortunately we haven't had much time to include some of that new information in the guide. The most update information should be in the help sections in each function.

I would also want to double check this with @WPringle as he performed the experiments for the paper you're mentioning.

[krober10nd]

Hey @Jiangchao3 I updated the user guide to clarify hopefully better the meaning of this fl parameter.

A low pass filter of the bathymetric gradient would be:

`fl`, [1e3, 0]

or

`fl`, 1e3

a band-pass filter of the bathymetric gradient would be:

`fl`, [1e3, 10e3]

Low-pass filtering the bathymetric gradient based on the LRR of Deformation dividing the LRR for example by 80 would be:

`fl`,-80

I hope this helps, by default all filtering is deactivated. I'm updating the user guide to reflect this fact.

PS I'm taking a look at the original question in this thread now.

[Jiangchao3]

Oh, that's great! Thank you very much for this updating. I think this makes it clearer!

I have been learning this software for a while and I have to say that it is really a great job. Thanks again!

[krober10nd]

Hey @Jiangchao3 thanks for the kind words!

I appeared to have fixed the issue you brought up in #116 . The bug had arisen because I built in functionality to not only bound above the maximum Courant number in the mesh but also bound it below by a constant. The latter is important for implicit ALE type numerical schemes that SCHISM users mainly rely on. In those cases, it's actually necessary to ensure a minimum Courant number.

However, the intended usage of the automatic timestep selection was originally for ADCIRC models which often run with a semi-explicit numerical scheme. In this case, I always thought it made more sense to see what's the maximum Courant number you can get away with given your finest nearshore resolution.

[WPringle]

Hi @krober10nd, I have another question about the low-pass filter (or the parameter of 'fl ')

In the userguide, it's said the local Rossby radius of deformation is used automatically for a low-pass filter cutoff:

And in the section of "edgefx: the mesh size function class", the 'fl' is described as

In the Example_4_PRVI, a negative value of 'fl' is used as the filter equal to Rossby radius divided by 50:

In the article of "Global Storm Tide Modeling with ADCIRC v55: Unstructured
Mesh Design and Performance", the parameter 'FL' is defined ('FL' is defined by the variable parameter alpha), and three value of 1/80, 1/20, 1/5 was applied to the experiments.

I am confused about the definition of the parameter 'fl' or 'FL'：

Is the 'fl' in the userguide and the 'FL' in the article the same parameter?

Compare the equation of FL (L) in the article with the equation of local Rossby radius in the userguide. The parameter alpha is very special.
My understanding：if I set fl=[1,0], the alpha equals to 1, the low-pass filter cutoff is the local Rossby radius. If the value of 'FL' equals to 1/80 (for example as in the article), I should set the fl = -80, this is a value equals to (local Rossby radius/80). Is that correct?

I wanted to use a system default option. However, when I only set a 'slp' parameter, don't set a 'fl' parameter, 'fl' equals to the default value 0 and it will prompt "INFO: Slope filter is off". This means the low-pass filter is off. If I want to use a system default option, should I use the local Rossby radius and set fl=[1,0]?

Please forgive me for writing so much！I am really confused about this parameter.

Hi @Jiangchao3, thanks for going through our user guide and the ADCIRC v55 paper in detail. What you have written and what @krober10nd confirmed is correct about the negative fl option where it is a fraction of Rossby radius length. The way I wrote the equation and the number in the article was alpha as a fraction, e.g., 1/80. But in reality that actual fl used as input to oceanmesh2d is fl = -80 for this case.

@krober10nd, I see where you updated the edgefx properties of fl, I also updated the text in Section 3.5 of user guide and added the latest user guide to the static version on ResearchGate.

[Jiangchao3]

@krober10nd Thanks for your work. I have merged the updating code into my branch forked from you. I tested the example 3 again, it displayed the same dt as you did. I will follow your suggestion that sees what's the maximum Courant number after using the automatic timestep for my case.

@WPringle Thanks for your reply. I have clearly understood the meaning of the parameters 'fl' or 'FL' under your and @krober10nd's explanation. Thanks for your discussion about the combination of parameter ‘slp' and 'fl'. I am setting my parameters according to your suggestions.

[Jiangchao3]

Hi @krober10nd, the automatic timestep was calculated as 5.1048s for my case, when I do the calculation of CalcCFL and BoundCr, the dt should be set to 5.1048, otherwise, there will be very strange results.
when I use the script:
CFL = CalcCFL(m);
min(CFL)
max(CFL)
the result was:

However, the userguide describe:

when I use the script:
dt = fh.dt;
CFL = CalcCFL(m,dt);
min(CFL)
max(CFL)
m = BoundCr(m,dt);
Everything is ok!

So I think it is necessary to specify a timestep (using the value of fh.dt) when checking the Cr.

[krober10nd]

Hm that's strange. Are you using the SRTMv2.1 for the DEM in the bigger box?

In any case, bounding the CFL at the mesh sizing stage normally isn't sufficient which is why msh.BoundCr has to exist. That will guarantee the Cr exist is enforced by iteratively adapting mesh resolution.

[Jiangchao3]

I use the topo15_compressed for the DEM. I also tested the Example3_ECGC, when I calculate the CFL by typing:

CFL = CalcCFL(m)

the same strange display again, the max CFL=2631, the min CFL=2.4366

However, the CFL should be < 0.5, isn't it?

I think the problem may not be about the DEM file and maybe here:

The nargin of the above strange case = 1, the code executes the 'else' branch.
If it is CalcCFL(m, fl.dt), the code will execute the first branch and the dt equals to the fl.dt. this case is ok.

So I think it is necessary to specify a timestep (using the value of fh.dt) when calculating the CFL, but why the user guide say: "If you don’t pass it a dt it will output the theoretical maximum possible stable dt in the domain given the mesh
size."

And a further question, in example_3, there are two bbox. the automatic timestep of the outer domain is 6.2979 seconds and the automatic timestep of PostSandyNCEI domain is 1.1134 seconds. How to perform the CalCFL and BoundCr? There are two timestep but only one msh file.

[krober10nd]

I use the topo15_compressed for the DEM. I also tested the Example3_ECGC, when I calculate the CFL by typing:

Yea, I would expect differences given that you're using a different DEM than I did (I used SRTMV2.1).

So I think it is necessary to specify a timestep (using the value of fh.dt) when calculating the CFL, but why the user guide say: "If you don’t pass it a dt it will output the theoretical maximum possible stable dt in the domain given the mesh
size."

Ah ok. Yes, you should specify the stable timestep that was output to the screen when you call msh.CalcCFL (as the second input argument). That second branch in msh.CalcCFL where it calculates the maximum theoretical timestep should be either improved or removed IMO @WPringle . This actually isn't the first time this issue has been brought up.

And a further question, in example_3, there are two bbox. the automatic timestep of the outer domain is 6.2979 seconds and the automatic timestep of PostSandyNCEI domain is 1.1134 seconds. How to perform the CalCFL and BoundCr? There are two timestep but only one msh file.

So then in this case with the high resolution zone around NY, you're numerically restricted by your finest resolution since the mesh is seamless. The multiscale meshing technique is to enable the usage of vastly different horizontal resolution DEMs to be utilized in the same meshing process at the same time.