Blowing up problem

[ywang152]

Dear all,

I am working on a case study of simulation the water current of estuary without wind forcing. I have created my own model grid(by using GridBuilder), boundary and initial condition(by using ROMS matlab tool). The model blew up after 1st step. The standard output showing this:

Code: Select all

 Basin information for Grid 01:

 Maximum grid stiffness ratios:  rx0 =   2.000000E-01 (Beckmann and Haidvogel)
                                 rx1 =   5.855392E+00 (Haney)

 Initial domain volumes:  TotVolume =  9.9844783620E+07 m3
                         MinCellVol =  2.0501004235E+02 m3
                         MaxCellVol =  4.0155433950E+03 m3
                            Max/Min =  1.9587057048E+01

 NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 000000000002 - 000000008641)


 TIME-STEP YYYY-MM-DD hh:mm:ss.ss  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
                     C => (i,j,k)       Cu            Cv            Cw         Max Speed

         1 0001-01-01 00:00:01.00  5.573662E-02  1.564449E+02  1.565007E+02  4.488730E+07
                     (013,017,01)  1.373405E-02  9.012729E-03  0.000000E+00  7.432341E-01
      DEF_HIS     - creating  history      file, Grid 01: ocean_his.nc
         2 0001-01-01 00:00:02.00           NaN           NaN           NaN           NaN
                     (049,043,01)  6.110394E-03  7.331021E-03  3.991203E-02  7.416617E-01
 Found Error: 01   Line: 298      Source: ROMS/Nonlinear/main3d.F
 Found Error: 01   Line: 298      Source: ROMS/Drivers/nl_ocean.h

Below is my CPP options:

Code: Select all

 ESTUARYNW               Estuary current simulation without wind force
 ANA_BSFLUX              Analytical kinematic bottom salinity flux
 ANA_BTFLUX              Analytical kinematic bottom temperature flux
 ANA_SMFLUX              Analytical kinematic surface momentum flux
 ANA_SSFLUX              Analytical kinematic surface salinity flux
 ANA_STFLUX              Analytical kinematic surface temperature flux
 ASSUMED_SHAPE           Using assumed-shape arrays
 AVERAGES                Writing out time-averaged nonlinear model fields
 !BOUNDARY_ALLGATHER     Using mpi_allreduce in mp_boundary routine
 CANUTO_A                Canuto A-stability function formulation
 !COLLECT_ALL...         Using mpi_isend/mpi_recv in mp_collect routine
 CURVGRID                Orthogonal curvilinear grid
 DEBUGGING               Internal debugging switch activated
 DIAGNOSTICS_UV          Computing and writing momentum diagnostic terms
 DOUBLE_PRECISION        Double precision arithmetic numerical kernel.
 GLS_MIXING              Generic Length-Scale turbulence closure
 HDF5                    Creating NetCDF-4/HDF5 format files
 INLINE_2DIO             Processing 3D IO level by level to reduce memory needs
 MASKING                 Land/Sea masking
 MPI                     MPI distributed-memory configuration
 NONLINEAR               Nonlinear Model
 !NONLIN_EOS             Linear Equation of State for seawater
 N2S2_HORAVG             Horizontal smoothing of buoyancy and shear
 OUT_DOUBLE              Double precision output fields in NetCDF files
 PARALLEL_IO             Parallel I/O processing
 PERFECT_RESTART         Processing perfect restart variables
 POWER_LAW               Power-law shape time-averaging barotropic filter
 PRSGRD31                Standard density Jacobian formulation (Song, 1998)
 PROFILE                 Time profiling activated
 K_C2ADVECTION           Second-order centered differences advection of TKE fields
 RADIATION_2D            Use tangential phase speed in radiation conditions
 READ_WATER              Reading data at water points only
 REDUCE_ALLGATHER        Using mpi_allgather in mp_reduce routine
 RI_SPLINES              Parabolic Spline Reconstruction for Richardson Number
 RHO_SURF                Include difference between rho0 and surface density
 !RST_SINGLE             Double precision fields in restart NetCDF file
 SOLVE3D                 Solving 3D Primitive Equations
 TS_C4HADVECTION         Fourth-order centered horizontal advection of tracers
 TS_C4VADVECTION         Fourth-order centered vertical advection of tracers
 UV_ADV                  Advection of momentum
 UV_COR                  Coriolis term
 UV_C2ADVECTION          Second-order centered differences advection of momentum
 UV_QDRAG                Quadratic bottom stress
 VAR_RHO_2D              Variable density barotropic mode
 VISC_3DCOEF             Horizontal, time-dependent 3D viscosity coefficient

