nrelax in dynamical sims

[Marios Chatzikos]

Hello,

I'm now looking at the coding of relaxation iterations in dynamical
sims.  As I mentioned previously, the check is written in several
different idioms, e.g.,

atom_leveln.cpp:381:                    iteration >
dynamics.n_initial_relax )
conv_init_solution.cpp:623:             if( iteration <=
dynamics.n_initial_relax )
conv_itercheck.cpp:326:                 if( iteration <=
dynamics.n_initial_relax+1 ||
dynamics.cpp:1054:              else if(iteration >
dynamics.n_initial_relax+1 )

and I'll like to rationalize them into a function, so as to keep policy
consistent and centralized.

Most of these calls can be refactored into a function, save for
conv_itercheck.cpp:326.  The comment in line 325 alludes to an e-mail
from Will Henney, but it is not clear if it refers to the iteration or
convergence portion of that if statement.  I would think that the +1 on
that line is not needed, but I'd like to confirm that before removing it.

The attached patch refactors all other like calls into a method, and
should be pretty straight-forward.  Notice that the call on
dynamics.cpp:1065 is correct in the context of the if statement (the
case iteration == dynamics.n_initial_relax+1 on that line is never
executed).

I've run the test suite with it, and got the following performances botches:

dynamics_orion_flow.out:  ChkMonitor botch>>itrz                      
0    27.3764 = 24.1900  -0.132   0.050
dynamics_orion_flow.out:  ChkMonitor botch>>zone                      
0   356.0000 = 419.0000   0.150   0.050             <<BIG (3sig) BOTCH!!

Let me know if you have any comments.

Marios

[Marios Chatzikos]

Please disregard the previous patch.  There was an error that caused the
reported botches.  The attached one passes clean.

Thanks,

Marios

[Marios Chatzikos]

Regarding conv_itercheck.cpp:326, note that removing the +1 triggers no
botches in the test suite.  If you could weigh in on its presence,
that'd be great.

Thanks,

Marios

[Robin Williams]

Hi Marios --

2002 is past the particle horizon of my e-mail.  I have a vague feeling the logic was intentional for steady-state structures, to ensure that at least one iteration was run with dynamical sources rather than having the code declare convergence & stop before it had actually switched all the physics on.

I think "initial" in "initial_relax" is important -- I'm not sure of a better name, though.  dynamics.inPrewash() certainly isn't one, though it kind of captures the spirit of the thing ;-)

The code in DynaIterEnd is about saving the values from the iteration just completed so it's ready to act as a source for the next iteration.  I *think* the comment in iter_startend.cpp just before IterRestart() suggests that the logic can be consistent with elsewhere because "iteration" has already been incremented by here.  Note also the additional "iteration == dynamics.n_initial_relax+1" tests in dynamics.cpp, one of which is followed by a somewhat detailed comment about the logic.  It might be less liable to fencepost errors if this whole thing was turned into a logical flag (toggle) [or state variable], which was flipped as required, rather than doing the repeated tests on "iteration".

  Robin

--
--
http://groups.google.com/group/cloudy-dev
---
You received this message because you are subscribed to the Google Groups "cloudy-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cloudy-dev+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/4cf75fe5-dc69-41ca-13b2-35ead96ca24c%40gmail.com.

[Robin Williams]

...of course, it's possible that using a single iteration loop to converge both the optical depth scale (for interation < n_initial_relax) and integrate forward the flow/time dependence (thereafter) is over-ambitious -- it does rather assume that the optical depths can look after themselves once the dynamics has been switched on.  The logic might be clearer & more robust if this were two nested convergence loops.

[Marios Chatzikos]

Hi Robin,

Thanks for the response.

Okay, I think I will reinstate the "+1" in conv_itercheck.cpp.  The orion sim worked just fine, but more extensive testing is needed to make sure this is correct.

I could set up a state variable to do the job without the repeated checks on the iteration value, but, of course, the variable would change after such a check against the iteration number.

Marios

To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAJBtPxce31jgjKr6i9sZRKqSEheev13OOMXfVNJnFpvF-BoM4w%40mail.gmail.com.

[Robin Williams]

The point with the state variable is that you need to think through the various behaviours you want, and switch the variable exactly when you want them to change.

The atomic elements of the behaviour seem to be A) iterating the optical depths to convergence; B) copying the old-state mesh data as reference for the dynamics; C) running with dynamical sources enabled; D) declaring convergence.

Obviously C) depends on B) having happened at the end of the last step, so B) needs to be enabled at the end of the last iteration of A).  D) requires at least one "dynamics" step to have run.

I think the logic can work out if you have a single toggle between burn-in and dynamics, but it's a little intricate:

burnin = true

dynamics = true

while (true) {

  if (dynamics && ! burnin)

    step(with dynamics sources)

  else

    step(no sources)

  if (converged?()) {

    if (dynamics && burnin)

      burnin = false

    else

      break // Declare convergence

  }

  if (dynamics && ! burnin)

    copystate()

}

  Robin

To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CA%2B2Qw8oHSRGTHuun%2B9Xiny6LYcqBCVQOBSWFHrFbPoE%3D2d_aCw%40mail.gmail.com.

[Marios Chatzikos]

Hi Robin,

I have changed the method name to isInitialRelaxIteration(), which hopefully clarifies its intent.  Your other suggestion seems useful, but I feel like it'd be more appropriate to consider when I start working on RT.

I'll push branch nrelax to nublado, and create a merge request.

Thanks,

Marios

To view this discussion on the web visit https://groups.google.com/d/msgid/cloudy-dev/CAJBtPxfC7y5cohyY5NvDnHstFsTy%2BTRbLDJAp-n6gWUPozC2aA%40mail.gmail.com.