problem in running my simulation of ISM with adding Radiation
Hesam
By the way, it have stopped me for some weeks ! I don't know what is the problem in spite of searching and looking at other samples such as solar atmosphere .
Frederick Gent
--
You received this message because you are subscribed to the Google Groups "pencil-code-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pencil-code-dis...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/pencil-code-discuss/fab2b34e-e850-428c-b395-5c19e3cecbbdn%40googlegroups.com.
Frederick Gent
Axel Brandenburg
The recommended way to share a run is to do: pc_tar_run
You do this while being in the run directory, and then it produces a
gzipped tar file in that same directory with a corresponding name.
Use that for sending.
Cheers, Axel
On Sat, May 17, 2025 at 05:57:09AM +0000, Frederick Gent wrote:
> Hi Hasam,
>
> Try sharing the files rather than screen shots and I'll be better able to help. (*.in and *.local)
>
> Cheers,
>
> Fred
>
>
> Hello,
> I am doing a simulation of interaction between matter and electromagnetic field in the interstellar medium and as interstellar sample in the pencil code, i considered a box with x and y periodic boundary condition and bcz=s,s,ouf,ouf,ism,ouf,ouf,s,ouf and an initial SN I explosion in the center of the box. It somehow works and i have set several parameters to have a good result and now i want to add radiation but it does not work and it gives "fatal error denergy_dt" and when i turn off SN, it works !
> By the way, it have stopped me for some weeks ! I don't know what is the problem in spite of searching and looking at other samples such as solar atmosphere .
> cparam:
> Make:
> [Make.png]
> start:
> [start.png]
> run:
> [run.png]
> --
> You received this message because you are subscribed to the Google Groups "pencil-code-discuss" group.
> To view this discussion visit https://groups.google.com/d/msgid/pencil-code-discuss/fab2b34e-e850-428c-b395-5c19e3cecbbdn%40googlegroups.com<https://groups.google.com/d/msgid/pencil-code-discuss/fab2b34e-e850-428c-b395-5c19e3cecbbdn%40googlegroups.com?utm_medium=email&utm_source=footer>.
> --
> You received this message because you are subscribed to the Google Groups "pencil-code-discuss" group.
> To view this discussion visit https://groups.google.com/d/msgid/pencil-code-discuss/CAM%2B36SuGoYFneNNrOeNBTesNcAPQc9cAYw4Yjet3O_%2Bk5eRvkw%40mail.gmail.com<https://groups.google.com/d/msgid/pencil-code-discuss/CAM%2B36SuGoYFneNNrOeNBTesNcAPQc9cAYw4Yjet3O_%2Bk5eRvkw%40mail.gmail.com?utm_medium=email&utm_source=footer>.
حسام نادری
Frederick Gent
Dear Hesam.
Having inspected the parameters, the grid resolution in each direction is 0.002, which necessitates a hyper_diffusion coefficient 2^5 times smaller than that used, which is derived for a grid spacing of 0.004. It is possible that this is part of the problem, but I identify a few other issues. It is unlikely that you can resolve the shock fronts without applying shock diffusion to density and to entropy (not magnetic). It is also unlikely that the system will resolve the grid scales without also applying hyper diffusion to each of entropy, viscosity and magnetic (not density). I have no experience yet working with ionization.
I do not know your objectives, but viscosity nu and resistivity eta seem unnecessarily high for the resolution.
It may also be that the default values for the SN injection do not suit your purposes or resolution, so I would consider specifying an explicit width_SN in kpc (code units) and if it is only 5 parsecs as likely with rfactor_SN=2.5 and dx=2 parsecs. I would check maxTT in the file data/sn_series.dat. If it is above 10⁸ K it is likely too hot and consider increasing the radius to bring it down.
I notice you are using the default timestep.f90, but the parameters are for timestep_rkf_lowsto.f90. I would recomment adding to src/Makefile.local TIMESTEP = timestep_rkf_lowsto and using itorder=4 in run.in. In any case dtmin=1e-11 rather than 1e-9 and if you stick with the Courant timestep module I would reduce the Courant numbers.
I cannot tell if you only want 1 initial SN or many. If the
former you need to set lSNII=F and lSNI=F and
l_persist_overwrite_lSNII=T and l_persist_overwrite_lSNI=T. This
is because of the order in which persistent variables are read
from the snapshot, it will overwrite values in run.in unless it is
disabled with the l_persist* switches.
What is the source of the init_ism file? Has it been derived using a convergent hydrostatic equilibrium to the 1D model for this resolution with same xyz0(3) and xyz1(3)? Otherwise it may be out unstable.
I have made recent adjustments, so would recommend pulling the
latest version of the code if not already done so.
Hopefully this helps and keep in touch if still struggling.
Cheers,
Fred
To view this discussion visit https://groups.google.com/d/msgid/pencil-code-discuss/CAENm1Dzm3nei_YDzp_%2BJ_7WWuC_VjSVdiJ3D71N7v_srCC3B1Q%40mail.gmail.com.
حسام نادری
To view this discussion visit https://groups.google.com/d/msgid/pencil-code-discuss/f496a179-7946-41c9-b8d9-7a273659d38c%40gmail.com.
Frederick Gent
Hi Hesam,
You are working at the frontier of the code with resepect to this combination of physics, so you inevitably face some teething problems, but not unsurmountable.
SN turbulence is fairly robust now with the ideal gas equation of state, but you now are introducing radiation and ionization.
While ionization in principle should be compatible with entropy, Axel pointed out to me last week that due to the complex relations between temperature and ionization it would be preferable to work with temperature_ionization directly instead of entropy. This would require a check that any functionality you are applying to the entropy module is correctly implemented in the temperature module, such as shock diffusion, hyper diffusion, etc. Given you were working successfully without SN, this maybe is not the root of the problem now, and can be deferred.
When referring to the values in init_ism.in this initial condition is a themo-hydrostatic equilibrium solved in 1D (pencil-code/samples/supernova-driven-turbulence/1Dhd_init_equidistant) for a given resolution and Lxzy(3)/heating/cooling/gravity combination. If it is used for a different profile in z or different heating and cooling function it will be unstable and if the perturbations are strong enough could lead to a crash. Replacing ideal gas with ionization eos might also alter the thermo-hydrostatic equilibrium. An easy test of the equilibrium would be to use your profile without SN and run a 1D simulation in z, to check whether the system is static. you might get some motion, but if the velocities are small (umax<<1) this would not explain the crash and your. If you do need to solve the 1D profile again for your setup I would recommend using nu and chi of order unity and then reducing both to the values planned for the 3D problem as the system approaches equilibrium. The is a script in pencil-code/samples/supernova-driven-turbulence/1Dhd_init_equidistant ism_save.py, which can be used while you are running the 1D simulation to monitor the approach to equilibrium and to create the file init_ism.in to add to your 3D simulation setup. This of course may not be the problem causing your crash, if the system works adequately without the SN.
I am suspicious that the values for lnrho_* in
data/pc_constant.pro appear high given the unit of density is
1.6737236e-24 g/cmcube. I would expect values for lnrho_H around
0.
unit_magnetic is too high, yielding mu0=1.6002754082279367E-014. This will require the code to reconcile huge numbers in the induction equation, which even with double precision is risky. Use lfix_unit_std=T in init_pars and this will ensure mu0=1 and unit_magnetic=4.4843788716874900E-007. It will also alter the unit of temperature from 1 to something higher, which might report smaller coefficients in pc_constants.pro, but I have no expereince yet with ionization. It was some surprise to me a few years ago when I discovered the sensitivity of the equations to having coefficients far from unity.
While it should not matter in the initial iteration the vertical boundaries will work fine with bcz=9*'s'.
The most likely cause of the crash is too hot an initial explosion. There are files data/sn_series.dat and data/sn_history.dat which record what is tried and applied in the explosion and I do not see these. While it should not be critical, I would try to ensure the explosion site lies on a grid point. Typically z=0.0 is between the nearest two points at the midplane due to the symmetry of the disc profile. So, set center_SN_z=\pm 0.0009765625, however in the periodic directions 0.0 does lie on the grid without lshift_origin=T.
initss="nothing" in entropy_init_pars is compatible with initial_conditions/ths_equilibrium_ism as is initaa='nothing' in magnetic_init_pars. initss="Ferriere" likely doubles up on the energy or overwrites the equilbrium in init_ism.in.
From your screen shot the time_step appears to be about 5e-8.
This is likely too high. However, I would try the RFK timestep
again, timestep_rkf_lowsto.f90, but set
timestep_scaling=8*"rel_err" not "cons_frac_err" and reduce
eps_rkf=0.001 and dt_ratio=1e-6. If you still use the Courant
timestep cdt=0.4 is way too high. Try cdt=0.01, cdtv=0.01 and
cdts=0.01. The timestep will drop after an explosion, but should
recover in between. If the code still crashes with either timestep
module try setting dt=1e-9 in run.in or lower and see if the code
runs.
In your tar files you still are not including hyper or shock diffusion, so I would expect the code to crash.
If you still have problems getting the code not to crash, we might meet over zoom and go through the steps.
Best regards,
Fred
Dear Fred,Thank you for your help , i have done all changes as you suggested , first i got an error for the new timstep that you said. So i switched to the previous one but with smaller dtc . Then it ran but after 1 step , again denergy_dt error. I turned off SN and again there was error and i turned off shear , there was no error and it ran ! I add SN in this case there was again error ( attached screenshot) and I am really confused.
Summary of my work: i have a box with an initial SN in the center of box and i want to see how it affects the ism in my box and how long does it take to dissolve. I had some simulations without radiation that the box got hotter and hotter as time goes on but now i want to turn on radiation in the box and i am unfortunately stuck ! Also i want to test a suitable boundary condition for this work.... .
Also i didn't understand your question about init_ism file , in fact I used the interstellar sample in the code and its files and worked on it for some months.I really appreciate your time and would appreciate your help in fixing this problem.
With best regards,Hesam