GAHM-based Perturbation Convergence

[FariborzDaneshvar-NOAA]

Default timelimit for prep.sbatch is 30 minutes. It timed out for Irma on Hercules. I increased the time limit to 2 hours, but it didn't help. This is what I got in slurm-***.setup.out

+ setup_ensemble --track-file /work2/noaa/nosofs/shared/nhc_hurricanes/irma_2017_aa7c2ce2-eaf9-46ec-8298-9e7562ab9ee4/nhc_track/hurricane-track.dat \
--output-directory /work2/noaa/nosofs/shared/nhc_hurricanes/irma_2017_aa7c2ce2-eaf9-46ec-8298-9e7562ab9ee4/setup/ensemble.dir/ \
--num-perturbations 29 \
--mesh-directory /work2/noaa/nosofs/shared/nhc_hurricanes/irma_2017_aa7c2ce2-eaf9-46ec-8298-9e7562ab9ee4/mesh/ \
--hires-region /work2/noaa/nosofs/shared/nhc_hurricanes/irma_2017_aa7c2ce2-eaf9-46ec-8298-9e7562ab9ee4/mesh/hires \
--sample-from-distribution \
--sample-rule korobov \
--date-range-file /work2/noaa/nosofs/shared/nhc_hurricanes/irma_2017_aa7c2ce2-eaf9-46ec-8298-9e7562ab9ee4/setup/dates.csv \
--nwm-file /work2/noaa/nos-surge/smani/data/nwm/NWM_v2.0_channel_hydrofabric/nwm_v2_0_hydrofabric.gdb \
--tpxo-dir /work2/noaa/nos-surge/smani/data/tpxo \
--variables cross_track along_track radius_of_maximum_winds max_sustained_wind_speed \
--pahm-model gahm irma 2017
/work2/noaa/nos-surge/daneshva/miniconda3/envs/stormworkflow/lib/python3.11/site-packages/pyschism/forcing/hycom/gofs.py:8: 
UserWarning: The seawater library is deprecated! Please use gsw instead.   
  import seawater as sw
