[mesa-users] Show the evolution of the model at each step of the solver iteration

[Ehsan Moravveji]

Dear Mathieu,
This is a very interesting question. I did not have a look at such yet, but maybe the following ideas would help.
To trace the convergence, you can add the following public variables to your profiles which store the residuals for each of the structure+hydro equations being solved:

      ! index definitions for the equations (= 0 if equation not in use)
         integer :: equP ! momentum conservation
         integer :: equT ! energy transport
         integer :: equR ! mass-volume-density relation
         integer :: equL ! energy conservation
         integer :: equv ! velocity = radius time derivative

To dump such columns into your profile, just add them to your how_many_extra_profile_columns and data_for_extra_profile_columns and call the public routine star_write_profile_info(). 
But, the question - as you point out - is from which other hook to call this.
Well, I am not a hook guru, but just "propose" the following:

Since v. 6871, Bill has implemented a new hook called s% other_after_set_mixing_info.
You may do the following to activate it:

      subroutine extras_controls(s, ierr)
         type (star_info), pointer :: s
         integer, intent(out) :: ierr
         ierr = 0
         s% other_after_set_mixing_info => check_convergence
      end subroutine extras_controls


      subroutine check_convergence(id, ierr)
         use const_def, only: dp
         use star_def
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type(star_info), pointer :: s
         integer :: k
         ierr = 0

         call star_ptr(id, s, ierr)
         if (ierr/=0) then
            write(*,*) 'Error: run_star_extras: check_convergence: star_ptr failed'
            return
         endif

         print*, 'M = ', s% model_number

      end subroutine check_convergence

Once you manage to compile and run it, the print statement shows that the Newton solver is being called 3 to 4 times, until the next top up in the s% model_number.
This "might be" where you would call your star_write_profile_info().

I hope this helps, but there are several gurus who use other hooks, and I am sure they can assist you.

Keep me posted please on how it goes (cause I'm super interested in this stuff).

Cheers
Ehsan.



On Sep 4, 2014, at 10:49 AM, Mathieu wrote:

> Hi everyone,
> 
> I was wondering if there is a way to look at the structure of a model at 
> each iteration of the newton solver, before the solution for the next 
> timestep is found.
> I'd like to make some plots of the structure before the convergence to a 
> solution, e.g. T(\rho), M(R), etc...
> The purpose of this is to show how the numerical procedure "shakes" the 
> star at every timestep to get it to the hydrostatic solution, in order 
> to illustrate how the numerics of MESA work.
> Of course, I would need it only for a few timesteps, just to have a few 
> examples.
> So, is there a way to save something like a profile*.data for each 
> iteration of the newton solver? In "evolve.f" the profile is saved at 
> the end of the step, so I guess I will not be able to get a true profile 
> for each iteration. I think what I need is what is stored in the 
> structure pointed by x in the routine newton in 
> "$MESA_DIR/star/private/star_newton.f":
> 
>     subroutine newton( &
>          s, nz, nvar, x, xold, sparse, &
>          lrd, rpar_decsol, lid, ipar_decsol, which_decsol, &
>          tol_correction_norm, &
>          xscale, equ, work, lwork, iwork, liwork, AF, &
>          lrpar, rpar, lipar, ipar, convergence_failure, ierr)
> 
>          type (star_info), pointer :: s
>          ! the primary variables
>          integer, intent(in) :: nz ! number of zones
>          integer, intent(in) :: nvar ! number of variables per zone
>          ! the total number of primary variables is neq
>          real(dp), pointer, dimension(:) :: x ! =(nvar,nz)
>          ! new vector of primaries
>          real(dp), pointer, dimension(:) :: xold ! =(nvar,nz)
>          ! old vector of primaries
> 
> So maybe I could use the routine newton_core_dump defined in star_newton.f .
> Am I correct to try to get the values pointed by x if I want to get the 
> variables for each zone at each iteration of the solver? Is there a 
> function in the run_star_extras which is called at each iteration of the 
> solver that I could use?
> 
> I appreciate any hint!
> Thanks,
> 
> Mathieu
> 
> 
> ------------------------------------------------------------------------------
> Slashdot TV.  
> Video for Nerds.  Stuff that matters.
> http://tv.slashdot.org/
> _______________________________________________
> mesa-users mailing list
> mesa-users at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/mesa-users