Inverse CFL handling

[WPringle]

Mainly for functionality for SCHISM, it seems that it would be good to include the inverse of the current 'dt' function in edgefx and CheckTimestep functionality that will make the mesh finer where the CFL is too small (They suggest something like CFL > 0.5 everywhere for SCHISM so that the Lagrangian advection scheme stays accurate).

Ideally we can should be able to set a band so that CFL should be at least 0.5 and less than say 6 or whatever.

For example in CheckTimestep we can add points somewhere close to where the CFL is too small in addition to deleting points where CFL is too high like we do now.

[krober10nd]

Well at the mesh size function level this is straightforward since it entails adding a new name value pair and enforcing a minimum Cr and a maximum Cr. At the post processing level it’s a bit more difficult but perhaps splitting bars would be the way to go and then using a local direct smooth

[ivicajan]

That would be perfect if you could implement Cr_min and Cr_max options during the creation phase. Default values in explicit time stepping models like ADCIRC could be set to Cr = Cr_max and Cr_min to whatever default small value you like (as max. triangle size length should be stronger condition limiting that Cr_min).

Thanks in advance for your time and possible modification/extension of OceanMesh2D, that would make your tool directly applicable to SCHISM model community.

Cheers Ivica

[krober10nd]

I started a pull request as you can see above. With this branch you pass a negative timestep value for the dt name value pair to operate. By default, the option will bound the minimum Cr to 0.5 in the mesh size function and thus alter the final mesh vertex distribution.

[krober10nd]

As an example, I ran Example_JBAY with dt=-2.0 seconds and produced the following mesh.

[WPringle]

Interesting! It is much finer. Btw, I believe that for SCHISM they suggest to try to use something on the order of 100 sec time step pretty much as a standard choice. So it won't have such a large impact in that scenario. It's a good test case to see the effect though.

[ivicajan]

William is right, We run schism with dt 100-300 depending on application and if running 2d or 3d case. So what would be cook book to create mesh for constraint of dt=200s and cfl in interval btw 0.5 and 6?

Cheers Ivica

Sent from my mobile, possible typos

---

