{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate mppnp flux information\n",
"\n",
"Warning: I am making the assumption of constant rho and temperature profiles for the whole simulation in the turning it on for a couple cycles. If that isn't true you will have to do this for a variety of cycles.\n",
"\n",
"\n",
"What we need to do this is the reaction rate information for all reactions in our network. To get that information from ```mppnp.f```, you need to make some edits to the existing code. \n",
"\n",
"#### Step 1: Locate where the subroutine ```printonecycle``` is defined in ```mppnp.f```\n",
"\n",
"#### Step 2. Edit ```printonecycle```\n",
"\n",
"This section inside ```printonecycle``` needs the additions in red to include the array of reaction rates ```v```\n",
"\n",
"\n",
" ! write out fluxes\n",
" if ((iplot_flux_type.eq.0).or.(iplot_flux_type.eq.2)) then\n",
" open(769,file='fluxes/flux_'//CMODELL//'_'//cprefix(1:\n",
" $ 5)//\".DAT\")\n",
" write(769,*)'Model = ',MODELL,', Zone = ',zone,\n",
" $ ', Age = ',dzeit\n",
" write(769,'(3x,a,2x,4(a,2x),6x,4(a,2x),2x,a,2(4x,a))')'#',\n",
" 1 'Z_k1','A_k1','Z_k3','A_k3','Z_k5','A_k5','Z_k7','A_k7',\n",
" 2 'flux [dY/dt]','energy [erg/(g*s)]','timescale reac [s]', 'v'\n",
" do i=1,nrnc1\n",
" write(769,'(i5,4(2x,i3),8x,4(2x,i3),3(8x,ES12.5))')i,\n",
" 1 int(znetw(k1(i))),int(anetw(k1(i))),\n",
" 2 int(znetw(k3(i))),int(anetw(k3(i))),\n",
" 3 int(znetw(k5(i))),int(anetw(k5(i))),\n",
" 4 int(znetw(k7(i))),int(anetw(k7(i))),\n",
" 5 max(1.0d-99,flux_to_print_0(i)),\n",
"! CR 6 max(1.0d-99,flux_to_print_0(i)*bind_energy_diff(i)),\n",
" 6 max(1.0d-99,flux_to_print_0(i))*bind_energy_diff(i),\n",
" 7 min(1.0d+90,time_scale_reac(i)),\n",
" 8 v(i)\n",
" end do\n",
" close(769)\n",
" end if\n",
"\n",
"\n",
"\n",
"#### Step 3. Allow ```printonecycle``` to execute\n",
"Also in the main code you will find the call to ```printonecycle```. You will want to make the edits I suggest so that it will actually call it.\n",
"\n",
"\n",
" !CR: flux plot implementation\n",
" !write out fluxes\n",
" call printonecycle(itag,icountmod,iolevel,dzeit,\n",
" $ T9,rho,y0,ZIS,considerisotope ,AN ,ZN,\n",
" $ isomeric_state,iplot_flux_option,v,\n",
" $ considerreaction,\n",
" $ reaction_rev_len,reaction_rev_num,reaction_rev_j,\n",
" $ iplot_flux_type,xm(itag),ipid)\n",
" if (iplot_flux_option.gt.0) then\n",
" ! every icountmodw,iplot_flux_tinterval cycles\n",
" if (modulo(icountmodw,iplot_flux_tinterval).eq.0) then\n",
" ! every iplot_flux_sinterval cycles\n",
" if (modulo(itag,iplot_flux_sinterval).eq.0) then\n",
" !write(*,*) 'print fluxfile, cycle: ',icountmodw,', zone: ',itag\n",
" call printonecycle(itag,icountmod,iolevel,dzeit,\n",
" $ T9,rho,y0,ZIS,considerisotope ,AN ,ZN,\n",
" $ isomeric_state,iplot_flux_option,v,\n",
" $ considerreaction,\n",
" $ reaction_rev_len,reaction_rev_num,reaction_rev_j,\n",
" $ iplot_flux_type,xm(itag),ipid)\n",
" !write(*,*) 'finished printonecycle'\n",
" end if\n",
" end if\n",
" end if\n",
" !CR\n",
"\n",
"\n",
"\n",
"#### Step 4. Get ready to run\n",
"Recompile ```mppnp.f```, create a directory called ```fluxes``` in the same directory as where you are running, and run with the ```ppn_frame.input``` parameter ```modstop = 2```. This will generate the flux files which you can then do what you want with them.\n",
"\n",
"Make sure that you undo these changes afterwards or have a separate version of ```mppnp.f``` with these edits. The flux files that are generated can be quite overbearing in terms of file size."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Ensure that the fluxes information is ```pandas``` friendly\n",
"\n",
"Since the way the mppnp flux plots are generated requires ```pandas```, it is important that you do this step since ```printonecycle``` doesn't format the files correctly (and I am unsure why at this moment). This step will basically just fix it so that the columns are good to go.\n",
"\n",
"#### Step 1. Create ```fluxes_fixed```\n",
"\n",
"In the same directory as your folder ```fluxes``` create another folder called ```fluxes_fixed```. This is where the fixed flux information from ```mppnp``` will be saved\n",
"\n",
"#### Step 2. Run the cells below\n",
"\n",
"Since this is parallelized, it will only take a handful of seconds."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"from pathlib import Path\n",
"from multiprocessing import Pool\n",
"from functools import partial"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def process_single_file(file, fixed_dir):\n",
" HEADER = '#\\tZ_k1\\tA_k1\\tZ_k3\\tA_k3\\tZ_k5\\tA_k5\\tZ_k7\\tA_k7\\tflux [dY/dt]\\tenergy [erg/(g*s)]\\ttimescale reac [s]'\n",
" \n",
" with open(file, mode='r') as f:\n",
" lines = f.readlines()\n",
" \n",
" fixed = [lines[0]] # Keep first line unchanged\n",
" con = ''\n",
" \n",
" for i, line in enumerate(lines):\n",
" if i == 0:\n",
" continue\n",
" \n",
" splits = [piece for piece in line.split(' ') if piece != '']\n",
" \n",
" if i == 1:\n",
" edit_line = HEADER\n",
" else:\n",
" if len(splits) > 1:\n",
" edit_line = '\\t'.join(splits)\n",
" else:\n",
" edit_line = splits[0]\n",
" \n",
" if i % 2 == 1:\n",
" con = edit_line\n",
" elif i % 2 == 0:\n",
" if i == 2:\n",
" con = con + '\\t' + str(edit_line)\n",
" else:\n",
" con = con[:-2] + '\\t' + str(edit_line)\n",
" fixed.append(con)\n",
" con = ''\n",
" \n",
" output_path = os.path.join(fixed_dir, os.path.basename(file))\n",
" with open(output_path, mode='w') as f:\n",
" f.writelines(fixed)\n",
"\n",
"def process_flux_files(flux_dir, fixed_dir, n_processes=None):\n",
" # Create output directory if it doesn't exist\n",
" os.makedirs(fixed_dir, exist_ok=True)\n",
" \n",
" # Get list of files to process\n",
" files = [os.path.join(flux_dir, f) for f in sorted(os.listdir(flux_dir)) \n",
" if not os.path.isdir(os.path.join(flux_dir, f))]\n",
" \n",
" # Create a partial function with fixed_dir\n",
" process_func = partial(process_single_file, fixed_dir=fixed_dir)\n",
" \n",
" # Use all available CPU cores if n_processes is None\n",
" with Pool(processes=n_processes) as pool:\n",
" # Process files in parallel\n",
" pool.map(process_func, files)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"flux_dir = '/scratch/f/fherwig/jissa/OC_RUNS/MLT_RUNS/flux_run/fluxes'\n",
"fixed_dir = '/scratch/f/fherwig/jissa/OC_RUNS/MLT_RUNS/flux_run/fluxes_fixed'\n",
"\n",
"process_flux_files(flux_dir, fixed_dir)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}