IPE Physics: IPE stability issue

[twfang]

Under certain WAM conditions, WAM-IPE coupled run produces strange feature in O+ density, velocity, and eventually molecular ions also become NaNs. To debug this issue, we will run the model with time-vary drivers (timeobs) with empirical E fields for the 2013 St. Patrick Day storm period using 20130316 WAM initial condition under /scratch3/NCEPDEV/swpc/noscrub/Adam.Kubaryk/ICS/WAM-IPE/2013031600

and IPE initial condition under /scratch3/NCEPDEV/swpc/noscrub/refactored_ipe_input_decks/IPE_State.apex.201303160000.h5

I will then report TEC, O+ density, O+ velocity from this run to provide information to debug. Everyone who is doing the same run on different systems feels free to report plots here.

[twfang]

It appears that the high O+ velocity starts to show up in IPE within the first hour of model run. After removing the line of code with if statement in IPE_Plasma_Class.F90 that Tim suggested,

if(grid % altitude(i,lp).ge.200000.0_prec) then

The problem at 200 km seems to be resolved. Below shows the comparisons of O+ density and velocity between the runs with and without the fix at the 1hr (03160100).

No fix: O plus density

O plus velocity

With new fix: O plus density

O plus velocity

It is clear that the blobs in O+ velocity at around 200km disappeared after the fix. The strange density cutoff also disappeared. These odd features were not there in the IPE initial condition. I have submitted a 7 day run for both cases and will update the results tomorrow.

[wanghj2015]

I think there are still a clear discontinuity near 150km in both O+ density and velocity. That may be related to the conditional 'transport_min_altitude'; you need to remove that too.

Better start with a short run test before longer runs; examine the other fields too.

Another clear 'feature'/discontinuity is in the 'vertical'; and I think that is related to this parameter in IPE_Plasma_Class.F90 (line 114):

INTEGER, PARAMETER , PRIVATE :: perp_transport_max_lp = 151

! perp_transport_max_lp - Maximum of the lp-loop in ! Cross_Flux_Tube_Transport. Plasma properties are ! only advected on magnetic flux tubes with lp ! between 1 and perp_transport_max_lp ( inclusive ).

[twfang]

Tubes 151 to 170 have apex height below 150 km. Transport should not be happening below 150 km. This setting is not a bug and is necessary, but can be revisited.

[wanghj2015]

No one said it is a bug. Transport stops at 200km is not a bug either, it's a new 'feature' of refactored ipe, unless it is a typo though.

Houjun Wang

On Mon, Sep 23, 2019 at 6:04 PM twfang notifications@github.com wrote:

Tubes 151 to 170 have apex height below 150 km. Transport should not be happening below 150 km. This setting is not a bug and is necessary, but can be revisited.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AC7MMJ6UFXV5N4VRRCHUZILQLFKQ7A5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7MUBXI#issuecomment-534331613, or mute the thread https://github.com/notifications/unsubscribe-auth/AC7MMJ6HKSRHSJYBESGCJTLQLFKQ7ANCNFSM4IZRF3VQ .

[twfang]

The plot is from 90 km.

[timfullerrowell]

Houjun, I think you are looking at a grid resolution change, not a discontinuity. It all looks good in Tzu-Wei's plot so far. Tim

On Mon, Sep 23, 2019 at 5:21 PM wanghj2015 notifications@github.com wrote:

I think there are still a clear discontinuity near 150km in both O+ density and velocity. That may be related to the conditional 'transport_min_altitude'; you need to remove that too.

Better start with a short run test before longer runs; examine the other fields too.

Houjun Wang

On Mon, Sep 23, 2019 at 5:01 PM twfang notifications@github.com wrote:

It appears that the high O+ velocity starts to show up in IPE in the first hour of model run. After removing the line of code with if statement in IPE_Plasma_Class.F90 that Tim suggested,

if(grid % altitude(i,lp).ge.200000.0_prec) then

The problem at 200 km seems to be resolved. Below shows the comparisons of O+ density and velocity between the runs with and without the fix at the 1hr (03160100).

No fix: O plus density [image: oplus_03160100_nofix] < https://user-images.githubusercontent.com/22968399/65469122-746c7f80-de23-11e9-8cd3-f2aa70734a49.png