From: William Pringle notifications@github.com Sent: Friday, February 7, 2020 9:12:38 AM To: CHLNDDEV/OceanMesh2D OceanMesh2D@noreply.github.com Cc: Ivica Janekovic ivica.jan@gmail.com; Comment comment@noreply.github.com Subject: Re: [CHLNDDEV/OceanMesh2D] Inverse CFL handling (#51)

Interesting! It is much finer. Btw, I believe that for SCHISM they suggest to try to use something on the order of 100 sec time step pretty much as a standard choice. So it won't have such a large impact in that scenario. It's a good test case to see the effect though.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMD5OFQSNFNL5AQXE7LRBSYQNA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELBL45I#issuecomment-583188085, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEN4TMBFQDLNSWHDNEUJAT3RBSYQNANCNFSM4KM53XRA.

[krober10nd]

I will modify the commit so it by default bounds a maximum Courant number for a passed timestep (to preserve default behavior) but the user can optionally pass a minimum and maximum Cr for the passed timestep.

[ivicajan]

Great, just let me know when I can test it.

Thanks a lot!

Regards Ivica

Sent from my mobile, possible typos

---

From: keith roberts notifications@github.com Sent: Saturday, February 8, 2020 6:32:46 PM To: CHLNDDEV/OceanMesh2D OceanMesh2D@noreply.github.com Cc: Ivica Janekovic ivica.jan@gmail.com; Comment comment@noreply.github.com Subject: Re: [CHLNDDEV/OceanMesh2D] Inverse CFL handling (#51)

I will modify the commit so it by default bounds a maximum Courant number for a passed timestep (to preserve default behavior) but the user can optionally pass a minimum and maximum Cr for the passed timestep.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMGASNAWRIWDD66P2U3RB2C45A5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFOZVI#issuecomment-583724245, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEN4TMF6ULEGST3TFZLV55DRB2C45ANCNFSM4KM53XRA.

[krober10nd]

Ok, I made some progress on it with the latest commit on the #55 . You can see the technical details there. However, bounding the Courant number to 0.5 to 6.0 for a 100 s to 300-s timestep for the desired 15 m minimum resolution in the Example_JBAY proves to be pretty challenging!

The 15-m minimum element size produces a maximum estimated Courant number of 72 for a 100-sec timestep. In contrast, the minimum Courant number falls nicely within your bounds with this example with a value of 0.62.

What kind of mesh resolution do you normally use with these 100 to 300 second models?

[ivicajan]

Well, depending on coastline features and depth etc, but not less than 100m or so. Usually, small triangles have sides 300m. Cheers I

Sent from my mobile, possible typos

---

From: keith roberts notifications@github.com Sent: Saturday, February 8, 2020 7:56:46 PM To: CHLNDDEV/OceanMesh2D OceanMesh2D@noreply.github.com Cc: Ivica Janekovic ivica.jan@gmail.com; Comment comment@noreply.github.com Subject: Re: [CHLNDDEV/OceanMesh2D] Inverse CFL handling (#51)

Ok, I made some progress on it with the latest commit on the #55https://github.com/CHLNDDEV/OceanMesh2D/pull/55 . You can see the technical details there. However, bounding the Courant number to 0.5 to 6.0 for a 100 s to 300-s timestep for the desired 15 m minimum resolution in the Example_JBAY proves to be pretty challenging!

The 15-m minimum element size produces a maximum estimated Courant number of 72 for a 100-sec timestep. In contrast, the minimum Courant number falls nicely within your bounds with this example with a value of 0.62.

What kind of mesh resolution do you normally use with these 100 to 300 second models?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMH2QOHONG442GTCHBTRB2MX5A5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFQGPA#issuecomment-583729980, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEN4TMFLAY66MC2R2FE3P5LRB2MX5ANCNFSM4KM53XRA.

[krober10nd]

That resolution solution makes a lot of sense with a timestep of 100 to 300-s. I would suggest to to try out this short example (modified Example_ECGC.m)

% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region clearvars; clc; addpath(genpath('utilities/')); addpath(genpath('datasets/')); addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. %% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [100,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data %% to be used later with other OceanMesh classes... dem = 'SRTM15+V2.nc'; coastline = 'GSHHS_f_L1'; gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,... 'bbox',bbox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,... 'fs',R,'wl',wl,'max_el',max_el,... 'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options %% and build the mesh... mshopts = meshgen('ef',fh1,'bou',gdat1... 'plot_on',1,'proj','lam'); mshopts = mshopts.build;

m = mshopts.grd; Cr=CalcCFL(m); min(Cr) max(Cr)

Note the dt option at the top.

[ivicajan]

Thanks a lot! That way limited grid would be big improvements to anything currently available, I am not aware of any meshing tool capable of doing that. I like the elegance how you implemented it with dt arg and additional cr_min, cr_max and keeping default as well. Schism is super sensitive to grid and getting right grid is 90% of job. I am at conference on Monday, but will pull latest updates and try asap one of my regions at the north of Australia which is super sensitive to internal tides etc.

Thanks once again!

Regards, Ivica

Sent from my mobile, possible typos

---

From: keith roberts notifications@github.com Sent: Saturday, February 8, 2020 8:09:56 PM To: CHLNDDEV/OceanMesh2D OceanMesh2D@noreply.github.com Cc: Ivica Janekovic ivica.jan@gmail.com; Comment comment@noreply.github.com Subject: Re: [CHLNDDEV/OceanMesh2D] Inverse CFL handling (#51)

That resolution solution makes a lot of sense with a timestep of 100 to 300-s. I would suggest to to try out this short example (modified Example_ECGC.m)

% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region % with a high resolution inset around New York clearvars; clc; addpath(genpath('utilities/')); addpath(genpath('datasets/')); addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. %% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [100,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data %% to be used later with other OceanMesh classes... dem = 'SRTM15+V2.nc'; coastline = 'GSHHS_f_L1'; gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,... 'bbox',bbox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,... 'fs',R,'wl',wl,'max_el',max_el,... 'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options %% and build the mesh... mshopts = meshgen('ef',fh1,'bou',gdat1... 'plot_on',1,'proj','lam'); mshopts = mshopts.build;

m = mshopts.grd; Cr=CalcCFL(m); min(Cr) max(Cr)

Note the dt option at the top.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHZZBNOFHCZHG3YFODRB2OJJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFQN2Y#issuecomment-583730923, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEN4TMHDANP2CYXMAVCDYLLRB2OJJANCNFSM4KM53XRA.

[krober10nd]

Thank you Ivica, I appreciate your sentiments. I will continue to work on this pull request and let's also see what William thinks. Please share the tool with others if you think it can benefit their work.

[ivicajan]

Hi Keith,

this is still implemented in your InverseCFL branch which is not yet merged with default Projection? In other words I should pull specific this branch? Cheers Ivica

On Sat, Feb 8, 2020 at 10:30 PM keith roberts notifications@github.com wrote:

Thank you Ivica, I appreciate your sentiments. I will continue to work on this pull request and let's also see what William thinks. Please share the tool with others if you think it can benefit their work.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHHYSBSNQ23DGRAESTRB26YVA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFTAKA#issuecomment-583741480, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEN4TMD6F4ECADHHEG7JDILRB26YVANCNFSM4KM53XRA .

[krober10nd]

You can clone the repo and then checkout the branch. Or, you could just pull this specify branch.

On Sun, 9 Feb 2020 at 9:12 AM Ivica Janekovic notifications@github.com wrote:

Hi Keith,

this is still implemented in your InverseCFL branch which is not yet merged with default Projection? In other words I should pull specific this branch? Cheers Ivica

On Sat, Feb 8, 2020 at 10:30 PM keith roberts notifications@github.com wrote:

Thank you Ivica, I appreciate your sentiments. I will continue to work on this pull request and let's also see what William thinks. Please share the tool with others if you think it can benefit their work.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub < https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHHYSBSNQ23DGRAESTRB26YVA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFTAKA#issuecomment-583741480 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEN4TMD6F4ECADHHEG7JDILRB26YVANCNFSM4KM53XRA

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEOBZ7GLAKZVY7VU7RNU5NTRB7XMVA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGKTCQ#issuecomment-583838090, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOBZ7C2VERRJ6CQKVOEJJDRB7XMVANCNFSM4KM53XRA .

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

[ivicajan]

Yeah, just did it. Will give it a try and let you know. Cheers I

On Sun, Feb 9, 2020 at 8:20 PM keith roberts notifications@github.com wrote:

You can clone the repo and then checkout the branch. Or, you could just pull this specify branch.

On Sun, 9 Feb 2020 at 9:12 AM Ivica Janekovic notifications@github.com wrote:

Hi Keith,

this is still implemented in your InverseCFL branch which is not yet merged with default Projection? In other words I should pull specific this branch? Cheers Ivica

On Sat, Feb 8, 2020 at 10:30 PM keith roberts notifications@github.com wrote:

Thank you Ivica, I appreciate your sentiments. I will continue to work on this pull request and let's also see what William thinks. Please share the tool with others if you think it can benefit their work.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub <

https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHHYSBSNQ23DGRAESTRB26YVA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFTAKA#issuecomment-583741480

, or unsubscribe <

https://github.com/notifications/unsubscribe-auth/AEN4TMD6F4ECADHHEG7JDILRB26YVANCNFSM4KM53XRA

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub < https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEOBZ7GLAKZVY7VU7RNU5NTRB7XMVA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGKTCQ#issuecomment-583838090 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEOBZ7C2VERRJ6CQKVOEJJDRB7XMVANCNFSM4KM53XRA

.

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMG74U4FUQDOBRYQ4BDRB7YJPA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGKYPA#issuecomment-583838780, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEN4TME6MECEMZQ7D7PFWODRB7YJPANCNFSM4KM53XRA .

[ivicajan]

Hi Keith, Just tried something similar for my test case at the north west Aus but didn't work. It took some time as I used high resolution coastline. If you see something obvious that I did wrong, let me know, here is the history I did (for bathy I used GEBCO for testing): %%%%%%%%%%%%%%%%% addpath(genpath('utilities/')) addpath(genpath('datasets/')) addpath(genpath('m_map/'))

% get polyline x0, y0 tmp = [... 113.8367 -24.4197 111.9918 -23.6754 109.5319 -22.5589 106.6621 -20.9463 103.9289 -18.9201 102.4257 -15.8602 101.8107 -13.2138 102.0157 -11.0636 103.4506 -9.5337 104.8172 -7.9210 106.6800 -7.2530 113.9966 -7.9056 114.0696 -7.5029 115.5356 -7.5581 116.5781 -7.7051 118.0440 -7.7235 119.2494 -7.7418 120.6828 -7.8154 122.3117 -7.7051 123.9080 -7.5397 125.2436 -7.3192 126.6770 -6.8965 127.9801 -6.5105 129.3158 -5.8857 129.9999 -5.5916 130.3908 -5.7387 131.6939 -6.1246 133.0622 -6.7127 134.2675 -7.5213 135.0494 -8.2564 135.4729 -9.5245 135.3752 -10.7190 134.7082 -12.0038 134.6910 -12.2995 113.8367 -24.4197]; x0=tmp(:,1),y0=tmp(:,2); min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [300,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature. gdat1 = geodata('shp','GSHHS_h_L1','bbox',[x0 y0],'dem',' GEBCO_2014_2D100.0-50.0_170.0_0.0.nc','h0',min_el); fh1 = edgefx('geodata',gdat1,'fs',R,'max_el_ns',1e3,'max_el',max_el,'g',grade,'dt',dt);

mshopts_new=meshgen('ef',{fh1},'bou',{gdat1},'plot_on',1); mshopts_new = mshopts_new.build; m = mshopts_new.grd; m = interp(m,gdat1); % need depth to compute U for CFL Cr=CalcCFL(m); min(Cr), max(Cr) % I got %>> min(Cr) % %ans = % % 14.7828 % %>> max(Cr) % %ans = % % 1.0162e+03 % %>> dt % %dt = % % 300.0000 0.5000 6.0000 % %>>

Thanks, Ivica

On Sat, Feb 8, 2020 at 8:09 PM keith roberts notifications@github.com wrote:

That resolution solution makes a lot of sense with a timestep of 100 to 300-s. I would suggest to to try out this short example (modified Example_ECGC.m)

% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region % with a high resolution inset around New York clearvars; clc; addpath(genpath('utilities/')); addpath(genpath('datasets/')); addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. %% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [100,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data %% to be used later with other OceanMesh classes... dem = 'SRTM15+V2.nc'; coastline = 'GSHHS_f_L1'; gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,... 'bbox',bbox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,... 'fs',R,'wl',wl,'max_el',max_el,... 'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options %% and build the mesh... mshopts = meshgen('ef',fh1,'bou',gdat1... 'plot_on',1,'proj','lam'); mshopts = mshopts.build;

m = mshopts.grd; Cr=CalcCFL(m); min(Cr) max(Cr)

Note the dt option at the top.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHZZBNOFHCZHG3YFODRB2OJJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFQN2Y#issuecomment-583730923, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEN4TMHDANP2CYXMAVCDYLLRB2OJJANCNFSM4KM53XRA .

[krober10nd]

Looks good to me the code. I somewhat expected that since I’m enforcing the bounds of the Courant number in the mesh size function not the connectivity. We normally have a second step called CheckTimestep that edits mesh connectivity after it’s finished to guarantee the bound is obeyed. That’s already working for bounding the maximum courant number already so maybe that will help you out here

another thing you can do is check how many violations there are and see where they are on the mesh. Normally the mesh size function bounds dramatically reduce the number of violations but generally we always have to use a method afterward to edit the connectivity

On Sun, 9 Feb 2020 at 11:14 AM Ivica Janekovic notifications@github.com wrote:

Hi Keith, Just tried something similar for my test case at the north west Aus but didn't work. It took some time as I used high resolution coastline. If you see something obvious that I did wrong, let me know, here is the history I did (for bathy I used GEBCO for testing): %%%%%%%%%%%%%%%%% addpath(genpath('utilities/')) addpath(genpath('datasets/')) addpath(genpath('m_map/'))

% get polyline x0, y0 tmp = [... 113.8367 -24.4197 111.9918 -23.6754 109.5319 -22.5589 106.6621 -20.9463 103.9289 -18.9201 102.4257 -15.8602 101.8107 -13.2138 102.0157 -11.0636 103.4506 -9.5337 104.8172 -7.9210 106.6800 -7.2530 113.9966 -7.9056 114.0696 -7.5029 115.5356 -7.5581 116.5781 -7.7051 118.0440 -7.7235 119.2494 -7.7418 120.6828 -7.8154 122.3117 -7.7051 123.9080 -7.5397 125.2436 -7.3192 126.6770 -6.8965 127.9801 -6.5105 129.3158 -5.8857 129.9999 -5.5916 130.3908 -5.7387 131.6939 -6.1246 133.0622 -6.7127 134.2675 -7.5213 135.0494 -8.2564 135.4729 -9.5245 135.3752 -10.7190 134.7082 -12.0038 134.6910 -12.2995 113.8367 -24.4197]; x0=tmp(:,1),y0=tmp(:,2); min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [300,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature. gdat1 = geodata('shp','GSHHS_h_L1','bbox',[x0 y0],'dem',' GEBCO_2014_2D100.0-50.0_170.0_0.0.nc','h0',min_el); fh1 =

edgefx('geodata',gdat1,'fs',R,'max_el_ns',1e3,'max_el',max_el,'g',grade,'dt',dt);

mshopts_new=meshgen('ef',{fh1},'bou',{gdat1},'plot_on',1); mshopts_new = mshopts_new.build; m = mshopts_new.grd; m = interp(m,gdat1); % need depth to compute U for CFL Cr=CalcCFL(m); min(Cr), max(Cr) % I got %>> min(Cr) % %ans = % % 14.7828 % %>> max(Cr) % %ans = % % 1.0162e+03 % %>> dt % %dt = % % 300.0000 0.5000 6.0000 % %>>

Thanks, Ivica

On Sat, Feb 8, 2020 at 8:09 PM keith roberts notifications@github.com wrote:

That resolution solution makes a lot of sense with a timestep of 100 to 300-s. I would suggest to to try out this short example (modified Example_ECGC.m)

% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region % with a high resolution inset around New York clearvars; clc; addpath(genpath('utilities/')); addpath(genpath('datasets/')); addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. %% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [100,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data %% to be used later with other OceanMesh classes... dem = 'SRTM15+V2.nc'; coastline = 'GSHHS_f_L1'; gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,... 'bbox',bbox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,... 'fs',R,'wl',wl,'max_el',max_el,... 'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options %% and build the mesh... mshopts = meshgen('ef',fh1,'bou',gdat1... 'plot_on',1,'proj','lam'); mshopts = mshopts.build;

m = mshopts.grd; Cr=CalcCFL(m); min(Cr) max(Cr)

Note the dt option at the top.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub < https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHZZBNOFHCZHG3YFODRB2OJJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFQN2Y#issuecomment-583730923 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEN4TMHDANP2CYXMAVCDYLLRB2OJJANCNFSM4KM53XRA

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEOBZ7AYRITY463W2GQW5I3RCAFTJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGNQPI#issuecomment-583850045, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOBZ7EQ7GLFOIW7AB5UB4TRCAFTJANCNFSM4KM53XRA .

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

[krober10nd]

I am actively working on modifying CheckTimestep for this use case.

On Sun, 9 Feb 2020 at 11:23 AM Keith Roberts Keith.J.Roberts.133@nd.edu wrote:

Looks good to me the code. I somewhat expected that since I’m enforcing the bounds of the Courant number in the mesh size function not the connectivity. We normally have a second step called CheckTimestep that edits mesh connectivity after it’s finished to guarantee the bound is obeyed. That’s already working for bounding the maximum courant number already so maybe that will help you out here

another thing you can do is check how many violations there are and see where they are on the mesh. Normally the mesh size function bounds dramatically reduce the number of violations but generally we always have to use a method afterward to edit the connectivity

On Sun, 9 Feb 2020 at 11:14 AM Ivica Janekovic notifications@github.com wrote:

Hi Keith, Just tried something similar for my test case at the north west Aus but didn't work. It took some time as I used high resolution coastline. If you see something obvious that I did wrong, let me know, here is the history I did (for bathy I used GEBCO for testing): %%%%%%%%%%%%%%%%% addpath(genpath('utilities/')) addpath(genpath('datasets/')) addpath(genpath('m_map/'))

% get polyline x0, y0 tmp = [... 113.8367 -24.4197 111.9918 -23.6754 109.5319 -22.5589 106.6621 -20.9463 103.9289 -18.9201 102.4257 -15.8602 101.8107 -13.2138 102.0157 -11.0636 103.4506 -9.5337 104.8172 -7.9210 106.6800 -7.2530 113.9966 -7.9056 114.0696 -7.5029 115.5356 -7.5581 116.5781 -7.7051 118.0440 -7.7235 119.2494 -7.7418 120.6828 -7.8154 122.3117 -7.7051 123.9080 -7.5397 125.2436 -7.3192 126.6770 -6.8965 127.9801 -6.5105 129.3158 -5.8857 129.9999 -5.5916 130.3908 -5.7387 131.6939 -6.1246 133.0622 -6.7127 134.2675 -7.5213 135.0494 -8.2564 135.4729 -9.5245 135.3752 -10.7190 134.7082 -12.0038 134.6910 -12.2995 113.8367 -24.4197]; x0=tmp(:,1),y0=tmp(:,2); min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [300,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature. gdat1 = geodata('shp','GSHHS_h_L1','bbox',[x0 y0],'dem',' GEBCO_2014_2D100.0-50.0_170.0_0.0.nc','h0',min_el); fh1 =

edgefx('geodata',gdat1,'fs',R,'max_el_ns',1e3,'max_el',max_el,'g',grade,'dt',dt);

mshopts_new=meshgen('ef',{fh1},'bou',{gdat1},'plot_on',1); mshopts_new = mshopts_new.build; m = mshopts_new.grd; m = interp(m,gdat1); % need depth to compute U for CFL Cr=CalcCFL(m); min(Cr), max(Cr) % I got %>> min(Cr) % %ans = % % 14.7828 % %>> max(Cr) % %ans = % % 1.0162e+03 % %>> dt % %dt = % % 300.0000 0.5000 6.0000 % %>>

Thanks, Ivica

On Sat, Feb 8, 2020 at 8:09 PM keith roberts notifications@github.com wrote:

That resolution solution makes a lot of sense with a timestep of 100 to 300-s. I would suggest to to try out this short example (modified Example_ECGC.m)

% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region % with a high resolution inset around New York clearvars; clc; addpath(genpath('utilities/')); addpath(genpath('datasets/')); addpath(genpath('m_map/'));

%% STEP 1: set mesh extents and set parameters for mesh. %% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox min_el = 300; % minimum resolution in meters. max_el = 50e3; % maximum resolution in meters. wl = 30; % 60 elements resolve M2 wavelength. dt = [100,0.5,6.0]; % Set timestep to respect Courant bounds. grade = 0.35; % mesh grade in decimal percent. R = 3; % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data %% to be used later with other OceanMesh classes... dem = 'SRTM15+V2.nc'; coastline = 'GSHHS_f_L1'; gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,... 'bbox',bbox);

%% STEP 3: create an edge function class fh1 = edgefx('geodata',gdat1,... 'fs',R,'wl',wl,'max_el',max_el,... 'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options %% and build the mesh... mshopts = meshgen('ef',fh1,'bou',gdat1... 'plot_on',1,'proj','lam'); mshopts = mshopts.build;

m = mshopts.grd; Cr=CalcCFL(m); min(Cr) max(Cr)

Note the dt option at the top.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub < https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHZZBNOFHCZHG3YFODRB2OJJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELFQN2Y#issuecomment-583730923 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEN4TMHDANP2CYXMAVCDYLLRB2OJJANCNFSM4KM53XRA

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEOBZ7AYRITY463W2GQW5I3RCAFTJA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGNQPI#issuecomment-583850045, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOBZ7EQ7GLFOIW7AB5UB4TRCAFTJANCNFSM4KM53XRA .

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

[krober10nd]

Ah, wait a minute check out the function [out1,barlen,bars] = CalcCFL(obj,dt,type) in msh.m.

It's hardcoded for a CFL of 1.0 and if you don't pass it a timestep, it will gather one based a CFL condition of 1.

I will need to modify CalcCFL to see how well the mesh size function did.

[ivicajan]

If I pass dt of 300s then it is in the range 0.3 - 20, and there are ~10% of bad;

Cr=CalcCFL(m,300); min(Cr)

ans =

0.2952

max(Cr)

ans =

20.2939

bad=find(Cr>6);size(bad)

ans =

    7816           1

bad=find(Cr<0.5);size(bad)

ans =

    1545           1

size(m.p)

ans =

   91117           2

On Sun, Feb 9, 2020 at 10:34 PM keith roberts notifications@github.com wrote:

Ah, wait a minute check out the function [out1,barlen,bars] = CalcCFL(obj,dt,type) in msh.m.

It's hardcoded for a CFL of 1.0 and if you don't pass it a timestep, it will gather one based a CFL condition of 1.

I will need to modify CalcCFL to see how well the mesh size function did.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMGQDRYANNEWVGZXSA3RCAH7HA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGOAHQ#issuecomment-583852062, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEN4TMGUOVC4EHUJOZXJXVLRCAH7HANCNFSM4KM53XRA .

[krober10nd]

I see. And where are those points located? Would it be possible to share the mesh through email (with the bathy on?)?

In the mean time you could do is overshoot the timestep by 5-10% (I.e., build a mesh for 330 s)...

[ivicajan]

Sure, I can send it on github or you want to your private email? This is grid I made with old version of OceanMesh2d and was able to run it with schism and dt 300s, though it is not optimized grid for CFL.

Cheers I

On Sun, Feb 9, 2020 at 11:00 PM keith roberts notifications@github.com wrote:

I see. And where are those points located? Would it be possible to share the mesh through email (with the bathy on?)?

In the mean time you could do is overshoot the timestep by 5-10% (I.e., build a mesh for 330 s)...

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHM547OXICYPHV5QODRCAK7BA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGOUTA#issuecomment-583854668, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEN4TMGRQUGR4DG6YY5KPHTRCAK7BANCNFSM4KM53XRA .

[krober10nd]

Here is fine otherwise krober@usp.br works for me

On Sun, 9 Feb 2020 at 12:09 PM Ivica Janekovic notifications@github.com wrote:

Sure, I can send it on github or you want to your private email? This is grid I made with old version of OceanMesh2d and was able to run it with schism and dt 300s, though it is not optimized grid for CFL.

Cheers I

On Sun, Feb 9, 2020 at 11:00 PM keith roberts notifications@github.com wrote:

I see. And where are those points located? Would it be possible to share the mesh through email (with the bathy on?)?

In the mean time you could do is overshoot the timestep by 5-10% (I.e., build a mesh for 330 s)...

— You are receiving this because you commented. Reply to this email directly, view it on GitHub < https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEN4TMHM547OXICYPHV5QODRCAK7BA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGOUTA#issuecomment-583854668 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEN4TMGRQUGR4DG6YY5KPHTRCAK7BANCNFSM4KM53XRA

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/CHLNDDEV/OceanMesh2D/issues/51?email_source=notifications&email_token=AEOBZ7DBGEX2U2LDH7MH7RDRCAMCRA5CNFSM4KM53XRKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELGO4ZY#issuecomment-583855719, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOBZ7AXIOLAIJXO4QD3T6DRCAMCRANCNFSM4KM53XRA .

--

Keith Roberts Ph.D. (Eng.) University of São Paulo krober@usp.br krober10@nd.edu

[krober10nd]

Note some progress has been made with the new pr #56

[krober10nd]

PR #56 has been merged into the main branch (Projection)