[mesa-users] error in angular momentum conservation

Bill Paxton paxton at kitp.ucsb.edu
Wed Dec 9 12:13:19 EST 2015

There is history column entry for 'total_angular_momentum' and another for 'log_total_angular_momentum'.

If there is no mass gain or loss, then things are of course simpler.

In that case errors can occur either during the evolution or during remeshing.  Check the remesh errors by turning off remeshing briefly (set okay_to_remesh to .false.).     If the problem persists, and you are feeling brave, go into the routine in star/private/solve_omega_mix.f90 and turn on some extra debugging code.  Here's a place to start (about line 400 or so):

            if (.not. (s% use_other_torque .or. s% use_other_torque_implicit .or. &
                  associated(s% binary_other_torque))) then
               ! check conservation for cases with no extra torque
               J_tot1 = total_angular_momentum(s) ! what we have

               !write(*,2) 'total_angular_momentum at end of do_solve_omega_mix', &
               !   s% model_number, J_tot1

               if (.false.) write(*,2) 'solve omega transport err', &
                  s% model_number, (J_tot0 - J_tot1)/J_tot0, J_tot0, J_tot1
               if (abs(J_tot0 - J_tot1) > 1d-6*J_tot0) then
                  s% retry_message = 'retry: failed to conserve angular momentum in mixing'
                  do_solve_omega_mix = retry
                  !stop 'do_solve_omega_mix'
               end if         
               if (dbg) then
                  write(*,2) 'final J_tot1', s% model_number, J_tot1
                  write(*,2) '(J_tot1 - J_tot0)/J_tot0', &
                     steps_used, (J_tot1 - J_tot0)/J_tot0, J_tot0, J_tot1
               end if
            end if

If you are adding or removing mass and want the angular momentum to change accordingly, things are trickier of course.  The 'accreted_material_j' test case checks angular momentum accretion -- you might look at the run_star_extras for that.   Here are the relevant lines:

         ! Check accretion of angular momentum
         if (s% model_number > 1 .and. s% mstar_dot > 0) then
            write(*,*) "Total accreted J should be:", s% accreted_material_j*s% mstar_dot*s% dt
            write(*,*) "Current J, old J, delta J:", s% total_angular_momentum, &
                s% total_angular_momentum_old, &
                s% total_angular_momentum - s% total_angular_momentum_old
            am_error = (s% accreted_material_j*s% mstar_dot*s% dt &
                - (s% total_angular_momentum - s% total_angular_momentum_old)) &
                / (s% accreted_material_j*s% mstar_dot*s% dt)
            write(*,*) "Relative diff in accreted J vs expected:", am_error
            if (abs(am_error) > 1d-5) then
                extras_check_model = terminate
                write(*,*) "Error in accreted J is too high!"
            end if
         end if

Please keep us informed -- we'll certainly want to fix any bugs you find.


On Dec 8, 2015, at 4:45 PM, Josiah Schwab wrote:

> Hi David,
>> In the history_columns.list, we can read
>> !error_in_energy_conservation
>>  !rel_error_in_energy_conservation
>> So I’m wandering if there is a (simple) way to track down relative
>> error in angular momentum conservation?
> No, I do not believe there is an an existing history column with that
> information.  You can add such a column yourself using run_star_extras.
> http://mesa.sourceforge.net/run_star_extras.html#toc-1-2
> Josiah
> ------------------------------------------------------------------------------
> _______________________________________________
> mesa-users mailing list
> mesa-users at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/mesa-users

More information about the Mesa-users mailing list