can't find the 'adcircpy' module
can't find the 'adcircpy' module
can't find the 'adcircpy' module
/work2/noaa/nos-surge/daneshva/miniconda3/envs/stormworkflow/lib/python3.11/site-packages/stormevents/nhc/track.py:183: 
UserWarning: It is recommended to specify the file_deck and/or advisories when reading from file 
   warnings.warn(
slurmstepd: error: *** JOB 2426246 ON hercules-02-26 CANCELLED AT 2024-08-26T14:31:51 DUE TO TIME LIMIT ***

[SorooshMani-NOAA]

@FariborzDaneshvar-NOAA thanks for pointing this issue out. My initial assessment points to the ensembleperturbation/ensembleperturbation/perturbation/atcf.py(785)find_parameter_from_GAHM_profile() as the reason behind the workflow getting stuck. It seems that in some cases the while loop in this function does not converge.

I'll share more details when I learn more.

@WPringle do you have any suggestion of where to look to find the source of issues with convergence?

[FariborzDaneshvar-NOAA]

@SorooshMani-NOAA thanks for looking into this.

[SorooshMani-NOAA]

It is getting stuck when calculating 'max_sustained_wind_speed'. This is one of the instances where it has so far run 250k iterations and hasn't converged.

(Pdb) pp write_kwargs
{'filename': PosixPath('/home/smani/workarea2/runs/irma_2017_aa7c2ce2-fariborz/setup/ensemble.dir/track_files/vortex_4_variable_korobov_5.22'),
 'perturbation': {'along_track': -0.9674215661017008,
                  'cross_track': -0.967421566101701,
                  'max_sustained_wind_speed': -0.9674215661017016,
                  'radius_of_maximum_winds': -0.9674215661017016},
 'variables': [CrossTrack(lower_bound=<Quantity(-inf, 'nautical_mile')>, upper_bound=<Quantity(inf, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nm]
0               4.98
12             16.16
24              23.1
36             28.95
48             38.03
72             56.88
96             92.95
120           119.67, '50-95kt':      mean error [nm]
0               2.89
12             11.58
24             16.83
36              21.1
48             27.76
72             47.51
96             68.61
120           103.45, '>95kt':      mean error [nm]
0               1.85
12              7.79
24             12.68
36             17.92
48             25.01
72             40.48
96             60.69
120            79.98}, default=None, unit=<Unit('nautical_mile')>),
               AlongTrack(lower_bound=<Quantity(-inf, 'nautical_mile')>, upper_bound=<Quantity(inf, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nm]
0               6.33
12             17.77
24             26.66
36             37.75
48             51.07
72             69.22
96            108.59
120           125.01, '50-95kt':      mean error [nm]
0               3.68
12             12.74
24             19.43
36             27.51
48             37.28
72             57.82
96             80.15
120           108.07, '>95kt':      mean error [nm]
0               2.35
12              8.57
24             14.64
36             23.36
48             33.59
72             49.26
96              70.9
120            83.55}, default=None, unit=<Unit('nautical_mile')>),
               RadiusOfMaximumWinds(lower_bound=<Quantity(5, 'nautical_mile')>, upper_bound=<Quantity(200, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nmi]
0                3.38
12              11.81
24              15.04
36              17.44
48              17.78
60              17.78
72              18.84
84              19.34
96              21.15
108             23.19
120             24.77, '50-95kt':      mean error [nmi]
0                3.38
12               6.67
24               9.51
36               10.7
48              11.37
60              12.27
72              12.67
84              13.32
96              13.97
108             14.45
120             15.79, '>95kt':      mean error [nmi]
0                3.38
12               3.66
24               4.98
36               6.26
48               7.25
60               8.57
72               9.22
84              10.69
96              11.35
108             11.57
120             11.38}, default=None, unit=<Unit('nautical_mile')>),
               MaximumSustainedWindSpeed(lower_bound=<Quantity(15, 'knot')>, upper_bound=<Quantity(175, 'knot')>, historical_forecast_errors={'<50kt':      mean error [kt]
0               1.45
12              4.01
24              6.17
36              8.42
48             10.46
72             14.28
96             18.26
120            19.91, '50-95kt':      mean error [kt]
0               2.26
12              5.75
24              8.54
36              9.97
48             11.28
72             13.11
96             13.46
120            12.62, '>95kt':      mean error [kt]
0                2.8
12              7.94
24             11.53
36             13.27
48             12.66
72             13.41
96             13.46
120            13.55}, default=None, unit=<Unit('knot')>)],
 'weight': 1.0}

Where perturb is being called for maximum sustained wind with the following args:

(Pdb) p dataframe
   basin storm_number            datetime advisory_number advisory  ... isowave_radius_for_SWQ  isowave_radius_for_NWQ  extra_values                    geometry  track_start_time
0     AL           11 2017-09-10 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-80.90000 23.40000)        2017-09-10
1     AL           11 2017-09-10 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-80.90000 23.40000)        2017-09-10
2     AL           11 2017-09-10 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-80.90000 23.40000)        2017-09-10
3     AL           11 2017-09-10 03:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.00000 23.50000)        2017-09-10
4     AL           11 2017-09-10 03:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.00000 23.50000)        2017-09-10
5     AL           11 2017-09-10 03:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.00000 23.50000)        2017-09-10
6     AL           11 2017-09-10 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.70000 24.70000)        2017-09-10
7     AL           11 2017-09-10 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.70000 24.70000)        2017-09-10
8     AL           11 2017-09-10 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-81.70000 24.70000)        2017-09-10
9     AL           11 2017-09-11 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-82.40000 26.80000)        2017-09-10
10    AL           11 2017-09-11 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-82.40000 26.80000)        2017-09-10
11    AL           11 2017-09-11 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-82.40000 26.80000)        2017-09-10
12    AL           11 2017-09-11 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-83.40000 29.50000)        2017-09-10
13    AL           11 2017-09-11 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-83.40000 29.50000)        2017-09-10
14    AL           11 2017-09-11 12:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-83.40000 29.50000)        2017-09-10
15    AL           11 2017-09-12 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-85.00000 32.20000)        2017-09-10
16    AL           11 2017-09-12 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-85.00000 32.20000)        2017-09-10
17    AL           11 2017-09-13 00:00:00              03     OFCL  ...                    NaN                     NaN          <NA>  POINT (-89.00000 35.30000)        2017-09-10