After reading some posts that having the same problem, I tried:
1)shorten my DT from 10 seconds to 1 second. The standard output is shown below. The maximum barotropic Courant Number is smaller than 1.

Code: Select all

 Minimum barotropic Courant Number =  1.27269811E-02
 Maximum barotropic Courant Number =  4.75991568E-02
 Maximum Coriolis   Courant Number =  9.50472208E-05

2)check if bathymetry requiring smoothing:
The standard output is shown below. It seems the rx0 and rx1 satisfying 0 < rx0 < 0.4 and 3 < rx1 < 7

Code: Select all

 Minimum X-grid spacing, DXmin =  4.01181852E-02 km
 Maximum X-grid spacing, DXmax =  4.01577087E-02 km
 Minimum Y-grid spacing, DYmin =  4.02366393E-02 km
 Maximum Y-grid spacing, DYmax =  4.02661984E-02 km
 Minimum Z-grid spacing, DZmin = -2.17046894E-01 m
 Maximum Z-grid spacing, DZmax =  2.48532885E+00 m
...
...
 Horizontal mixing scaled by grid size, GRDMAX =  4.02119170E-02 km

 Minimum horizontal viscosity coefficient =  0.00000000E+00 ▒▒▒
 Maximum horizontal viscosity coefficient =  0.00000000E+00 ▒▒▒
...
...
...
 Basin information for Grid 01:

[b] Maximum grid stiffness ratios:  rx0 =   2.000000E-01 (Beckmann and Haidvogel)
                                 rx1 =   5.855392E+00 (Haney)[/b]

3)check if boundary/initial/grid nc file has NAN values:
I didn't find nan value.

I tried to plot the restart file and come up with some questions..
1) The figures for temperature and salinity looks ok, same as my initial conditions. However, figures for u and v looks weird, it seems the mask doesn't work. I attached the screenshot of it.
The temperature and salinity in restart files is NaN in land area, while u and v is 0 on land area. Since I actived "READ_WATER" option, which said 'use if only reading water points data', all the value on land should be NaN in restart file (am I right?). I am a little confused about it.
2) I mentioned the dimension of some variables (temp, salt, u and v) in ocean_rst.nc file, is 'Dimensions: xi_v,eta_v,s_rho,two,ocean_time'. What does the 'two' mean? Is this because that I set NAT=2?
3) I saw in one post, people suggested to check the forcing file if ROMS blow up at 1st step. I don't have forcing file right now, only set five analytical options in header file. They are:

Code: Select all

#define ANA_SMFLUX
#define ANA_STFLUX
#define ANA_SSFLUX
#define ANA_BTFLUX
#define ANA_BSFLUX

Could this be the reason of blowing up? Do I have to create a ocean_frc.nc file for ROMS, and set the wind field to be zero?

Thanks for your attention and help!

Best regards
Yifan

[kate]

Minimum Z-grid spacing, DZmin = -2.17046894E-01 m

This would worry me. It means you need to be using WET_DRY or make the minimum depth deeper.

1)shorten my DT from 10 seconds to 1 second. The standard output is shown below. The maximum barotropic Courant Number is smaller than 1.

Code: Select all

 Minimum barotropic Courant Number =  1.27269811E-02
 Maximum barotropic Courant Number =  4.75991568E-02
 Maximum Coriolis   Courant Number =  9.50472208E-05

Seems OK.

