CHLNDDEV / OceanMesh2D

A two-dimensional triangular mesh generator with pre- and post-processing utilities written in pure MATLAB (no toolboxes required) designed specifically to build models that solve shallow-water equations or wave equations in a coastal environment (ADCIRC, FVCOM, WaveWatch3, SWAN, SCHISM, Telemac, etc.).
https://github.com/sponsors/krober10nd
GNU General Public License v3.0
178 stars 64 forks source link

Get stuck in using Cal_IT_Fric. #278

Closed Jiangchao3 closed 1 year ago

Jiangchao3 commented 1 year ago

Hi @WPringle and @krober10nd, I am testing my new study in Bangladesh, I am stuck in using the Cal_IT_Fric.

When I call the function of Cal_IT_Fric, the code just get stuck, no error is printed and can't run to the end, like this:

image

it seems stuck in the following line in the function (Gridded_to_Mesh_SeaBed_DepthAveraged) of Cal_IT_Fric: image

image

Actually, I also rerun my previous study's code to calculate the IT friction, it works very well, the function of Gridded_to_Mesh_SeaBed_DepthAveraged also works very well.

It is so odd.

Could you please give some suggestion to help me solve this problem? Many thanks.

The following is my msh file.

m_bangladesh.zip

krober10nd commented 1 year ago

Dear @Jiangchao3,

Perhaps you could pause the code when it appears "stuck" and see where it is? How long did you wait? This routine can become very costly.

Jiangchao3 commented 1 year ago

Hi @krober10nd , thanks for your kindly reply.

This time, I wait for about 5 minutes, pause it and enter the debug mode, find the code stuck in the following functions:

image

Calc_IT_Fric: image

Gridded_to_Mesh_SeaBed_DepthAveraged: image

fillmissing: image

I will continue to test the code.

Jiangchao3 commented 1 year ago

Hi @krober10nd, I think I found what the error lies:

in the function of Gridded_to_Mesh_SeaBed_DepthAveraged, it stuck in the iteration of zvalue = 82;

image

J is calculated as the following: image

the result from the 54 line is a 12*1 NaN vector: image image

it seems that line 58 and 59 can not work when the input is a full NaN vector: image

nothing can be fillmissing.

That may be why the Calc_IT_Fract stuck and can not move forward.

I also test the iteration zvalue = 81;

J: image

image

It's clear that the N_interp{zvalue}(J) is not a full NaN vector, so fillmissing in line 59 can work.

Hope the above description is clear to present the problem, thanks.

Finally, maybe the depth for each node is critical to calculate the IT fraction.

krober10nd commented 1 year ago

Thanks Jiang. I think it's clear what the problem is now. When the interpolation returns back all missing values, it becomes stuck in an infinite loop. It would be helpful if @WPringle took a look at this since he is the creator of this script.

Jiangchao3 commented 1 year ago

Yeah, waiting for @WPringle have a look!

WPringle commented 1 year ago

@Jiangchao3 I see, so I guess the problem is partly with N data being passed to the function, that there is too many NaNs at certain zcontour values. Perhaps then we need to fill NaNs in the N data in the vertical direction too which could be done just using nearest neighbor from the zcontour value above. Can you try this?

The reason for this can be that the depth of the N data at a certain location say goes to 1000 m but the mesh actually is at 1100 m, so at 1100 m the data has a NaN value but we want a value at the 1100 m contour on the mesh. A fix then is to extrapolate all the N values down the column using the nearest non-NaN value above. I think this may work better.

Jiangchao3 commented 1 year ago

@WPringle, thanks very much for your kindly reply. I am out of my office in recent days, I will check it after I come back!

Jiangchao3 commented 1 year ago

Dear @WPringle, today, I followed your suggestion to fix my Gridded_N_values.mat dataset. This is my code:

clearvars; clc; load('H:\QJC\8-Github\OceanMesh2D_Jch\datasets\Gridded_N_values.mat')

%%% bangladesh % lon from 60 to 110, (i from 960 to 1160) % lat from 0 to 30, (j from 361 to 501)

for i = 960:1160 for j = 361:501 temp_vertical = N(i,j,:); temp_vertical_new = fillmissing(temp_vertical,'previous'); N(i,j,:) = temp_vertical_new; end end

o_name = 'Gridded_N_values_new'; save([o_name '.mat'],'lon','lat','z','N');

Now, the Cal_IT_Fric function seems able to work well.

WPringle commented 1 year ago

Nice, thanks for letting me know! I would welcome a PR to add this change into the function

Jiangchao3 commented 1 year ago

Sure, I will create a function called Calc_IT_Fric_fix_N.m and then add it into the unitities folder and add a note into the Calc_IT_Fric.m to remind users who meet with similar problem to know how to fix it.

WPringle commented 1 year ago

Could you just make the change directly to `Gridded_to_Mesh_SeaBed_DepthAveraged'?

Jiangchao3 commented 1 year ago

`Gridded_to_Mesh_SeaBed_DepthAveraged'?

Ok, let me have a try.

Jiangchao3 commented 1 year ago

Hi @WPringle , I have created a PR https://github.com/CHLNDDEV/OceanMesh2D/pull/280, please review.