O plus velocity [image: 03160100_nofix] < https://user-images.githubusercontent.com/22968399/65469130-7c2c2400-de23-11e9-977e-9709e00b3ad4.png

With new fix: O plus density [image: oplus_03160100_fix] < https://user-images.githubusercontent.com/22968399/65469152-8817e600-de23-11e9-8a3a-14e77caffd93.png

O plus velocity [image: 03160100_fix] < https://user-images.githubusercontent.com/22968399/65469160-9108b780-de23-11e9-8a2c-d6d57db53762.png

It is clear that the blobs in O+ velocity disappeared after the fix. The strange density cutoff also disappeared. These odd features were not there in the IPE initial condition. I have submitted a 7 day run for both cases and will update the results tomorrow.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub < https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AC7MMJ2Y6UBYUFDAUC5CZC3QLFDFNA5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7MQVDQ#issuecomment-534317710 , or mute the thread < https://github.com/notifications/unsubscribe-auth/AC7MMJ2NKCXN4AHMPVYDYPDQLFDFNANCNFSM4IZRF3VQ

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AH5BFOB2FOGXAHXDDDTWKM3QLFFQLA5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7MR2EQ#issuecomment-534322450, or mute the thread https://github.com/notifications/unsubscribe-auth/AH5BFOFL3KUMHVEPSAHRKH3QLFFQLANCNFSM4IZRF3VQ .

[wanghj2015]

OK. But there is still another feature near 120km which is still discernible in ion velocity fields, which could become more problematic in time.

[twfang]

After implementing the fix (removing the if statement in IPE_Plasma_Class.F90), no NaN in molecular ions appears after a 4-day run. However, different problems already show up in IPE within the first few hour run.

The E region O+ density becomes very large:

The field-aline O+ velocity increases at some tubes below 120km (12 UT is shown here):

Very large NO+ density shows up below 120km (12 UT is shown here):

The large O+ ion velocity and molecular ion density didn't show up in the default run (with no fix):

A clear cut off at 120 km is shown in ionospheric parameters, but is not there in thermospheric parameters.

[wanghj2015]

@twfang would you send me your rundir/rotdir so that I can take a look at the output also? Thanks.

[twfang]

They are under /scratch4/NCEPDEV/stmp3/Tzu-Wei.Fang/development_ipe_refactor_merge /scratch4/NCEPDEV/stmp3/Tzu-Wei.Fang/WAMIPE_TECvalidation_withfix

[wanghj2015]

@twfang You may have missed, but there are a lot 'nan's:

cd /scratch4/NCEPDEV/stmp3/Tzu-Wei.Fang/WAMIPE_TECvalidation_withfix [Houjun.Wang@hfe09 WAMIPE_TECvalidation_withfix]$ h5dump -d "/apex/no_plus_density" IPE_State.apex.201303200000.h5 | grep nan (16,26,14): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,24): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,34): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,44): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,54): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,64): -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, -nan, (16,26,74): -nan, -nan, -nan, -nan, -nan, -nan, -nan, 0, 0, 0, 0, 0, 0, 0,

[twfang]

You are right, I have not looked at results beyond 18th since there are problems in the first day already.

[naomimaruyama]

@twfang ,

great news that no more NaNs are generated.

the increase in O+ @100km appears at high latitude at nighttime at 3UT. Is aurora included in your run? it would be interesting to see what is happening in O+ @100km between 0UT and 3UT.

[wanghj2015]

@IonospherePlasmasphereElectrodynamics There are 'NaN' == 'nan's.

[twfang]

No aurora is included in this run.

[timfullerrowell]

It would be nice to find the place in the code that refers to 120 or 125 km, which ever the discontinuity is. I couldn't find it in the cross tube transport. Anyone know where this altitude crops up in the code? We didn't have this issue with old IPE did we Naomi, after you redid the chemical equilibrium calculation, after the transport?

In addition, we need to introduce a vertical profile in the cross tube drift for the initial condition for the O+ after transport. We have the information to calculate the actual plasma velocity. Then calculate the ratio of actual perpendicular velocity Vi to ExB/B**2, say R. Then initialize the O+ after transport to: O+ with ExB transport (R) + O+ profile without transport (1-R) We can test with a simple linear change from 1 to zero.

[wanghj2015]

