[mesa-users] controlling timestep using diffusion timescale
Alina Istrate
aistrate at astro.uni-bonn.de
Sat Mar 7 12:35:25 EST 2015
Hi Bill,
Thanks a lot for your answer.
I will play with the new suggestion and I'll get back when I have some
results .
Cheers,
Alina
On 03/06/2015 07:28 PM, Bill Paxton wrote:
> Hi Alina,
>
> On Mar 6, 2015, at 6:46 AM, Alina Istrate wrote:
>
>> Hi Bill,
>>
>>
>>> really?? can you provide me with inlist etc so I can see this in action.
>>> when I look at the code, I find subroutine do_solve_diffusion in star/private/diffusion.f90 doing this to save the number of steps taken with isolve:
>>> steps_used = iwork(17) ! number of accepted steps
>>> and the caller of do_solve_diffusion (in star/private/element_diffusion) does this
>>> s% num_diffusion_solver_steps = steps_used
>>> in star/private/timestep, the routine check_diffusion_steps_limit then uses s% num_diffusion_solver_steps to enforce the limits set by the controls diffusion_steps_limit and diffusion_steps_hard_limit.
>>>
>>> so if this isn't working, please give me details so we can fix it.
>> The value of diffusion_dt_div_timescale seems only to be used in the
>> do_step_loop routine, when using isolve its not taken into account at
>> all. The attached inlist shows this, modifying
>> diffusion_dt_div_timescale between 1d99 and 1d-99 makes no difference at
> first, that isn't addressing my question which was about your statement that timesteps limits based on diffusion steps used were not working when you selected isolve as the solver for diffusion. isn't that what you said in your previous email? that's what I thought you said, and that's what I was asking about and hoping you'd supply details.
>
> your reply isn't about the steps_used question, but instead is concerned about the effect of setting the diffusion_dt_div_timescale. so let's consider that.
>
> there are 2 different choices for solvers: the 1st order "step solver" and the higher order "isolve".
> the step solver uses diffusion_dt_div_timescale in its operation to pick sizes for substeps.
> isolve picks substep sizes in a much more sophisticated manner making use of embedded lower order solutions.
>
> so you are correct, isolve doesn't use the value of diffusion_dt_div_timescale -- it has no need for it.
>
>
>> all.
>>
>>>
>>>> Also, would it be possible to have the value computed by
>>>> get_timescale(...) in star/private/mod_diffusion.f90 accesible through a
>>>> public interface, and provide a timestep limit associated to diffusion
>>>> using this timescale (something like dt_div_diffusion_timescale_limit)?
>>> sure. but keep in mind that when you use isolve for this, it can take multiple substeps, each of which can be much smaller than the timescale calculated by this routine. this routine is for use in case you are not using isolve. and it is called at each substep taken by the diffusion solver, and the value changes a lot as the initial transients are smoothed out allowing larger substeps to be taken. so it isn't at all clear what you can usefully do with a particular value of this. we can of course save the final value of get_timescale in a star_data variable so you can inspect it from your run_star_extras.
>>>
>>> You can give that a try for yourself by making a patch. Just after the call on get_timescale, add this line
>>>
>>> s% xtra1 = timescale
>>>
>>> Then you can access s% xtra1 in your run_star_extras.
>>>
>>> If this turns out to be a useful way to set timesteps, let me know and we'll add something so you won't have to do a patch for future releases.
>>
>> The idea here is to use this diffusion timescale to reduce the
>> evolutionary step.
> which diffusion timescale? the timescale you get back from get_timescale before doing the solve will reflect small scale jumps in abundances introduced since the last step. as a result that starting timescale can be extremely small reflecting the need to remove transients. as the transients go away, later substeps can be much longer: the timescale can be orders of magnitude longer.
>
>
>
>> So, if I have the timescale stored in xtra1, I can
>> add to extras_finish step in run_star_extras.f : s% dt_next = min(s%
>> dt_next, s%xtra1*0.01), and that would give me the desired limit. I
>> modified mod_diffusion to store the computed timescale into xtra1,
>> however, since for isolve the timescale is not computed (could not find
>> more calls to get_timescale), I added this extra code to the routine
>>
>> call update_coeffs((kmax_rad_accel > 0 .and. steps_used >
>> 1), ierr)
>> call get_timescale( &
>> s, nz, nzlo, nzhi, m, nc, -s%
>> diffusion_min_X_hard_limit*0.5d0, v_advection_face, &
>> upwind_limit, X, X_face, C_div_X, SIG_face, GT_face,
>> AD_face, cell_dm, r_face, &
>> steps_used, total_num_iters, dbg, iter_dbg, j_dbg, k_dbg, &
>> j, k, class_chem_id, total_time, timescale)
>>
>> s% xtra1 = timescale
>>
>> I don't know exactly what update_coeffs does, so I might be breaking
>> something, but without it no timescale is computed (I essentially copied
>> this from do_step_loop). I attach run_star_extras.f and also
>> mod_diffusion.f ; the run_star_extras has the code to limit timestep
>> with xtra1, and mod_diffusion has the changes to store the timescale in
>> xtra1.
>>
>> This seems to work well, the computed timescale for the 1.2 solar mass
>> star changes more or less continuously, and demanding the timestep to be
>> around 1d-2 or 1d-3 of this value gives me the value for timestep I had
>> before when using max_years_for_timestep.
>>
> in my experience, making code changes you don't understand and hoping the results are okay doesn't work very well in the long term. ;D
>
> you have added calls to routines from a location in the code where they are not design to work and can be changing things that are not meant to be changed at the time and place from which you call them. will it work? maybe. maybe not. but please don't ask me to fix it!
>
> It might be a useful exercise for you to decide how you'd calculate a diffusion timescale yourself rather than pulling some number out of the bowels of the code. You have access to abundances, and other relevant information such as the following (from public/star_data.inc)
>
> ! element diffusion
> real(dp), pointer :: diffusion_D_self(:,:) ! self coeff for chem diffusion. [cm^2/sec]
> real(dp), pointer :: edv(:,:) ! velocity [cm/s] (species,nz) -- excludes chem diffusion.
> real(dp), pointer :: v_rad(:,:) ! velocity from radiative levitation [cm/s] (species,nz)
> real(dp), pointer :: g_rad(:,:) ! acceleration from radiative levitation [cm/s^2] (species,nz)
> real(dp), pointer :: typical_charge(:,:) ! (species,nz)
> real(dp), pointer :: diffusion_dX(:,:) ! (species,nz)
> real(dp), pointer :: E_field(:) ! (nz)
> real(dp), pointer :: g_field_element_diffusion(:) ! (nz)
>
> Is that enough for you to calculate a diffusion timescale? If not, what else do you need?
>
> Let us know what you learn.
>
> Cheers,
> Bill
>
>
>
>
>
>
>> Cheers,
>>
>>
>> Alina
>>
>>
>>> Cheers,
>>> Bill
>>>
>>>
>>>
>>>
>>>> Cheers!
>>>>
>>>> Alina
>>>>
>>>>
>>>> --
>>>> Alina Istrate
>>>> PhD student
>>>>
>>>>
>>>> Argelander Institute fur Astronomie
>>>> Auf dem Hugel 71
>>>> D-53121 Bonn
>>>> Germany
>>>> Room 2.024
>>>>
>>>>
>>>> ------------------------------------------------------------------------------
>>>> Dive into the World of Parallel Programming The Go Parallel Website, sponsored
>>>> by Intel and developed in partnership with Slashdot Media, is your hub for all
>>>> things parallel software development, from weekly thought leadership blogs to
>>>> news, videos, case studies, tutorials and more. Take a look and join the
>>>> conversation now. http://goparallel.sourceforge.net/
>>>> _______________________________________________
>>>> mesa-users mailing list
>>>> mesa-users at lists.sourceforge.net
>>>> https://lists.sourceforge.net/lists/listinfo/mesa-users
>> --
>> Alina Istrate
>> PhD student
>>
>>
>> Argelander Institute fur Astronomie
>> Auf dem Hugel 71
>> D-53121 Bonn
>> Germany
>> Room 2.024
>>
>> <inlist1.txt><mod_diffusion.f><run_star_extras.f>
--
Alina Istrate
PhD student
Argelander Institute fur Astronomie
Auf dem Hugel 71
D-53121 Bonn
Germany
Room 2.024
More information about the Mesa-users
mailing list