Abort criterion oplus

[rainbowsend]

The current implementation uses a fixed number of iterations to update O+ number concentration. This number is controlled by the nstep_sub input option. Instead, I propose to use nstep_sub as the maximum number of iterations and break the loop once a certain threshold/epsilon has been reached. This makes the model more versatile, as it automatically adapts the number of iterations to the current situation.

[hzfywhn]

This change looks good to me. It might help saving time during spin-up and other situations when O+ does not change drastically. What do you think? @phamkh @jvilaperez

[jvilaperez]

It looks good to me too. It can save some iterations in some cases.

[rainbowsend]

calc_terms is called only at the last iteration

tiegcm/src/oplus.F

Lines 105 to 106 in a0d5e33

	if (istep == nstep_sub)
	| call calc_terms(xnmbar,diffj,Fe,Fn,

Thus, when the loop is exited before the last iteration, it is never called.

It may be sufficient to move the call outside the loop. Since calc_terms only updates Fe and Fn, and neither variables occur as an argument in any call within the loop, this should be safe.

Sorry for overlooking this.

[hzfywhn]

calc_terms is called only at the last iteration

tiegcm/src/oplus.F

Lines 105 to 106 in a0d5e33

	if (istep == nstep_sub)
	| call calc_terms(xnmbar,diffj,Fe,Fn,

Thus, when the loop is exited before the last iteration, it is never called.
It may be sufficient to move the call outside the loop. Since calc_terms only updates Fe and Fn, and neither variables occur as an argument in any call within the loop, this should be safe.

Sorry for overlooking this.

That is right. We all missed this as well. Just post a quick fix to correct this.

[rainbowsend]

I conducted some tests for the runtime, I will also investigate the differences in O+ later.
For the test, I forwarded a 1 hour IMF run using 60 s integration time and 5° resolution.

The larger the maximal number of iterations, the more runtime is saved when using the abort criterion.

NSTEP_SUB	runtime (fixed, 006bb4e)	runtime (abort, 674d446)
200	49 s	46 s
400	79 s	75 s
600	118 s	102 s
800	148 s	117 s
1000	178 s	120 s

Using a less strict abort criterion could further reduce runtime. However, we need to ensure the differences in O+ are reasonable.

[hzfywhn]

I conducted some tests for the runtime, I will also investigate the differences in O+ later. For the test, I forwarded a 1 hour IMF run using 60 s integration time and 5° resolution.

The larger the maximal number of iterations, the more runtime is saved when using the abort criterion.

NSTEP_SUB runtime (fixed, 006bb4e) runtime (abort, 674d446)
200 49 s 46 s
400 79 s 75 s
600 118 s 102 s
800 148 s 117 s
1000 178 s 120 s
Using a less strict abort criterion could further reduce runtime. However, we need to ensure the differences in O+ are reasonable.

Thanks a lot! This benchmark is very helpful.

[rainbowsend]

I checked how the number of iterations influences O+. For this, I ran the same experiment, varying only the number of iterations (10, 50, 100, 200, 600, 1000, 2000) and not using the abort criterion. The following figure shows the difference in O+ relative to the run with 2000 iterations.
.

However, when using the aboart criterion, O+ differences increased with the number of iterations.
So I conducted an experiment in which the loop always exits just one iteration before the maximum number of iterations. With this setup, the differences to the run without an abort criterion are also unrealistically high (above 1E+6 cm-3 in the case of 1000 iterations. For context: Even when using only 10 iterations without abort, the differences in the figure above do not exceed 1000.)

Therefore, the abort criterion should not be used unless the cause of this behaviour is identified. I will check the code and search for dependencies of nsub.

[rainbowsend]

It turned out that nstep_sub is not only used to control the number of iterations but is also used in the calculation of O+:

• in the subroutine prep_oplus, which is called before the iteration
  

  tiegcm/src/oplus.F

  Lines 568 to 570 in 006bb4e

  	q_coeff(k,i,lat) =
  	| diffq(k,i,lat)+driftq(k,i,lat)+windq(k,i,lat)-
  	| dtx2inv*nstep_sub-op_loss(k,i,lat)

• in the subroutine smooth_oplus, which is called at the beginning of each iteration
  

  tiegcm/src/oplus.F

  Lines 674 to 678 in 006bb4e

  	optm1_smooth1(k,i,lat) = optm1(k,i,lat)-
  	| shapiro/nstep_sub*
  	| (optm1(k,i,lat+2)+optm1(k,i,lat-2)-
  	| 4.*(optm1(k,i,lat+1)+optm1(k,i,lat-1))+
  	| 6.*optm1(k,i,lat))

  
  

  tiegcm/src/oplus.F

  Lines 688 to 692 in 006bb4e

  	optm1_smooth(k,i,lat) = optm1_smooth1(k,i,lat)-
  	| shapiro/nstep_sub*
  	| (optm1_smooth1(k,i+2,lat)+optm1_smooth1(k,i-2,lat)-
  	| 4.*(optm1_smooth1(k,i+1,lat)+optm1_smooth1(k,i-1,lat))+
  	| 6.*optm1_smooth1(k,i,lat))

• in the subroutine iterate_oplus, which is called in each iteration
  

  tiegcm/src/oplus.F

  Lines 902 to 905 in 006bb4e

  	explicit(k,i,lat) =
  	| vdotn_h(k,i,lat)+bdotdh_bvel(k,i,lat)-
  	| diffexp(k,i,lat)-op_prod(k,i,lat)-
  	| optm1_smooth(k,i,lat)*dtx2inv*nstep_sub

• in the subroutine calc_terms, which is called in the last iteration
  

  tiegcm/src/oplus.F

  Lines 1162 to 1163 in 006bb4e

  	dopdt(k,i,lat) = dtx2inv*nstep_sub*
  	| (opout(k,i,lat)-optm1_smooth(k,i,lat))

This explains why the results are off when the loop is exited before nstep_sub is reached. To exit the loop, the calculations must be performed without depending on nstep_sub.