In flip/RSLPSD.f

 ZLBDY=120.        !.. Lower boundary altitude

But this may not be the reason for NaNs in molecular ions.

[twfang]

Anyone has old IPE output with O+ velocity? Want to see if the 120 km boundary was there before. The setting Houjun showed is also there in the old IPE.

[akubaryk]

Don't think ion velos were part of the traditional legacy output, but I could compile and run that on development if it'd help.

[twfang]

@IonospherePlasmasphereElectrodynamics Have you looked at O+ ion velocity in old IPE before? Anything that we can have a quick look? Either coupled run or standalone IPE output will work.

[akubaryk]

Actually it does look like the ion velocities are in IPE output.

from scipy.io import FortranFile

b = FortranFile('/scratch4/NCEPDEV/stmp3/Naomi.Maruyama/v1.varf107kp.v2.6sn3/ipe_grid_plasma_params.201303180000')
counter = 0 ; a = []
while True:
  try: a.append(b.read_record(dtype='f4')) ; counter +=1
  except: break

print counter
print a[14].max()

Gives 16 fields (9 ion densities, electron temperature, ion temperature, and 4(?) ion velocities): the max of the 15th field is ~8438, which seems reasonable. If we can confirm the order of those last four fields, I can rig up the post for this unless @wanghj2015 completed that work for the legacy IPE post already -- let me know.

[wanghj2015]

I download some output from runs with old ipe:

