Barotropic model Blowing-up with high courant and dt = 1.

Message

asantana · #1 Unread post by **asantana** » Fri Aug 23, 2019 10:31 am

Hello,

This is my first post here.
So, I would like to say thank you for accepting me and for all ROMS developers for the tool and support.
I am modeling my first problem with ROMS.
It is a 2D model (barotropic only) with a resolution of 10m.

[Problem]
The model is blowing-up in step 9 with high barotropic velocities.

Code: Select all

      STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
            C => (i,j)       Cu            Cv           C Max       Max Speed

      0 18707 00:00:00  0.000000E+00  0.000000E+00  0.000000E+00  5.811092E+07
             (000,000)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
      DEF_HIS   - creating history file, Grid 01: /mnt/waverider/data/BRIFS_Tests/Outputs/ROMS/20190811_op_wrf3.6/roms_BRIFS_child_20190811_op_his_0001.nc
      WRT_HIS   - wrote history  fields (Index=1) into time record = 0000001
      1 18707 00:00:01  0.000000E+00  0.000000E+00  0.000000E+00  5.811092E+07
             (000,000)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
      2 18707 00:00:02  1.084924E-14  0.000000E+00  1.084924E-14  5.811092E+07
             (001,001)  8.706092E-08  9.115097E-08  1.782119E-07  2.079068E-06
      3 18707 00:00:03  1.969264E-13  1.809879E-13  3.779143E-13  5.811092E+07
             (002,002)  3.749301E-07  3.928720E-07  7.678021E-07  4.664422E-06
      4 18707 00:00:04  5.893866E-13  7.733326E-13  1.362719E-12  5.811092E+07
             (004,004)  5.819446E-07  6.100008E-07  1.191945E-06  8.985014E-06
      5 18707 00:00:05  7.659377E-12  9.666165E-12  1.732554E-11  5.811092E+07
             (001,003)  1.448769E-06  2.626366E-06  4.075135E-06  4.763931E-05
      6 18707 00:00:06  4.201891E-10  4.859953E-10  9.061843E-10  5.811092E+07
             (002,004)  9.736066E-06  1.712538E-05  2.686145E-05  7.233986E-04
      7 18707 00:00:07  5.394890E-08  4.766217E-08  1.016111E-07  5.811092E+07
             (001,003)  5.420863E-04  4.928241E-05  5.913688E-04  1.775107E-02
      8 18707 00:00:08  5.679739E-05  3.656752E-05  9.336491E-05  5.811093E+07
             (001,002)  1.399511E-02  3.388924E-03  1.738403E-02  5.952867E-01
      9 18707 00:00:09  7.606368E-02  4.657086E-02  1.226345E-01  5.811088E+07
             (001,004)  5.894658E-01  6.879632E-02  6.582621E-01  2.162523E+01

Blowing-up: Saving latest model state into  RESTART file

Looking at the log files it is clear that I have a problem with the Courant numbers:

Code: Select all

 Minimum barotropic Courant Number =  1.19183864E+00
 Maximum barotropic Courant Number =  3.12685140E+00
 Maximum Coriolis   Courant Number =  9.37217488E-05

I could reduce my dt to control it ([1], [2], and [3]), but the problem is that I already have a dt = 1.

[Question]
How can I handle the Courant number in this case?
Or any ideas to identify where is my problem?

[More Information]
To identify the problem I did some tests.
First, to simplify my problem, I turn off the ATM_PRESS flag.
I have open boundaries in the South and West.
Also, I use Chapman implicit for free-surface and Flather for ubar and vbar.
I checked my boundary conditions and the velocities are okay:

Code: Select all

    GET_NGFLD   - free-surface western boundary condition,   t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -1.94486187E-04 Max = -1.64477533E-04)
    GET_NGFLD   - free-surface southern boundary condition,  t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -2.58851418E-04 Max = -1.64477533E-04)
    GET_NGFLD   - 2D u-momentum western boundary condition,  t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -8.54891041E-05 Max = -8.04100491E-05)
    GET_NGFLD   - 2D v-momentum western boundary condition,  t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -9.36471333E-05 Max = -9.36471333E-05)
    GET_NGFLD   - 2D u-momentum southern boundary condition, t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -8.54891041E-05 Max = -7.47929153E-05)
    GET_NGFLD   - 2D v-momentum southern boundary condition, t = 18707 00:02:00
                   (Rec=0000002, Index=1, File: roms_BRIFS_child_bry_20190811_op.nc)
                   (Tmin=      18707.0000 Tmax=      18709.0000)
                   (Min = -9.36471333E-05 Max = -8.32194564E-05)

This is my header with cppflags definition:

Code: Select all

# define UV_ADV
# define UV_QDRAG
# define UV_VIS2
# define UV_COR
# define RADIATION_2D

# define MASKING

# define VAR_RHO_2D
# define ASSUMED_SHAPE      
# define DOUBLE_PRECISION    
# define MPI                
# define NONLINEAR           
# define POWER_LAW           

# define ANA_INITIAL
# define ANA_DIAG
# define ANA_SMFLUX
# define ANA_STFLUX
# define ANA_SSFLUX

# define MIX_S_UV
# define MIX_GEO_TS
# define DEBUGGING

Some figures to ilustrate the problem:
Bathymetry

Zeta at step 9

V-bar at step 9

U-bar at step 9

Running the model with #ATM_PRESS, zeta at step 9.

Thanks in advance,

kate · #2 Unread post by **kate** » Fri Aug 23, 2019 1:58 pm

Why not try dt=0.5 or shorter? Your horizontal grid spacing is quite a bit smaller than anything I've tried outside of test problems, so it makes sense that you'd need a short timestep.

asantana · #3 Unread post by **asantana** » Fri Aug 23, 2019 3:09 pm

I did not try 0.5, because I just didn't know that it was possible.
I thought dt parameter was an integer.

I tested it with 0.5 and the max barotropic courant was higher than 1 (the model crashed - step 23).
But with dt <= 0.3 the Courant numbers are okay, so I adopted 0.25.
Now, I will let the model run during the weekend and evaluate the results on Monday.

Thanks for your answer kate!

ymamoutos · #4 Unread post by **ymamoutos** » Fri Aug 23, 2019 4:02 pm

Greetings,

how many barotropic steps you set (ndtfast)?

if ndtdast=1 try more i.e 10 for beginning and
check the value of maximum barotropic Courant number.
In any case - if i recall correctly - for stability
reasons it cannot exceed 0.8

Giannis

kate · #5 Unread post by **kate** » Fri Aug 23, 2019 4:18 pm

The OP said it was a barotropic model which requires NDTFAST = 1.

ymamoutos · #6 Unread post by **ymamoutos** » Fri Aug 23, 2019 4:25 pm

oops I missed that

I am apologizing. Please ignore my post

Giannis

asantana · #7 Unread post by **asantana** » Mon Aug 26, 2019 9:39 am

Updating, with dt equals 0.25 the model run without a problem until in an advanced where it blowed-up.
But this is not related to the original problem of this post.
Now, I need to adjust better the grid and bathymetry.

So, the original problem was solved!
Thanks!

Ocean Modeling Discussion

Barotropic model Blowing-up with high courant and dt = 1.

Barotropic model Blowing-up with high courant and dt = 1.

Re: Barotropic model Blowing-up with high courant and dt = 1

Re: Barotropic model Blowing-up with high courant and dt = 1

Re: Barotropic model Blowing-up with high courant and dt = 1

Re: Barotropic model Blowing-up with high courant and dt = 1

Re: Barotropic model Blowing-up with high courant and dt = 1

Re: Barotropic model Blowing-up with high courant and dt = 1