[18 rows x 38 columns]
(Pdb) p perturbed_values
<Quantity([ -3.39488706  -3.39488706  -3.39488706  -4.95289773  -4.95289773
  -4.95289773  -9.62692973  -9.62692973  -9.62692973 -13.97965993
 -13.97965993 -13.97965993 -16.08933974 -16.08933974 -16.08933974
 -15.34973935 -15.34973935 -16.2590841 ], 'knot')>
(Pdb) p self.validation_times
0    0 days 00:00:00
1    0 days 00:00:00
2    0 days 00:00:00
3    0 days 03:00:00
4    0 days 03:00:00
5    0 days 03:00:00
6    0 days 12:00:00
7    0 days 12:00:00
8    0 days 12:00:00
9    1 days 00:00:00
10   1 days 00:00:00
11   1 days 00:00:00
12   1 days 12:00:00
13   1 days 12:00:00
14   1 days 12:00:00
15   2 days 00:00:00
16   2 days 00:00:00
17   3 days 00:00:00
Name: datetime, dtype: timedelta64[ns]

[SorooshMani-NOAA]

The issue happens calling find_parameter_from_GAHM_profile when

> /work2/noaa/nos-surge/smani/sandbox/ensembleperturbation/ensembleperturbation/perturbation/atcf.py(878)perturb()
-> isotach_rad_new = self.find_parameter_from_GAHM_profile(
(Pdb) p B
<Quantity([1.06355669 1.09404263 1.20258902 1.09446344 1.15370529 1.30797113
 1.02714703 1.21740032 1.54407791 0.84498247 1.19780192 2.1694523
 0.5840552  1.04473255 3.85651383 0.35032164 2.53739058 0.78149583], 'dimensionless')>
(Pdb) p Vr_new
<Quantity([33.47134769 51.17436185 66.52492666 35.98759579 53.71332687 69.19205481
 32.50896769 50.13027531 65.43037453 31.41401263 48.92169327 64.07228195
 29.36840886 46.55199002 61.48399824 28.48452583 44.78573501 30.91440999], 'knot')>
(Pdb) p Vmax_new
<Quantity([114.73681956 114.73681956 114.73681956 119.79958963 119.79958963
 119.79958963 125.94262436 125.94262436 125.94262436 123.77365661
 123.77365661 123.77365661  95.63147867  95.63147867  95.63147867
  65.99886785  65.99886785  36.75596061], 'knot')>
(Pdb) p f_new
<Quantity([5.80771564e-05 5.80771564e-05 5.80771564e-05 5.84144580e-05
 5.84144580e-05 5.84144580e-05 6.14717060e-05 6.14717060e-05
 6.14717060e-05 6.65791593e-05 6.65791593e-05 6.65791593e-05
 7.30774232e-05 7.30774232e-05 7.30774232e-05 7.94474118e-05
 7.94474118e-05 8.70895327e-05], '1 / second')>
(Pdb) p Roinv_new
<Quantity([0.0348013  0.0348013  0.0348013  0.03718388 0.03718388 0.03718388
 0.04118302 0.04118302 0.04118302 0.05623155 0.05623155 0.05623155
 0.09515583 0.09515583 0.09515583 0.18110131 0.18110131 0.8289199 ], 'dimensionless')>
(Pdb) p Bg
<Quantity([1.10107311 1.13260627 1.24488685 1.13571391 1.19713075 1.35707186
 1.07016126 1.26814073 1.60814624 0.89402072 1.26624446 2.29205274
 0.64498395 1.14727364 4.224365   0.43518941 3.00070197 1.48580393], 'dimensionless')>
(Pdb) p phi
<Quantity([1.03054375 1.02969338 1.02701523 1.03156676 1.02994728 1.02641777
 1.03696085 1.0311906  1.02459606 1.05954884 1.04204393 1.02322717
 1.13471327 1.07573427 1.02056828 1.35225848 1.0510989  1.30490981], 'dimensionless')>
(Pdb) p Rmax_new
<Quantity([19.09811367 19.09811367 19.09811367 21.18298584 21.18298584 21.18298584
 23.43760237 23.43760237 23.43760237 29.03804913 29.03804913 29.03804913
 34.5899975  34.5899975  34.5899975  41.79033257 41.79033257 97.1788781 ], 'nautical_mile')>
(Pdb) p Rrat
<Quantity([0.11234185 0.19098114 0.31830189 0.1246058  0.21182986 0.35304976
 0.10653456 0.21306911 0.39062671 0.0967935  0.24198374 0.58076098
 0.10809374 0.34589997 0.98828564 0.1816971  1.04475831 0.84107202], 'dimensionless')>

[SorooshMani-NOAA]

The issue is one item in the test for difference between velocity profiles:

(Pdb) p Vr
<Quantity([33.47134769 51.17436185 66.52492666 35.98759579 53.71332687 69.19205481
 32.50896769 50.13027531 65.43037453 31.41401263 48.92169327 64.07228195
 29.36840886 46.55199002 61.48399824 28.48452583 44.78573501 30.91440999], 'knot')>
(Pdb) p Vr_test
<Quantity([33.47134769 51.17436185 66.52492666 35.98759579 53.71332687 69.19205481
 32.50896769 50.13027531 65.43037453 31.41401263 48.92169327 64.07228195
 29.36840886 46.55199002 61.48399824 26.28496988 44.78573501 30.91440999], 'knot')>
(Pdb) p abs(Vr_test - Vr) > tol
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False,  True, False, False])