/scratch4/NCEPDEV/stmp4/Houjun.Wang/prwam.rad/*ipe_grid_plasma_params* and *ipe_grid_neutral_params*

I think ion velocities are there too:

j_loop1: DO jth=1,ISTOT !=(ISPEC+3+ISPEV)

mp_loop1:DO mp=1,mpstop
  lp_loop1:DO lp=1,nlp

    IN=JMIN_IN(lp)
    IS=JMAX_IS(lp)
    npts = IS-IN+1 
    dumm(JMIN_ING(lp):JMAX_ISG(lp),mp) = plasma_3d(IN:IS,lp,mp,jth)

  ENDDO lp_loop1!lp
ENDDO mp_loop1!mp

! ghgm - also write the 16 plasma files to a single unit.....
write (unit=5999) (dumm(:,mp),mp=1,mpstop)

ENDDO j_loop1

Not sure Adam can handle this by modifying your post.py? George's old post-proc produce only tec etc.

[naomimaruyama]

@twfang ,

please let me have a quick look at o+ velocity in some old ipe results. will keep you updated shortly.

[wanghj2015]

@akubaryk If you could handle the old output using your post.py, that would be great. I haven't done anything along that line.

[akubaryk]

Thanks for the update.

[akubaryk]

With @IonospherePlasmasphereElectrodynamics's help, I have a legacy post going.

Doesn't take long to run these on Hera, so I picked an arbitrary stmp3 directory of Naomi's that runs from 2013031600-2013031800: /scratch4/NCEPDEV/stmp4/Adam.Kubaryk/opsdat/v1.varf107kp.v2.6sn3/netcdf

If we want to do fresh runs of coupled development, that can be done, or if I picked the wrong directory to convert, let me know.

We should probably put the legacy post in the development branch and tag it prior to NOAA-SWPC/WAM-IPE#336 coming in just so it exists and people can easily access it.

[naomimaruyama]

@akubaryk,

Thanks for converting the old output. I am going to take a look. But I am really not sure if this run is the appropriate one to be a reference. It might be one of the old failed runs. I might ask for your help again to convert different runs later. Thanks.

[twfang]

I plotted out the old IPE netcdf file that Adam provided yesterday. Even though the run might not be the best one and several longitudes showed extremely large O+ velocity, but the 120 km discontinuity can be clearly seen. Thus, this feature does not only exist in the refactored IPE.

Results from /scratch4/NCEPDEV/stmp4/Adam.Kubaryk/opsdat/v1.varf107kp.v2.6sn3/netcdf/IPE_Params.geo.201303160100.nc4

Plots show oplus velocity in mp=30, mp=50, mp=70

[twfang]

There is no boundary at 120 km in O+ density, only in the velocity.

[twfang]

After implementing a hyperbolic tangent function for gradually decreasing ExB impact on atomic ion density/velocity, ion temperature and electron density below 200km, the O+ velocity seems to be a lot more reasonable. However, the implementation still shows the incorrect impact on the E-region O+ density and molecular ion density. Results from 6UT are shown here (03160600)

O+ density

O+ velocity

NO+ density

TEC and O+ density

[timfullerrowell]

One step forward.

It would be good to understand the 120 km boundary in the plots. Has anyone found any other reference apart from this one in Subroutine RSLPSD.f

dbg20121130 ZLBDY=115.!ZLBDY_flip !.. Lower boundary altitude ZLBDY=120. !.. Lower boundary altitude

IF(Z(JEQ).LE.ZLBDY+1.0) THEN DO J=JMIN,JMAX CALL HOEQ(FLDIM,J,N,TI) XIONN(1,J)=N(1,J) XIONN(2,J)=N(2,J) XIONV(1,J)=0.0 XIONV(2,J)=0.0 ENDDO

Does anyone know what HOEQ does? I am assuming it sets O+ and H+ to chemical equilibrium. Bot isn't there another place in the code that does that. Naomi, can you point be to the rcode that has the call to chemical equilibrium after the call to cross field transport.

I'm not sure I understand the 2D plots of O+ at 100 km, are they bad or not? That are the real numbers here? The O+ latitude/height plot appears to show very small O+ at 100 km.

Tim

On Thu, Sep 26, 2019 at 3:28 PM twfang notifications@github.com wrote:

After implementing a hyperbolic tangent function for gradually decreasing ExB impact on atomic ion density/velocity, ion temperature and electron density below 200km, the O+ velocity seems to be a lot more reasonable. However, the implementation still shows the incorrect impact on the E-region O+ density and molecular ion density. Results from 6UT are shown here (03160600)

O+ density [image: oplus_03160600_newfix] https://user-images.githubusercontent.com/22968399/65726314-09fd4e80-e072-11e9-96db-7f288c8e11a6.png

O+ velocity [image: voplus_03160600_newfix] https://user-images.githubusercontent.com/22968399/65726334-1386b680-e072-11e9-97d8-5e96df1be8b1.png

NO+ density [image: noplus_03160600_newfix] https://user-images.githubusercontent.com/22968399/65726439-50eb4400-e072-11e9-9228-32f2538f4cc5.png

TEC and O+ density [image: 0316_newfix] https://user-images.githubusercontent.com/22968399/65726449-58125200-e072-11e9-98d3-5e0bc85900ab.png

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AH5BFOA6L2K6QVFAE54UR63QLUSRRA5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7XAPUA#issuecomment-535693264, or mute the thread https://github.com/notifications/unsubscribe-auth/AH5BFOGXF3FZEUCMWU4BWETQLUSRRANCNFSM4IZRF3VQ .

[ZhuxiaoLi66]

Just a minor thing in the above code, IF(Z(JEQ).LE.ZLBDY+1.0) THEN should be IF(Z(JEQ).LE.(ZLBDY+1.0)) THEN ? Zhuxiao

[akubaryk]

Would be more readable like many potential changes to that source code but the result would be the same.

[naomimaruyama]

I was comparing RSLPSD.f between the master and development_ipe_refactor_merge branch. I found a major difference:

left: master right: branch

line 61--68 are commented out.

what this means is that the code does not do local equilibrium below 200km as stated in the comment line: "if the flux tube apex is below 200km, use local equilibrium".

The documentation for the solver RSLPSD.f can be found here

HOEQ is used to set densities at and below the lower boundary, before solving the densities.

[naomimaruyama]

@twfang ,

sw_lce=f was the setting in the old ipe: https://github.com/SWPC-IPE/WAM-IPE/blob/master/IPELIB/run/IPE.inp#L24

therefore, this part of the code was actually not used. I am sorry for the confusion.

[timfullerrowell]

Isn't there something similar with the 120 km reference calling HOEQ. Could check that as well.

On Fri, Sep 27, 2019 at 12:08 PM Naomi Maruyama notifications@github.com wrote:

@twfang https://github.com/twfang ,

sw_lce=f was the setting in the old ipe: https://github.com/SWPC-IPE/WAM-IPE/blob/master/IPELIB/run/IPE.inp#L24

therefore, this part of the code was actually not used. I am sorry for the confusion.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AH5BFOFWPAXUEPJEFXEEFI3QLZDY5A5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7ZVZNA#issuecomment-536042676, or mute the thread https://github.com/notifications/unsubscribe-auth/AH5BFOGCC5JFLKRXT2F55KTQLZDY5ANCNFSM4IZRF3VQ .

[naomimaruyama]

I keep working to compare FLIP routines between the master and development_ipe_refactor_merge branch.

CMINOR.f (chemistry routine to execute the chemical equilibrium below 200km):

ZNE is the electron density. It seems that this line could help the electron density going negative.

[akubaryk]

@gmillward do you remember what kind of values you were getting when you put that line in? Could do with a better solution if zne is significantly below 0.000001 -- would still be negative after that correction.

[naomimaruyama]

@gmillward , @akubaryk , I do not know what kind of values were expected because I never had problems with this part of the FLIP code.

[naomimaruyama]

I keep working further to compare FLIP routines between the master and development_ipe_refactor_merge branch.

RSPE2B.f (photoelectron routine)

IE is the index of the energy bin of the photoelectron energy, "E". This would get out of the energy loop quicker(IE starts from the IE max lowering down in the loop), resulting in changing the photoelectron spectra (that is used to heat plasma temperatures, and generate secondary ionization).

[twfang]

Did a test run with command out the if statement in CMINOR.f

! if(zne.lt.0.000001) zne = zne + 0.000001

Not seeing any changes in O+ density, velocity, and NO+ density within 6 hours of model run. We can try to revert some of these differences together for a test later on.

[naomimaruyama]

keep working further to compare FLIP routines

RSLPST.f

the difference of the temperature boundary condition between 100km vs 101km should not make much difference.

[naomimaruyama]

RSPRIM.f

The code appeared to look different, but both version turned out to use the same values, i.e. FNFAC=1.0

FNFAC is used to stabilize the O+ solution below 200 km.

[twfang]

Not sure if these altitude boundary conditions for density and temperature are necessary for WAM-IPE. They look the same in both version don't necesary mean these are correct. Someone should check with Phil.

[twfang]

Some updates from runs from Hera run. This run uses the time-varying drivers (i.e. timeobs) with empirical E field. All ion species are included in these TEC plots.

0316

0317

0318

0319

0320

[twfang]

With Adam's help, a run with the same version of WAM but with legacy IPE code (development) is done. This run uses the time-varying drivers (i.e. timeobs) with empirical E field. Essentially, same drivers and thermosphere for IPE, but use old IPE code. All ion species are included in these TEC plots.

0316

0317

0318

0319

0320

[twfang]

Conclusions from the previous two posts:

1. Refactored IPE appears to be a lot more stable compared to old IPE.
2. The yellow plots from old IPE are not associated with molecular ions. They are high values from H+.
3. Refactored IPE has problems on molecular ions from 03/20, which is not there in the old IPE

[timfullerrowell]

It is interesting that the first NaN for the molecular ions is at the start of a new day, or is it the last time-step of the previous day? Might be worth checking the driver file at 0UT on the 20th for odd values. When they first appear they are conjugate, but only above 120 km, which is where that boundary is. They are not conjugate below 120 km altitude, a clue there. Still need to find the source of the 120 km discontinuity.

On Fri, Oct 11, 2019 at 10:32 AM twfang notifications@github.com wrote:

Conclusions from the previous two posts:

1. Refactored IPE appears to be a lot more stable compared to old IPE.
2. The yellow plots from old IPE are not associated with molecular ions. They are high values from H+.
3. Refactored IPE has problems on molecular ions from 03/20, which is not there in the old IPE

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/SWPC-IPE/WAM-IPE/issues/348?email_source=notifications&email_token=AH5BFOH6FQ4QDSIBMPLQXJTQOCTB5A5CNFSM4IZRF3V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEBAQ6PY#issuecomment-541134655, or unsubscribe https://github.com/notifications/unsubscribe-auth/AH5BFOE4W6BVCB6RLC4CZGDQOCTB5ANCNFSM4IZRF3VQ .

[wanghj2015]

Why still 'yellow page'? Everywhere constant values?

[twfang]

Outputs from legacy IPE run can be found on Theia /scratch1/NCEPDEV/stmp2/Adam.Kubaryk/legacy_IPE_5day/netcdf