I tried to plot the restart file and come up with some questions..
1) The figures for temperature and salinity looks ok, same as my initial conditions. However, figures for u and v looks weird, it seems the mask doesn't work. I attached the screenshot of it.
The temperature and salinity in restart files is NaN in land area, while u and v is 0 on land area. Since I actived "READ_WATER" option, which said 'use if only reading water points data', all the value on land should be NaN in restart file (am I right?). I am a little confused about it.

I think your u and v are fine. You don't want to mask with the land mask because river source can bring in non-zero flow from the land mask. The READ_WATER option is for only reading and writing water points, which I've never used because your plotting software needs to know how to handle that too.

2) I mentioned the dimension of some variables (temp, salt, u and v) in ocean_rst.nc file, is 'Dimensions: xi_v,eta_v,s_rho,two,ocean_time'. What does the 'two' mean? Is this because that I set NAT=2?

'two' means the number two, as in the number of time records of some fields in the restart file for perfect restarts.

3) I saw in one post, people suggested to check the forcing file if ROMS blow up at 1st step. I don't have forcing file right now, only set five analytical options in header file. They are:

Code: Select all

#define ANA_SMFLUX
#define ANA_STFLUX
#define ANA_SSFLUX
#define ANA_BTFLUX
#define ANA_BSFLUX

Could this be the reason of blowing up? Do I have to create a ocean_frc.nc file for ROMS, and set the wind field to be zero?

Nope, not the problem.

[ywang152]

Hi Kate,

Thanks for you reply! The minimum depth is the problem. I correct it to 2 meter instead and now ROMS is running!

Best regards
Yifan

[crperezz]

Hello Kate, I have been trying to simulate with roms but I always get a blowup in the first time steps, in my case it is the density, I have tried reducing the time step, the grid metrics seem to be ok, the courant numbers are also fine, Here is the metrics of my vertical levels, I do not know if I have a problem there, thank you very much:

Minimum X-grid spacing, DXmin = 9.79433866E-01 km
Maximum X-grid spacing, DXmax = 1.04156164E+00 km
Minimum Y-grid spacing, DYmin = 1.05177613E+00 km
Maximum Y-grid spacing, DYmax = 1.15377342E+00 km
Minimum Z-grid spacing, DZmin = 6.41100902E-02 m
Maximum Z-grid spacing, DZmax = 5.73437391E+02 m

Minimum barotropic Courant Number = 5.82491996E-03
Maximum barotropic Courant Number = 3.16808226E-01
Maximum Coriolis Courant Number = 1.89801340E-03

[kate]

That Courant number could still be high. Does the model run longer with a shorter step?

Did you smooth the bathymetry? Did you look at where the blow-up happens?

[crperezz]

Hello Kate, the model does not seem to run more longer with a smaller time step, I have smoothed the bathymetry several times, I think the model is blowingup near the west coast of the domain, here is an picture of the earth mask, maybe that is the problem ?, thank you very much.

[crperezz]

Excuse me the picture was not uploaded correctly, here I attach it, thank you very much

[kate]

I can't tell from that picture if the mask is a problem or not. You need to zoom in and edit the mask so that there are no enclosed lakes. Also, near the boundary, you want no bits of land leaving just one grid of water at the boundary, like:

bitmap.png (1.66 KiB) Viewed 5975 times

The edge of this grid is the edge of the computational domain, outside halo points not shown. if you do something like this, the advection terms are not well-behaved. Anyway, a picture of some field blowing up might be educational.

[crperezz]

Hello Kate, I have revised my earth mask and I have some problems, I am solving them, here I attach an picture from the restart file, in this case, the density anomalies, which was the blowup problem, I do not know if I indicate something to you, thank you very much.

PD: The image was obtained with PANOPLY, the NASA software to visualize netcdf, in reality the domain is longer and narrower

[kate]

I assume you are showing us bottom density, which really only tells something about the bottom depth. It looks like you still have lakes. Note that the blue square here is a lake, not a bay (red being land):

bitmap2.png (1.76 KiB) Viewed 5968 times