@WPringle do you have any suggestions for how to address it?

[WPringle]

I had come across this a few times that it didn’t converge but then I addressed it so all the tests were converging. It requires some debugging to work it out properly. Can you let me know how to repeat your example, like what storm and how many ensembles? Thanks!

Get Outlook for iOShttps://aka.ms/o0ukef

---

From: Soroosh Mani @.> Sent: Tuesday, August 27, 2024 10:36:41 AM To: oceanmodeling/ondemand-storm-workflow @.> Cc: Pringle, William @.>; Mention @.> Subject: Re: [oceanmodeling/ondemand-storm-workflow] prep.sbatch timed out on Hercules (Issue #75)

The issue is one item in the test for difference between velocity profiles: (Pdb) p Vr <Quantity([33. 47134769 51. 17436185 66. 52492666 35. 98759579 53. 71332687 69. 19205481 32. 50896769 50. 13027531 65. 43037453 31. 41401263 48. 92169327 64. 07228195 ZjQcmQRYFpfptBannerStart This Message Is From an External Sender This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

The issue is one item in the test for difference between velocity profiles:

(Pdb) p Vr <Quantity([33.47134769 51.17436185 66.52492666 35.98759579 53.71332687 69.19205481 32.50896769 50.13027531 65.43037453 31.41401263 48.92169327 64.07228195 29.36840886 46.55199002 61.48399824 28.48452583 44.78573501 30.91440999], 'knot')> (Pdb) p Vr_test <Quantity([33.47134769 51.17436185 66.52492666 35.98759579 53.71332687 69.19205481 32.50896769 50.13027531 65.43037453 31.41401263 48.92169327 64.07228195 29.36840886 46.55199002 61.48399824 26.28496988 44.78573501 30.91440999], 'knot')> (Pdb) p abs(Vr_test - Vr) > tol array([False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False])

@WPringlehttps://urldefense.us/v3/__https://github.com/WPringle__;!!G_uCfscf7eWS!buyz6VzvNJV-8Oew9jpFLQ0F6JkdHA-WTR3C6MVlL9OJaYss3OIw7dYRYdLI8uzqAuVVGFAW-yc8Dzh_j8S-bWwP$ do you have any suggestions for how to address it?

— Reply to this email directly, view it on GitHubhttps://urldefense.us/v3/__https://github.com/oceanmodeling/ondemand-storm-workflow/issues/75*issuecomment-2312900312__;Iw!!G_uCfscf7eWS!buyz6VzvNJV-8Oew9jpFLQ0F6JkdHA-WTR3C6MVlL9OJaYss3OIw7dYRYdLI8uzqAuVVGFAW-yc8Dzh_j0vpR_hX$, or unsubscribehttps://urldefense.us/v3/__https://github.com/notifications/unsubscribe-auth/AFBHFHQYWNHA6Q5YQVEDRZLZTSMITAVCNFSM6AAAAABNGKHUVSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJSHEYDAMZRGI__;!!G_uCfscf7eWS!buyz6VzvNJV-8Oew9jpFLQ0F6JkdHA-WTR3C6MVlL9OJaYss3OIw7dYRYdLI8uzqAuVVGFAW-yc8Dzh_jwIbjngt$. You are receiving this because you were mentioned.Message ID: @.***>

[SorooshMani-NOAA]

Shouldn't this line change as the following?

-  Vr_test[Vr_test < tol] = numpy.nan  # no solution
+  Vr_test[abs(Vr_test - Vr) < tol] = numpy.nan  # no solution

[SorooshMani-NOAA]

Can you let me know how to repeat your example

We're having trouble with perturbing Irma 2017 for 29 ensembles, using regression RMW and GAHM perturbation method

[WPringle]

No I don’t think so. The reason to add that line was that sometimes the Vr_test would just go to zero without converging because there is no possible solution.

Get Outlook for iOShttps://aka.ms/o0ukef

---

From: Soroosh Mani @.> Sent: Tuesday, August 27, 2024 10:46:13 AM To: oceanmodeling/ondemand-storm-workflow @.> Cc: Pringle, William @.>; Mention @.> Subject: Re: [oceanmodeling/ondemand-storm-workflow] prep.sbatch timed out on Hercules (Issue #75)

Shouldn't this line change as the following? - Vr_test[Vr_test < tol] = numpy. nan # no solution + Vr_test[abs(Vr_test - Vr) < tol] = numpy. nan # no solution — Reply to this email directly, view it on GitHub, or unsubscribe. You are ZjQcmQRYFpfptBannerStart This Message Is From an External Sender This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd

Shouldn't this line change as the following?

• Vr_test[Vr_test < tol] = numpy.nan # no solution
• Vr_test[abs(Vr_test - Vr) < tol] = numpy.nan # no solution

— Reply to this email directly, view it on GitHubhttps://urldefense.us/v3/__https://github.com/oceanmodeling/ondemand-storm-workflow/issues/75*issuecomment-2312922978__;Iw!!G_uCfscf7eWS!az4eKj54x23DPcJM7LER0In6T4lCidBS-q2FPoQhrHEGKrQQTduQe25Ypw4rvjYZVjq61pGHcIUVUCXAgXpIdoyW$, or unsubscribehttps://urldefense.us/v3/__https://github.com/notifications/unsubscribe-auth/AFBHFHSHE5UMTMW7XHHPRBTZTSNMLAVCNFSM6AAAAABNGKHUVSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJSHEZDEOJXHA__;!!G_uCfscf7eWS!az4eKj54x23DPcJM7LER0In6T4lCidBS-q2FPoQhrHEGKrQQTduQe25Ypw4rvjYZVjq61pGHcIUVUCXAgbfdoStL$. You are receiving this because you were mentioned.Message ID: @.***>

[FariborzDaneshvar-NOAA]

Can you let me know how to repeat your example, like what storm and how many ensembles?

@WPringle It was Hurricane Irma with 48 hour leadtime, 29 ensemble members, using GAHM and regression for RMW.

[SorooshMani-NOAA]

These are the args to perturb_tracks:

(Pdb) args
perturbations = 29
directory = PosixPath('/home/smani/workarea2/runs/irma_2017_aa7c2ce2-fariborz/setup/ensemble.dir/track_files')
storm = PosixPath('/home/smani/workarea2/runs/irma_2017_aa7c2ce2-fariborz/setup/ensemble.dir/track_to_perturb.dat')
variables = [CrossTrack(lower_bound=<Quantity(-inf, 'nautical_mile')>, upper_bound=<Quantity(inf, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nm]
0               4.98
12             16.16
24              23.1
36             28.95
48             38.03
72             56.88
96             92.95
120           119.67, '50-95kt':      mean error [nm]
0               2.89
12             11.58
24             16.83
36              21.1
48             27.76
72             47.51
96             68.61
120           103.45, '>95kt':      mean error [nm]
0               1.85
12              7.79
24             12.68
36             17.92
48             25.01
72             40.48
96             60.69
120            79.98}, default=None, unit=<Unit('nautical_mile')>), AlongTrack(lower_bound=<Quantity(-inf, 'nautical_mile')>, upper_bound=<Quantity(inf, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nm]
0               6.33
12             17.77
24             26.66
36             37.75
48             51.07
72             69.22
96            108.59
120           125.01, '50-95kt':      mean error [nm]
0               3.68
12             12.74
24             19.43
36             27.51
48             37.28
72             57.82
96             80.15
120           108.07, '>95kt':      mean error [nm]
0               2.35
12              8.57
24             14.64
36             23.36
48             33.59
72             49.26
96              70.9
120            83.55}, default=None, unit=<Unit('nautical_mile')>), RadiusOfMaximumWinds(lower_bound=<Quantity(5, 'nautical_mile')>, upper_bound=<Quantity(200, 'nautical_mile')>, historical_forecast_errors={'<50kt':      mean error [nmi]
0                3.38
12              11.81
24              15.04
36              17.44
48              17.78
60              17.78
72              18.84
84              19.34
96              21.15
108             23.19
120             24.77, '50-95kt':      mean error [nmi]
0                3.38
12               6.67
24               9.51
36               10.7
48              11.37
60              12.27
72              12.67
84              13.32
96              13.97
108             14.45
120             15.79, '>95kt':      mean error [nmi]
0                3.38
12               3.66
24               4.98
36               6.26
48               7.25
60               8.57
72               9.22
84              10.69
96              11.35
108             11.57
120             11.38}, default=None, unit=<Unit('nautical_mile')>), MaximumSustainedWindSpeed(lower_bound=<Quantity(15, 'knot')>, upper_bound=<Quantity(175, 'knot')>, historical_forecast_errors={'<50kt':      mean error [kt]
0               1.45
12              4.01
24              6.17
36              8.42
48             10.46
72             14.28
96             18.26
120            19.91, '50-95kt':      mean error [kt]
0               2.26
12              5.75
24              8.54
36              9.97
48             11.28
72             13.11
96             13.46
120            12.62, '>95kt':      mean error [kt]
0                2.8
12              7.94
24             11.53
36             13.27
48             12.66
72             13.41
96             13.46
120            13.55}, default=None, unit=<Unit('knot')>)]
sample_from_distribution = True
sample_rule = 'korobov'
quadrature = False
quadrature_order = 3
quadrature_rule = 'Gaussian'
sparse_quadrature = False
start_date = datetime.datetime(2017, 9, 10, 0, 0)
end_date = datetime.datetime(2017, 9, 13, 12, 0)
file_deck = 'a'
advisories = ['OFCL']
overwrite = True
parallel = False

Do you want me to also provide the track file being used?

[WPringle]

@SorooshMani-NOAA @FariborzDaneshvar-NOAA Thank you, I'll try it out.

[SorooshMani-NOAA]

Note to self, thanks to @WPringle: https://wiki.adcirc.org/Generalized_Asymmetric_Holland_Model

[SorooshMani-NOAA]

@WPringle I think I don't understand your bisection method. What specifically my understanding is that bisection method works on intervals and then choosing the new interval boundary based on having a 0 intersect between the new point and one of the old points. Here we're doing something a bit different. Can you please explain what is going on?

[WPringle]

@SorooshMani-NOAA You are referring to lines below especially to update alpha right.

        """ 
        Use root finding algorithm to return a desired parameter 
        based on the GAHM profile of velocity that matches the input Vr
        Input B (guess), isotach_rad and Rrat to get back the actual B
        Input Bg, phi, Rmax, and Rrat to get back isotach_rad 
        """
        if B is not None:
            # initial guesses when trying to find Bg, phi (and new B)
            Bg, phi = self.find_GAHM_parameters(B, Ro_inv)
            rfo2 = 0.5 * isotach_rad * f
            alpha = Rrat ** Bg
        Vr_test = 1e6 * MaximumSustainedWindSpeed.unit
        tol = 1e-2 * MaximumSustainedWindSpeed.unit
        while any(abs(Vr_test - Vr) > tol):
            if B is None:
                # updates for when trying to find isotach_rad
                isotach_rad = Rmax / Rrat
                rfo2 = 0.5 * isotach_rad * f
                alpha = Rrat ** Bg
            Vr_test = (
                numpy.sqrt(
                    Vmax ** 2 * (1 + Ro_inv) * numpy.exp(phi * (1 - alpha)) * alpha + rfo2 ** 2
                )
                - rfo2
            )
            Vr_test[Vr_test < tol] = numpy.nan  # no solution
            # bi-section method
            alpha[Rrat <= 1] *= 0.5 * (1 + (Vr / Vr_test)[Rrat <= 1] ** 2)
            alpha[Rrat > 1] *= 0.5 * (1 + (Vr_test / Vr)[Rrat > 1] ** 2)
            if B is not None:
                # update to and Bg, phi
                Bg = numpy.log(alpha) / numpy.log(Rrat)
                phi = 1 + Ro_inv / (Bg * (1 + Ro_inv))
                phi[phi < 1] = 1
            else:
                # update to Rrat
                Rrat = numpy.exp(numpy.log(alpha) / Bg)
        if B is not None:
            return Bg * phi / ((1 + Ro_inv) * numpy.exp(phi - 1))  # B
        else:
            return Rmax / Rrat  # isotach_rad

Bi-section method terminology is used loosely here but essentially we are updating alpha to be mid-way between the previous alpha and an estimate of the new one. The estimate of the new one is just a square of the ratio of the target Vr and Vr_test (Vg from eq (13) given alpha) multiplied by the previous alpha.

[SorooshMani-NOAA]

@WPringle yes I'm referring to this part of code about alpha; I want to see how to improve the implemented method, and I don't see where the half of square root of ratio factor addition comes from. Also what is the bracket we're considering. Is (0, alpha_ini] our bracket for root?

[WPringle]

@SorooshMani-NOAA The square is coming from the equation where alpha could be considered approximately proportional to Vr ** 2. So how much to change alpha is approximated to be how much deviation in Vr ** 2. The bracket is there because at Rrat = 1 we have a local maximum of Vr and so the root finding must be constrained to not jump between the inside curve (Rrat < 1) and the outside curve (Rrat > 1). In general this root finding seems to be working very well, so I'm not sure this is the part that needs to be changed...

[SorooshMani-NOAA]

This is just a test code to reproduce the issue:

import pdb
import tempfile
from datetime import datetime
from pathlib import Path

from stormevents.nhc.const import RMWFillMethod
from stormevents.nhc.track import VortexTrack
from ensembleperturbation.perturbation.atcf import (
    MaximumSustainedWindSpeed,
    perturb_tracks,
)

if __name__ == '__main__':

    irma_track_path = Path(tempfile.gettempdir()) / 'irma.dat'
    irma_pert_path = Path(tempfile.gettempdir()) / 'irma_perturb_dir'

    if not irma_track_path.exists():
        track_irma = VortexTrack.from_storm_name(
            'irma',
            2017,
            file_deck='a',
            advisories=['OFCL'],
            rmw_fill=RMWFillMethod.regression_penny_2023
        )
        track_irma.to_file(irma_track_path)

    pdb.runcall(
        perturb_tracks,
        39,
        variables=[MaximumSustainedWindSpeed],
        directory=irma_pert_path,
        storm=irma_track_path,
        sample_from_distribution=True,
        sample_rule='korobov',
        file_deck='a',
        advisories=['OFCL'],
        start_date=datetime(2017,9,10),
        end_date=datetime(2017,9,13,12),
        parallel=False,
    )

Note that I used pdb.runcall and parallel=False to be able to debug it and step into the functions.

[SorooshMani-NOAA]

In this specific case I see the value of calculated Vr_test for one of the rows keep jumping:

[SorooshMani-NOAA]

@WPringle I feel like the approach we're taking is not necessarily a sound one. The condition of the while loop checks $V_{g-test}(r)$ against the original track's $Vg(r)$, but then we modify the $\alpha$ and $\frac{R{max}}{r}$ without recalculating the test velocity for the next check, it might still work fine but the condition is always one loop behind.

By the way I still don't understand why we invert the Vr/Vr_test ratio based on Rrat? What assumptions do we have about Vr/Vr_test ratio?

[SorooshMani-NOAA]

@WPringle I gave it some more thought and I think I understand now what is happening, please correct me if I'm wrong. As I understand, to adjust alpha, you add half of alpha and half of a "correction factor", in this case Vr/Vr_test or its inverse (I still don't know why the inversion based on Rrat). And I also see a flip-flop in values.

That gave me the idea of using a "correcting factor" type of a fix to this. By that I mean to reduce the size of correction (centered around value of 0.5) such that the more times we go around the loop, the slower it tries to correct. This avoids the flip flop in the Irma case:

(Pdb) p i
1004
(Pdb) p Vr_test
<Quantity([33.47134769 51.24912547 66.6089342  36.2756504  54.04941179 69.52463968
 32.76801585 50.52182348 65.77775212 31.72150613 49.39156366 64.52820714
 29.92941066 47.29759352 62.24800562 29.4024715  46.05569623 31.47934678], 'knot')>
(Pdb) p Vr
<Quantity([33.47134769 51.24912547 66.6089342  36.2756504  54.04941179 69.52463968
 32.76801585 50.52182348 65.77775212 31.72150613 49.39156366 64.52820714
 29.92941066 47.29759352 62.24800562 29.39962196 46.05569623 31.47934678], 'knot')>

where i in this case is the number of iterations. This is my current fix:

rate = 1 / 2**(i//1000)
alpha[Rrat <= 1] *= 1 + (0.5 * (1 + (Vr / Vr_test)[Rrat <= 1] ** 2) - 1) * rate
alpha[Rrat > 1] *= 1 + (0.5 * (1 + (Vr_test / Vr)[Rrat > 1] ** 2) - 1) * rate     

Or simply

rate = 1 / 2**(i//1000)
alpha[Rrat <= 1] *= 1 + (0.5 * (Vr / Vr_test)[Rrat <= 1] ** 2 - 0.5) * rate
alpha[Rrat > 1] *= 1 + (0.5 * (Vr_test / Vr)[Rrat > 1] ** 2 - 0.5) * rate    

or actually why all the trouble, why not simply:

rate = 1 / 2**(i//1000)
alpha[Rrat <= 1] *= 1 + ((Vr / Vr_test)[Rrat <= 1] ** 2 - 1) * rate
alpha[Rrat > 1] *= 1 + ((Vr_test / Vr)[Rrat > 1] ** 2 - 1) * rate

but this last one doesn't work!! I think I'm rushing it ... but do you agree with my general approach? At least the first one that I know worked?!

[SorooshMani-NOAA]

I think I got my answer about the last one ... I should start from a factor of $\frac{1}{2}$ instead of 1 (of course!): rate = 1 / 2**(i//1000 + 1)

Now the question is, is this the right approch? This approach assumes that the higher the iteration the closer to the answer we are, so we reduce the rate to avoid jumping back and forth or diverging.

@WPringle can you send me the non-converging cases you had to try with this? If they converge and the current tests pass as well, I think we can go with this for now, right?

[WPringle]

@SorooshMani-NOAA Let me check this in more detail, I haven't got to it yet.

[SorooshMani-NOAA]

You can see my code in https://github.com/noaa-ocs-modeling/EnsemblePerturbation/tree/enhance/gahm_converge

[WPringle]

@SorooshMani-NOAA I went through and yes I agree the method may not be ideal but now seems to work. I made some small changes but can't push to the EnsemblePerturbation main directly. I work on adding this branch to my fork.

[SorooshMani-NOAA]

@WPringle thanks, please create a PR once you're ready!

[SorooshMani-NOAA]

@WPringle it still fails to converge for Idalia 2023, adding the following in the test cases fails:

('idalia', 2023, datetime(2023, 8, 29, 0), datetime(2023, 9, 7, 18)),

with

RuntimeError: GAHM function could not converge

[WPringle]

@SorooshMani-NOAA OK thanks for finding this! It's a stubborn problem.