diff --git a/examples/Feedback/feedback.pro b/examples/Feedback/feedback.pro new file mode 100644 index 0000000000000000000000000000000000000000..93305120e00d54a377bf0a0ee75bcc2fe7586ab0 --- /dev/null +++ b/examples/Feedback/feedback.pro @@ -0,0 +1,22 @@ +base = 'Feedback' +inf = 'Feedback_003.hdf5' + +blast = [0.5,0.5,0.5] ; location of blast +pos = h5rd(inf,'PartType0/Coordinates') +vel = h5rd(inf,'PartType0/Velocities') +rho = h5rd(inf,'PartType0/Density') +utherm = h5rd(inf,'PartType0/InternalEnergy') + +; shift to centre +for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic] + +;; distnace from centre +dist = fltarr(n_elements(rho)) +for ic=0,2 do dist = dist + pos[ic,*]^2 +dist = sqrt(dist) + +; radial velocity +vr = fltarr(n_elements(rho)) +for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*] +vr = vr / dist +end diff --git a/examples/Feedback/feedback.yml b/examples/Feedback/feedback.yml new file mode 100644 index 0000000000000000000000000000000000000000..74d2445973db9be6a955efe9be5be86705c25bc1 --- /dev/null +++ b/examples/Feedback/feedback.yml @@ -0,0 +1,44 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1 # Grams + UnitLength_in_cgs: 1 # Centimeters + UnitVelocity_in_cgs: 1 # Centimeters per second + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 5e-2 # The end time of the simulation (in internal units). + dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: Feedback # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-3 # Time between statistics output + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours. + max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./Feedback.hdf5 # The file to read + +# Parameters for feedback + +SN: + time: 0.001 + energy: 1.0 + x: 0.5 + y: 0.5 + z: 0.5 diff --git a/examples/Feedback/getGlass.sh b/examples/Feedback/getGlass.sh new file mode 100755 index 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c --- /dev/null +++ b/examples/Feedback/getGlass.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5 diff --git a/examples/Feedback/log b/examples/Feedback/log new file mode 100644 index 0000000000000000000000000000000000000000..94225c1110423b6f2b4ba6fed47f33ec4743a4b6 --- /dev/null +++ b/examples/Feedback/log @@ -0,0 +1,468 @@ + Welcome to the cosmological hydrodynamical code + ______ _________________ + / ___/ | / / _/ ___/_ __/ + \__ \| | /| / // // /_ / / + ___/ /| |/ |/ // // __/ / / + /____/ |__/|__/___/_/ /_/ + SPH With Inter-dependent Fine-grained Tasking + + Version : 0.3.0 + Revision: v0.3.0-639-g352f3fa0-dirty, Branch: TTFeedback + Webpage : www.swiftsim.com + + Compiler: ICC, Version: 16.0.20160415 + HDF5 library version: 1.8.14 + +[00000.0] main: CPU frequency used for tick conversion: 1994987486 Hz +[00000.0] main: sizeof(struct part) is 128 bytes. +[00000.0] main: sizeof(struct xpart) is 32 bytes. +[00000.0] main: sizeof(struct gpart) is 96 bytes. +[00000.0] main: sizeof(struct task) is 96 bytes. +[00000.0] main: sizeof(struct cell) is 640 bytes. +[00000.0] main: Reading runtime parameters from file 'feedback.yml' +[00000.0] source_terms_print: Single SNe of energy= 1.000000e+00 will explode at time= 1.000000e-03 at location (5.000000e-01,5.000000e-01,5.000000e-01) +[00000.0] main: Reading ICs from file './Feedback.hdf5' +[00000.0] read_ic_single: IC and internal units match. No conversion needed. +[00000.0] main: Reading initial conditions took 22.088 ms. +[00000.0] main: Read 32768 gas particles and 0 gparts from the ICs. +[00000.0] main: space_init took 2.955 ms. +[00000.0] main: space dimensions are [ 1.000 1.000 1.000 ]. +[00000.0] main: space is periodic. +[00000.0] main: highest-level cell dimensions are [ 10 10 10 ]. +[00000.0] main: 32768 parts in 1000 cells. +[00000.0] main: 0 gparts in 1000 cells. +[00000.0] main: maximum depth is 0. +[00000.0] engine_init: no processor affinity used +[00000.1] engine_policy: engine policies are [ steal keep numa_affinity hydro sourceterms ] +[00000.1] hydro_props_print: Adiabatic index gamma: 1.666667. +[00000.1] hydro_props_print: Hydrodynamic scheme: Gadget-2 version of SPH (Springel 2005) in 3D. +[00000.1] hydro_props_print: Hydrodynamic kernel: Cubic spline (M4) with 48.00 +/- 0.10 neighbours (eta=1.234800). +[00000.1] hydro_props_print: Hydrodynamic integration: CFL parameter: 0.1000. +[00000.1] hydro_props_print: Hydrodynamic integration: Max change of volume: 2.00 (max|dlog(h)/dt|=0.231049). +[00000.1] engine_init: Absolute minimal timestep size: 1.862645e-10 +[00000.1] engine_init: Minimal timestep size (on time-line): 9.536743e-08 +[00000.1] engine_init: Maximal timestep size (on time-line): 9.765625e-05 +[00000.1] main: engine_init took 6.579 ms. +[00000.1] main: Running on 32768 gas particles and 0 DM particles from t=0.000e+00 until t=5.000e-02 with 1 threads and 1 queues (dt_min=1.000e-07, dt_max=1.000e-04)... +[00000.1] engine_init_particles: Running initialisation fake time-step. +# Step Time Time-step Updates g-Updates Wall-clock time [ms] + 0 0.000000e+00 0.000000e+00 0 0 1048.558 + 1 9.765625e-05 9.765625e-05 32768 0 955.300 + 2 1.953125e-04 9.765625e-05 32768 0 1260.160 + 3 2.929688e-04 9.765625e-05 32768 0 1261.074 + 4 3.906250e-04 9.765625e-05 32768 0 1240.943 + 5 4.882812e-04 9.765625e-05 32768 0 1204.562 + 6 5.859375e-04 9.765625e-05 32768 0 1173.909 + 7 6.835938e-04 9.765625e-05 32768 0 1146.211 + 8 7.812500e-04 9.765625e-05 32768 0 1130.690 + 9 8.789063e-04 9.765625e-05 32768 0 1136.483 + 10 9.765625e-04 9.765625e-05 32768 0 1113.689 +[00013.9] do_supernova_feedback: u_old= 1.500000e-06 u_new= 3.276800e+04 S_old= 9.973413e-07 S_new= 2.178725e+04 rho= 1.004002e+00 +[00013.9] do_supernova_feedback: p->s 2.178725e+04 + 11 1.074219e-03 9.765625e-05 32768 0 1122.928 + 12 1.123047e-03 4.882813e-05 32768 0 1166.416 + 13 1.171875e-03 4.882813e-05 57 0 26.710 + 14 1.220703e-03 4.882813e-05 32768 0 1099.410 + 15 1.269531e-03 4.882813e-05 57 0 27.786 + 16 1.318359e-03 4.882813e-05 32768 0 1092.165 + 17 1.367188e-03 4.882813e-05 57 0 29.624 + 18 1.416016e-03 4.882813e-05 32768 0 1083.173 + 19 1.464844e-03 4.882813e-05 81 0 28.283 + 20 1.513672e-03 4.882813e-05 32768 0 1080.139 + 21 1.562500e-03 4.882813e-05 81 0 89.744 + 22 1.611328e-03 4.882813e-05 32768 0 1119.688 + 23 1.660156e-03 4.882813e-05 80 0 30.407 + 24 1.708984e-03 4.882813e-05 32768 0 1066.095 + 25 1.757813e-03 4.882813e-05 80 0 28.911 + 26 1.806641e-03 4.882813e-05 32768 0 1174.870 + 27 1.855469e-03 4.882813e-05 92 0 34.274 + 28 1.904297e-03 4.882813e-05 32768 0 1063.577 + 29 1.953125e-03 4.882813e-05 92 0 30.317 + 30 2.001953e-03 4.882813e-05 32768 0 1063.683 + 31 2.050781e-03 4.882813e-05 108 0 26.311 + 32 2.099609e-03 4.882813e-05 32768 0 1188.527 + 33 2.148438e-03 4.882813e-05 114 0 27.428 + 34 2.197266e-03 4.882813e-05 32768 0 1080.263 + 35 2.246094e-03 4.882813e-05 42 0 19.459 + 36 2.294922e-03 4.882813e-05 32768 0 1074.380 + 37 2.343750e-03 4.882813e-05 42 0 18.298 + 38 2.392578e-03 4.882813e-05 32768 0 1165.363 + 39 2.441406e-03 4.882813e-05 42 0 21.901 + 40 2.490234e-03 4.882813e-05 32768 0 1079.358 + 41 2.539063e-03 4.882813e-05 42 0 17.172 + 42 2.587891e-03 4.882813e-05 32768 0 1073.910 + 43 2.636719e-03 4.882813e-05 42 0 15.637 + 44 2.685547e-03 4.882813e-05 32768 0 1068.533 + 45 2.734375e-03 4.882813e-05 12 0 42.965 + 46 2.783203e-03 4.882813e-05 32768 0 1123.024 + 47 2.832031e-03 4.882813e-05 12 0 15.465 + 48 2.880859e-03 4.882813e-05 32768 0 1066.808 + 49 2.929688e-03 4.882813e-05 12 0 13.545 + 50 2.978516e-03 4.882813e-05 32768 0 1058.231 + 51 3.027344e-03 4.882813e-05 12 0 13.778 + 52 3.076172e-03 4.882813e-05 32768 0 1058.188 + 53 3.125000e-03 4.882813e-05 12 0 41.869 + 54 3.173828e-03 4.882813e-05 32768 0 1116.429 + 55 3.222656e-03 4.882813e-05 12 0 15.152 + 56 3.271484e-03 4.882813e-05 32768 0 1065.698 + 57 3.320313e-03 4.882813e-05 12 0 13.500 + 58 3.369141e-03 4.882813e-05 32768 0 1054.867 + 59 3.417969e-03 4.882813e-05 12 0 13.381 + 60 3.466797e-03 4.882813e-05 32768 0 1052.531 + 61 3.515625e-03 4.882813e-05 6 0 13.089 + 62 3.564453e-03 4.882813e-05 32768 0 1232.052 + 63 3.613281e-03 4.882813e-05 6 0 15.675 + 64 3.662109e-03 4.882813e-05 32768 0 1145.922 + 65 3.710938e-03 4.882813e-05 6 0 14.843 + 66 3.759766e-03 4.882813e-05 32768 0 1144.235 + 67 3.808594e-03 4.882813e-05 6 0 12.273 + 68 3.857422e-03 4.882813e-05 32768 0 1160.155 + 69 3.906250e-03 4.882813e-05 6 0 13.384 + 70 3.955078e-03 4.882813e-05 32768 0 1201.139 + 71 4.003906e-03 4.882813e-05 6 0 15.017 + 72 4.052734e-03 4.882813e-05 32768 0 1342.790 + 73 4.101563e-03 4.882813e-05 6 0 49.967 + 74 4.199219e-03 9.765625e-05 32768 0 1265.088 + 75 4.296875e-03 9.765625e-05 32768 0 1202.305 + 76 4.394531e-03 9.765625e-05 32768 0 1143.659 + 77 4.492188e-03 9.765625e-05 32768 0 1169.572 + 78 4.589844e-03 9.765625e-05 32768 0 1162.901 + 79 4.687500e-03 9.765625e-05 32768 0 1147.462 + 80 4.785156e-03 9.765625e-05 32768 0 1374.074 + 81 4.882812e-03 9.765625e-05 32768 0 1276.302 + 82 4.980469e-03 9.765625e-05 32768 0 1266.279 + 83 5.078125e-03 9.765625e-05 32768 0 1266.334 + 84 5.175781e-03 9.765625e-05 32768 0 1270.059 + 85 5.273438e-03 9.765625e-05 32768 0 1264.018 + 86 5.371094e-03 9.765625e-05 32768 0 1255.714 + 87 5.468750e-03 9.765625e-05 32768 0 1329.432 + 88 5.566406e-03 9.765625e-05 32768 0 1240.819 + 89 5.664063e-03 9.765625e-05 32768 0 1236.233 + 90 5.761719e-03 9.765625e-05 32768 0 1234.927 + 91 5.859375e-03 9.765625e-05 32768 0 1263.383 + 92 5.957031e-03 9.765625e-05 32768 0 1232.714 + 93 6.054688e-03 9.765625e-05 32768 0 1236.793 + 94 6.152344e-03 9.765625e-05 32768 0 1239.297 + 95 6.250000e-03 9.765625e-05 32768 0 1314.640 + 96 6.347656e-03 9.765625e-05 32768 0 1240.433 + 97 6.445313e-03 9.765625e-05 32768 0 1259.359 + 98 6.542969e-03 9.765625e-05 32768 0 1292.100 + 99 6.640625e-03 9.765625e-05 32768 0 1267.954 + 100 6.738281e-03 9.765625e-05 32768 0 1256.223 + 101 6.835938e-03 9.765625e-05 32768 0 1291.486 + 102 6.933594e-03 9.765625e-05 32768 0 1263.898 + 103 7.031250e-03 9.765625e-05 32768 0 1290.375 + 104 7.128906e-03 9.765625e-05 32768 0 1362.429 + 105 7.226563e-03 9.765625e-05 32768 0 1227.329 + 106 7.324219e-03 9.765625e-05 32768 0 1234.036 + 107 7.421875e-03 9.765625e-05 32768 0 1249.199 + 108 7.519531e-03 9.765625e-05 32768 0 1258.591 + 109 7.617188e-03 9.765625e-05 32768 0 1253.583 + 110 7.714844e-03 9.765625e-05 32768 0 1250.293 + 111 7.812500e-03 9.765625e-05 32768 0 1265.752 + 112 7.910156e-03 9.765625e-05 32768 0 1247.696 + 113 8.007813e-03 9.765625e-05 32768 0 1608.997 + 114 8.105469e-03 9.765625e-05 32768 0 1505.135 + 115 8.203125e-03 9.765625e-05 32768 0 1508.400 + 116 8.300781e-03 9.765625e-05 32768 0 1524.896 + 117 8.398437e-03 9.765625e-05 32768 0 1481.166 + 118 8.496094e-03 9.765625e-05 32768 0 1520.336 + 119 8.593750e-03 9.765625e-05 32768 0 1494.457 + 120 8.691406e-03 9.765625e-05 32768 0 1517.346 + 121 8.789062e-03 9.765625e-05 32768 0 1514.266 + 122 8.886719e-03 9.765625e-05 32768 0 1466.870 + 123 8.984375e-03 9.765625e-05 32768 0 1535.258 + 124 9.082031e-03 9.765625e-05 32768 0 1512.334 + 125 9.179688e-03 9.765625e-05 32768 0 1514.996 + 126 9.277344e-03 9.765625e-05 32768 0 1499.448 + 127 9.375000e-03 9.765625e-05 32768 0 1530.981 + 128 9.472656e-03 9.765625e-05 32768 0 1551.520 + 129 9.570313e-03 9.765625e-05 32768 0 1515.622 + 130 9.667969e-03 9.765625e-05 32768 0 1507.469 + 131 9.765625e-03 9.765625e-05 32768 0 1493.654 + 132 9.863281e-03 9.765625e-05 32768 0 1510.786 + 133 9.960938e-03 9.765625e-05 32768 0 1489.244 + 134 1.005859e-02 9.765625e-05 32768 0 1530.756 + 135 1.015625e-02 9.765625e-05 32768 0 1684.936 + 136 1.025391e-02 9.765625e-05 32768 0 1445.483 + 137 1.035156e-02 9.765625e-05 32768 0 1457.336 + 138 1.044922e-02 9.765625e-05 32768 0 1485.821 + 139 1.054688e-02 9.765625e-05 32768 0 1470.178 + 140 1.064453e-02 9.765625e-05 32768 0 1475.452 + 141 1.074219e-02 9.765625e-05 32768 0 1480.344 + 142 1.083984e-02 9.765625e-05 32768 0 1494.713 + 143 1.093750e-02 9.765625e-05 32768 0 1491.871 + 144 1.103516e-02 9.765625e-05 32768 0 1490.843 + 145 1.113281e-02 9.765625e-05 32768 0 1566.565 + 146 1.123047e-02 9.765625e-05 32768 0 1482.148 + 147 1.132813e-02 9.765625e-05 32768 0 1510.317 + 148 1.142578e-02 9.765625e-05 32768 0 1522.768 + 149 1.152344e-02 9.765625e-05 32768 0 1483.167 + 150 1.162109e-02 9.765625e-05 32768 0 1496.622 + 151 1.171875e-02 9.765625e-05 32768 0 1439.106 + 152 1.181641e-02 9.765625e-05 32768 0 1512.526 + 153 1.191406e-02 9.765625e-05 32768 0 1502.762 + 154 1.201172e-02 9.765625e-05 32768 0 1533.897 + 155 1.210938e-02 9.765625e-05 32768 0 1544.583 + 156 1.220703e-02 9.765625e-05 32768 0 1603.579 + 157 1.230469e-02 9.765625e-05 32768 0 1482.876 + 158 1.240234e-02 9.765625e-05 32768 0 1464.410 + 159 1.250000e-02 9.765625e-05 32768 0 1453.033 + 160 1.259766e-02 9.765625e-05 32768 0 1439.090 + 161 1.269531e-02 9.765625e-05 32768 0 1441.893 + 162 1.279297e-02 9.765625e-05 32768 0 1438.807 + 163 1.289063e-02 9.765625e-05 32768 0 1474.443 + 164 1.298828e-02 9.765625e-05 32768 0 1488.145 + 165 1.308594e-02 9.765625e-05 32768 0 1484.710 + 166 1.318359e-02 9.765625e-05 32768 0 1482.595 + 167 1.328125e-02 9.765625e-05 32768 0 1524.337 + 168 1.337891e-02 9.765625e-05 32768 0 1430.587 + 169 1.347656e-02 9.765625e-05 32768 0 1429.264 + 170 1.357422e-02 9.765625e-05 32768 0 1466.449 + 171 1.367188e-02 9.765625e-05 32768 0 1502.286 + 172 1.376953e-02 9.765625e-05 32768 0 1455.993 + 173 1.386719e-02 9.765625e-05 32768 0 1435.804 + 174 1.396484e-02 9.765625e-05 32768 0 1444.008 + 175 1.406250e-02 9.765625e-05 32768 0 1428.647 + 176 1.416016e-02 9.765625e-05 32768 0 1436.682 + 177 1.425781e-02 9.765625e-05 32768 0 1433.351 + 178 1.435547e-02 9.765625e-05 32768 0 1438.288 + 179 1.445313e-02 9.765625e-05 32768 0 1428.455 + 180 1.455078e-02 9.765625e-05 32768 0 1508.945 + 181 1.464844e-02 9.765625e-05 32768 0 1422.767 + 182 1.474609e-02 9.765625e-05 32768 0 1418.819 + 183 1.484375e-02 9.765625e-05 32768 0 1433.828 + 184 1.494141e-02 9.765625e-05 32768 0 1423.189 + 185 1.503906e-02 9.765625e-05 32768 0 1426.707 + 186 1.513672e-02 9.765625e-05 32768 0 1426.201 + 187 1.523438e-02 9.765625e-05 32768 0 1445.366 + 188 1.533203e-02 9.765625e-05 32768 0 1444.843 + 189 1.542969e-02 9.765625e-05 32768 0 1438.595 + 190 1.552734e-02 9.765625e-05 32768 0 1435.407 + 191 1.562500e-02 9.765625e-05 32768 0 1434.259 + 192 1.572266e-02 9.765625e-05 32768 0 1440.294 + 193 1.582031e-02 9.765625e-05 32768 0 1429.044 + 194 1.591797e-02 9.765625e-05 32768 0 1446.566 + 195 1.601563e-02 9.765625e-05 32768 0 1132.670 + 196 1.611328e-02 9.765625e-05 32768 0 1006.622 + 197 1.621094e-02 9.765625e-05 32768 0 991.213 + 198 1.630859e-02 9.765625e-05 32768 0 1013.067 + 199 1.640625e-02 9.765625e-05 32768 0 1001.624 + 200 1.650391e-02 9.765625e-05 32768 0 1024.164 + 201 1.660156e-02 9.765625e-05 32768 0 1025.695 + 202 1.669922e-02 9.765625e-05 32768 0 1012.893 + 203 1.679687e-02 9.765625e-05 32768 0 999.274 + 204 1.689453e-02 9.765625e-05 32768 0 1010.005 + 205 1.699219e-02 9.765625e-05 32768 0 1011.801 + 206 1.708984e-02 9.765625e-05 32768 0 998.205 + 207 1.718750e-02 9.765625e-05 32768 0 1003.814 + 208 1.728516e-02 9.765625e-05 32768 0 998.934 + 209 1.738281e-02 9.765625e-05 32768 0 997.772 + 210 1.748047e-02 9.765625e-05 32768 0 1006.868 + 211 1.757812e-02 9.765625e-05 32768 0 1011.088 + 212 1.767578e-02 9.765625e-05 32768 0 1088.980 + 213 1.777344e-02 9.765625e-05 32768 0 988.543 + 214 1.787109e-02 9.765625e-05 32768 0 983.371 + 215 1.796875e-02 9.765625e-05 32768 0 984.862 + 216 1.806641e-02 9.765625e-05 32768 0 981.652 + 217 1.816406e-02 9.765625e-05 32768 0 987.319 + 218 1.826172e-02 9.765625e-05 32768 0 993.242 + 219 1.835938e-02 9.765625e-05 32768 0 983.838 + 220 1.845703e-02 9.765625e-05 32768 0 996.515 + 221 1.855469e-02 9.765625e-05 32768 0 947.203 + 222 1.865234e-02 9.765625e-05 32768 0 1029.318 + 223 1.875000e-02 9.765625e-05 32768 0 992.890 + 224 1.884766e-02 9.765625e-05 32768 0 998.988 + 225 1.894531e-02 9.765625e-05 32768 0 1027.067 + 226 1.904297e-02 9.765625e-05 32768 0 1001.001 + 227 1.914063e-02 9.765625e-05 32768 0 1028.753 + 228 1.923828e-02 9.765625e-05 32768 0 1136.533 + 229 1.933594e-02 9.765625e-05 32768 0 1022.059 + 230 1.943359e-02 9.765625e-05 32768 0 1015.718 + 231 1.953125e-02 9.765625e-05 32768 0 1013.532 + 232 1.962891e-02 9.765625e-05 32768 0 1022.330 + 233 1.972656e-02 9.765625e-05 32768 0 1019.488 + 234 1.982422e-02 9.765625e-05 32768 0 1026.366 + 235 1.992188e-02 9.765625e-05 32768 0 1019.859 + 236 2.001953e-02 9.765625e-05 32768 0 991.519 + 237 2.011719e-02 9.765625e-05 32768 0 1176.608 + 238 2.021484e-02 9.765625e-05 32768 0 991.395 + 239 2.031250e-02 9.765625e-05 32768 0 994.193 + 240 2.041016e-02 9.765625e-05 32768 0 996.700 + 241 2.050781e-02 9.765625e-05 32768 0 997.162 + 242 2.060547e-02 9.765625e-05 32768 0 1042.537 + 243 2.070313e-02 9.765625e-05 32768 0 1048.016 + 244 2.080078e-02 9.765625e-05 32768 0 1120.329 + 245 2.089844e-02 9.765625e-05 32768 0 1022.509 + 246 2.099609e-02 9.765625e-05 32768 0 1014.947 + 247 2.109375e-02 9.765625e-05 32768 0 1011.888 + 248 2.119141e-02 9.765625e-05 32768 0 951.152 + 249 2.128906e-02 9.765625e-05 32768 0 1017.106 + 250 2.138672e-02 9.765625e-05 32768 0 1013.423 + 251 2.148438e-02 9.765625e-05 32768 0 989.139 + 252 2.158203e-02 9.765625e-05 32768 0 1010.243 + 253 2.167969e-02 9.765625e-05 32768 0 993.511 + 254 2.177734e-02 9.765625e-05 32768 0 1011.160 + 255 2.187500e-02 9.765625e-05 32768 0 1065.467 + 256 2.197266e-02 9.765625e-05 32768 0 1022.221 + 257 2.207031e-02 9.765625e-05 32768 0 1003.495 + 258 2.216797e-02 9.765625e-05 32768 0 1010.803 + 259 2.226563e-02 9.765625e-05 32768 0 1039.478 + 260 2.236328e-02 9.765625e-05 32768 0 1101.938 + 261 2.246094e-02 9.765625e-05 32768 0 986.570 + 262 2.255859e-02 9.765625e-05 32768 0 993.218 + 263 2.265625e-02 9.765625e-05 32768 0 989.205 + 264 2.275391e-02 9.765625e-05 32768 0 983.149 + 265 2.285156e-02 9.765625e-05 32768 0 991.438 + 266 2.294922e-02 9.765625e-05 32768 0 1007.181 + 267 2.304688e-02 9.765625e-05 32768 0 1035.828 + 268 2.314453e-02 9.765625e-05 32768 0 1002.614 + 269 2.324219e-02 9.765625e-05 32768 0 1041.212 + 270 2.333984e-02 9.765625e-05 32768 0 993.185 + 271 2.343750e-02 9.765625e-05 32768 0 1025.189 + 272 2.353516e-02 9.765625e-05 32768 0 1013.712 + 273 2.363281e-02 9.765625e-05 32768 0 1003.209 + 274 2.373047e-02 9.765625e-05 32768 0 1010.710 + 275 2.382813e-02 9.765625e-05 32768 0 1025.100 + 276 2.392578e-02 9.765625e-05 32768 0 1057.269 + 277 2.402344e-02 9.765625e-05 32768 0 1124.723 + 278 2.412109e-02 9.765625e-05 32768 0 1025.496 + 279 2.421875e-02 9.765625e-05 32768 0 1025.231 + 280 2.431641e-02 9.765625e-05 32768 0 1024.641 + 281 2.441406e-02 9.765625e-05 32768 0 994.931 + 282 2.451172e-02 9.765625e-05 32768 0 998.965 + 283 2.460938e-02 9.765625e-05 32768 0 1040.216 + 284 2.470703e-02 9.765625e-05 32768 0 1022.758 + 285 2.480469e-02 9.765625e-05 32768 0 1032.167 + 286 2.490234e-02 9.765625e-05 32768 0 1028.403 + 287 2.500000e-02 9.765625e-05 32768 0 1038.900 + 288 2.509766e-02 9.765625e-05 32768 0 1001.048 + 289 2.519531e-02 9.765625e-05 32768 0 1029.682 + 290 2.529297e-02 9.765625e-05 32768 0 1047.342 + 291 2.539062e-02 9.765625e-05 32768 0 1036.171 + 292 2.548828e-02 9.765625e-05 32768 0 1030.387 + 293 2.558594e-02 9.765625e-05 32768 0 1052.380 + 294 2.568359e-02 9.765625e-05 32768 0 1096.088 + 295 2.578125e-02 9.765625e-05 32768 0 1015.682 + 296 2.587891e-02 9.765625e-05 32768 0 1003.744 + 297 2.597656e-02 9.765625e-05 32768 0 990.037 + 298 2.607422e-02 9.765625e-05 32768 0 1019.735 + 299 2.617188e-02 9.765625e-05 32768 0 1024.186 + 300 2.626953e-02 9.765625e-05 32768 0 994.061 + 301 2.636719e-02 9.765625e-05 32768 0 1020.561 + 302 2.646484e-02 9.765625e-05 32768 0 1005.136 + 303 2.656250e-02 9.765625e-05 32768 0 1001.353 + 304 2.666016e-02 9.765625e-05 32768 0 1000.047 + 305 2.675781e-02 9.765625e-05 32768 0 1007.191 + 306 2.685547e-02 9.765625e-05 32768 0 1009.866 + 307 2.695313e-02 9.765625e-05 32768 0 1008.549 + 308 2.705078e-02 9.765625e-05 32768 0 1046.824 + 309 2.714844e-02 9.765625e-05 32768 0 1010.715 + 310 2.724609e-02 9.765625e-05 32768 0 1042.475 + 311 2.734375e-02 9.765625e-05 32768 0 1137.615 + 312 2.744141e-02 9.765625e-05 32768 0 1031.113 + 313 2.753906e-02 9.765625e-05 32768 0 1028.340 + 314 2.763672e-02 9.765625e-05 32768 0 1031.166 + 315 2.773438e-02 9.765625e-05 32768 0 999.628 + 316 2.783203e-02 9.765625e-05 32768 0 1021.724 + 317 2.792969e-02 9.765625e-05 32768 0 1020.103 + 318 2.802734e-02 9.765625e-05 32768 0 1017.440 + 319 2.812500e-02 9.765625e-05 32768 0 1002.805 + 320 2.822266e-02 9.765625e-05 32768 0 1020.247 + 321 2.832031e-02 9.765625e-05 32768 0 1020.481 + 322 2.841797e-02 9.765625e-05 32768 0 1022.169 + 323 2.851563e-02 9.765625e-05 32768 0 1015.567 + 324 2.861328e-02 9.765625e-05 32768 0 1014.845 + 325 2.871094e-02 9.765625e-05 32768 0 1057.012 + 326 2.880859e-02 9.765625e-05 32768 0 1054.731 + 327 2.890625e-02 9.765625e-05 32768 0 1078.085 + 328 2.900391e-02 9.765625e-05 32768 0 1135.718 + 329 2.910156e-02 9.765625e-05 32768 0 1026.678 + 330 2.919922e-02 9.765625e-05 32768 0 1005.376 + 331 2.929688e-02 9.765625e-05 32768 0 1002.337 + 332 2.939453e-02 9.765625e-05 32768 0 996.464 + 333 2.949219e-02 9.765625e-05 32768 0 1010.166 + 334 2.958984e-02 9.765625e-05 32768 0 1013.649 + 335 2.968750e-02 9.765625e-05 32768 0 1015.403 + 336 2.978516e-02 9.765625e-05 32768 0 1008.705 + 337 2.988281e-02 9.765625e-05 32768 0 1026.790 + 338 2.998047e-02 9.765625e-05 32768 0 939.948 + 339 3.007813e-02 9.765625e-05 32768 0 897.359 + 340 3.017578e-02 9.765625e-05 32768 0 1195.256 + 341 3.027344e-02 9.765625e-05 32768 0 1090.396 + 342 3.037109e-02 9.765625e-05 32768 0 1102.880 + 343 3.046875e-02 9.765625e-05 32768 0 1099.729 + 344 3.056641e-02 9.765625e-05 32768 0 1082.379 + 345 3.066406e-02 9.765625e-05 32768 0 1163.097 + 346 3.076172e-02 9.765625e-05 32768 0 1087.996 + 347 3.085938e-02 9.765625e-05 32768 0 1079.661 + 348 3.095703e-02 9.765625e-05 32768 0 1104.227 + 349 3.105469e-02 9.765625e-05 32768 0 1101.164 + 350 3.115234e-02 9.765625e-05 32768 0 1120.536 + 351 3.125000e-02 9.765625e-05 32768 0 1088.789 + 352 3.134766e-02 9.765625e-05 32768 0 1068.192 + 353 3.144531e-02 9.765625e-05 32768 0 1071.413 + 354 3.154297e-02 9.765625e-05 32768 0 1089.774 + 355 3.164062e-02 9.765625e-05 32768 0 1078.393 + 356 3.173828e-02 9.765625e-05 32768 0 1109.767 + 357 3.183594e-02 9.765625e-05 32768 0 1106.486 + 358 3.193359e-02 9.765625e-05 32768 0 1109.993 + 359 3.203125e-02 9.765625e-05 32768 0 1108.227 + 360 3.212891e-02 9.765625e-05 32768 0 1109.371 + 361 3.222656e-02 9.765625e-05 32768 0 1115.127 + 362 3.232422e-02 9.765625e-05 32768 0 1123.371 + 363 3.242188e-02 9.765625e-05 32768 0 992.571 + 364 3.251953e-02 9.765625e-05 32768 0 1040.697 + 365 3.261719e-02 9.765625e-05 32768 0 1025.617 + 366 3.271484e-02 9.765625e-05 32768 0 1006.040 + 367 3.281250e-02 9.765625e-05 32768 0 1008.713 + 368 3.291016e-02 9.765625e-05 32768 0 1020.951 + 369 3.300781e-02 9.765625e-05 32768 0 1015.800 + 370 3.310547e-02 9.765625e-05 32768 0 1037.745 + 371 3.320312e-02 9.765625e-05 32768 0 1040.416 + 372 3.330078e-02 9.765625e-05 32768 0 1030.467 + 373 3.339844e-02 9.765625e-05 32768 0 1049.299 + 374 3.349609e-02 9.765625e-05 32768 0 1059.101 + 375 3.359375e-02 9.765625e-05 32768 0 1052.807 + 376 3.369141e-02 9.765625e-05 32768 0 1027.797 + 377 3.378906e-02 9.765625e-05 32768 0 1037.017 + 378 3.388672e-02 9.765625e-05 32768 0 1043.442 + 379 3.398438e-02 9.765625e-05 32768 0 1052.545 + 380 3.408203e-02 9.765625e-05 32768 0 1125.092 + 381 3.417969e-02 9.765625e-05 32768 0 1012.641 + 382 3.427734e-02 9.765625e-05 32768 0 1005.484 + 383 3.437500e-02 9.765625e-05 32768 0 1009.903 + 384 3.447266e-02 9.765625e-05 32768 0 1018.198 + 385 3.457031e-02 9.765625e-05 32768 0 1003.874 + 386 3.466797e-02 9.765625e-05 32768 0 1028.012 + 387 3.476563e-02 9.765625e-05 32768 0 1021.830 + 388 3.486328e-02 9.765625e-05 32768 0 1022.930 + 389 3.496094e-02 9.765625e-05 32768 0 1010.233 + 390 3.505859e-02 9.765625e-05 32768 0 1039.573 + 391 3.515625e-02 9.765625e-05 32768 0 1020.382 + 392 3.525391e-02 9.765625e-05 32768 0 1044.217 + 393 3.535156e-02 9.765625e-05 32768 0 1043.863 + 394 3.544922e-02 9.765625e-05 32768 0 1030.654 + 395 3.554687e-02 9.765625e-05 32768 0 1024.336 + 396 3.564453e-02 9.765625e-05 32768 0 1048.685 + 397 3.574219e-02 9.765625e-05 32768 0 1070.174 + 398 3.583984e-02 9.765625e-05 32768 0 1103.300 + 399 3.593750e-02 9.765625e-05 32768 0 1116.885 + 400 3.603516e-02 9.765625e-05 32768 0 1025.508 + 401 3.613281e-02 9.765625e-05 32768 0 1024.324 + 402 3.623047e-02 9.765625e-05 32768 0 1045.641 + 403 3.632813e-02 9.765625e-05 32768 0 1052.298 + 404 3.642578e-02 9.765625e-05 32768 0 1049.651 + 405 3.652344e-02 9.765625e-05 32768 0 1039.583 + 406 3.662109e-02 9.765625e-05 32768 0 927.192 + 407 3.671875e-02 9.765625e-05 32768 0 952.513 + 408 3.681641e-02 9.765625e-05 32768 0 1048.747 + 409 3.691406e-02 9.765625e-05 32768 0 1053.193 + 410 3.701172e-02 9.765625e-05 32768 0 1032.866 + 411 3.710938e-02 9.765625e-05 32768 0 1057.847 + 412 3.720703e-02 9.765625e-05 32768 0 1083.408 + 413 3.730469e-02 9.765625e-05 32768 0 1113.171 + 414 3.740234e-02 9.765625e-05 32768 0 1116.059 + 415 3.750000e-02 9.765625e-05 32768 0 1097.219 + 416 3.759766e-02 9.765625e-05 32768 0 1097.206 + 417 3.769531e-02 9.765625e-05 32768 0 1091.543 diff --git a/examples/Feedback/makeIC.py b/examples/Feedback/makeIC.py new file mode 100644 index 0000000000000000000000000000000000000000..d494a99e456d11019aae1450233125ffaf07f7ad --- /dev/null +++ b/examples/Feedback/makeIC.py @@ -0,0 +1,108 @@ +############################################################################### + # This file is part of SWIFT. + # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk), + # Matthieu Schaller (matthieu.schaller@durham.ac.uk) + # + # This program is free software: you can redistribute it and/or modify + # it under the terms of the GNU Lesser General Public License as published + # by the Free Software Foundation, either version 3 of the License, or + # (at your option) any later version. + # + # This program is distributed in the hope that it will be useful, + # but WITHOUT ANY WARRANTY; without even the implied warranty of + # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + # GNU General Public License for more details. + # + # You should have received a copy of the GNU Lesser General Public License + # along with this program. If not, see <http://www.gnu.org/licenses/>. + # + ############################################################################## + +import h5py +import sys +from numpy import * + +# Generates a swift IC file containing a cartesian distribution of particles +# at a constant density and pressure in a cubic box + +# Parameters +periodic= 1 # 1 For periodic box +boxSize = 1. +L = int(sys.argv[1]) # Number of particles along one axis +rho = 1. # Density +P = 1.e-6 # Pressure +gamma = 5./3. # Gas adiabatic index +eta = 1.2349 # 48 ngbs with cubic spline kernel +fileName = "Feedback.hdf5" + +#--------------------------------------------------- +numPart = L**3 +mass = boxSize**3 * rho / numPart +internalEnergy = P / ((gamma - 1.)*rho) + +#-------------------------------------------------- + +#File +file = h5py.File(fileName, 'w') + +# Header +grp = file.create_group("/Header") +grp.attrs["BoxSize"] = boxSize +grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0] +grp.attrs["Time"] = 0.0 +grp.attrs["NumFilesPerSnapshot"] = 1 +grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +grp.attrs["Flag_Entropy_ICs"] = 0 + +#Runtime parameters +grp = file.create_group("/RuntimePars") +grp.attrs["PeriodicBoundariesOn"] = periodic + +#Units +grp = file.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = 1. +grp.attrs["Unit mass in cgs (U_M)"] = 1. +grp.attrs["Unit time in cgs (U_t)"] = 1. +grp.attrs["Unit current in cgs (U_I)"] = 1. +grp.attrs["Unit temperature in cgs (U_T)"] = 1. + +#Particle group +grp = file.create_group("/PartType0") + +v = zeros((numPart, 3)) +ds = grp.create_dataset('Velocities', (numPart, 3), 'f') +ds[()] = v +v = zeros(1) + +m = full((numPart, 1), mass) +ds = grp.create_dataset('Masses', (numPart,1), 'f') +ds[()] = m +m = zeros(1) + +h = full((numPart, 1), eta * boxSize / L) +ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f') +ds[()] = h +h = zeros(1) + +u = full((numPart, 1), internalEnergy) +ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f') +ds[()] = u +u = zeros(1) + + +ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1)) +ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L') +ds[()] = ids + 1 +x = ids % L; +y = ((ids - x) / L) % L; +z = (ids - x - L * y) / L**2; +coords = zeros((numPart, 3)) +coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L) +coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L) +coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L) +ds = grp.create_dataset('Coordinates', (numPart, 3), 'd') +ds[()] = coords + +file.close() diff --git a/examples/main.c b/examples/main.c index efd40c3da9802ebfab72a51eb841c035c8490eba..d4a920db0a56066c7f89e172d15094dbfbfd4d28 100644 --- a/examples/main.c +++ b/examples/main.c @@ -73,6 +73,8 @@ void print_help_message() { "Overwrite the CPU frequency (Hz) to be used for time measurements"); printf(" %2s %8s %s\n", "-g", "", "Run with an external gravitational potential"); + printf(" %2s %8s %s\n", "-F", "", + "Run with feedback "); printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity"); printf(" %2s %8s %s\n", "-n", "{int}", "Execute a fixed number of time steps"); @@ -143,6 +145,7 @@ int main(int argc, char *argv[]) { int nsteps = -2; int with_cosmology = 0; int with_external_gravity = 0; + int with_sourceterms = 0; int with_self_gravity = 0; int with_hydro = 0; int with_fp_exceptions = 0; @@ -153,7 +156,7 @@ int main(int argc, char *argv[]) { /* Parse the parameters */ int c; - while ((c = getopt(argc, argv, "acdef:gGhn:st:v:y:")) != -1) switch (c) { + while ((c = getopt(argc, argv, "acdef:gFGhn:st:v:y:")) != -1) switch (c) { case 'a': with_aff = 1; break; @@ -178,6 +181,9 @@ int main(int argc, char *argv[]) { break; case 'G': with_self_gravity = 1; + break; + case 'F': + with_sourceterms = 1; break; case 'h': if (myrank == 0) print_help_message(); @@ -335,6 +341,11 @@ int main(int argc, char *argv[]) { if (with_external_gravity) potential_init(params, &us, &potential); if (with_external_gravity && myrank == 0) potential_print(&potential); + /* Initialise the feedback properties */ + struct sourceterms sourceterms; + if (with_sourceterms) source_terms_init(params, &us, &sourceterms); + if (with_sourceterms && myrank == 0) source_terms_print(&sourceterms); + /* Read particles and space information from (GADGET) ICs */ char ICfileName[200] = ""; parser_get_param_string(params, "InitialConditions:file_name", ICfileName); @@ -440,6 +451,7 @@ int main(int argc, char *argv[]) { if (with_hydro) engine_policies |= engine_policy_hydro; if (with_self_gravity) engine_policies |= engine_policy_self_gravity; if (with_external_gravity) engine_policies |= engine_policy_external_gravity; + if (with_sourceterms) engine_policies |= engine_policy_sourceterms; if (with_cosmology) engine_policies |= engine_policy_cosmology; /* Initialize the engine with the space and policies. */ @@ -447,7 +459,7 @@ int main(int argc, char *argv[]) { struct engine e; engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff, engine_policies, talking, &us, &prog_const, &hydro_properties, - &potential); + &potential, &sourceterms); if (myrank == 0) { clocks_gettime(&toc); message("engine_init took %.3f %s.", clocks_diff(&tic, &toc), diff --git a/src/Makefile.am b/src/Makefile.am index 4d6a67c7569464486a017d8306c3e54730ebb3b2..e1f0235cab7bb6111ab42bf69292fabc88f63fd8 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -42,7 +42,8 @@ endif include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \ common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \ - physical_constants.h physical_constants_cgs.h potentials.h version.h hydro_properties.h + physical_constants.h physical_constants_cgs.h potentials.h version.h hydro_properties.h \ + sourceterms.h # Common source files AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ @@ -50,7 +51,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ units.c common_io.c single_io.c multipole.c version.c map.c \ kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \ physical_constants.c potentials.c hydro_properties.c \ - runner_doiact_fft.c + runner_doiact_fft.c sourceterms.c # Include files for distribution, not installation. nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ @@ -60,6 +61,7 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h gravity.h gravity_io.h \ gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \ + sourceterms.h \ hydro.h hydro_io.h \ hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \ hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \ diff --git a/src/cell.h b/src/cell.h index 150718eeb4bd3857e37d517718fe53661033a330..a37f5d9e1bfa7108d9908686a966925f751dba97 100644 --- a/src/cell.h +++ b/src/cell.h @@ -136,6 +136,9 @@ struct cell { /* Task for external gravity */ struct task *grav_external; + /* Task for source terms */ + struct task *sourceterms; + /* Number of tasks that are associated with this cell. */ int nr_tasks; diff --git a/src/const.h b/src/const.h index 6710fe9e52dfdbede468282d84d4441b923a59fe..7f882f828b5df491f6d2f88467b2acc383d8684d 100644 --- a/src/const.h +++ b/src/const.h @@ -73,6 +73,9 @@ #define EXTERNAL_POTENTIAL_POINTMASS //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL +/* Source terms */ +#define SN_FEEDBACK + /* Are we debugging ? */ //#define SWIFT_DEBUG_CHECKS diff --git a/src/engine.c b/src/engine.c index a980bcbc9c2f42127fb43b30c7d5573a8e9957ef..91a82bab529d250811bf270c5c3d6f587d108e4d 100644 --- a/src/engine.c +++ b/src/engine.c @@ -68,7 +68,7 @@ #include "units.h" #include "version.h" -const char *engine_policy_names[13] = {"none", +const char *engine_policy_names[14] = {"none", "rand", "steal", "keep", @@ -80,7 +80,8 @@ const char *engine_policy_names[13] = {"none", "hydro", "self_gravity", "external_gravity", - "cosmology_integration"}; + "cosmology_integration", + "sourceterms"}; /** The rank of the engine as a global variable (for messages). */ int engine_rank; @@ -163,6 +164,7 @@ void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c, if (is_with_external_gravity) c->grav_external = scheduler_addtask( s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0); + } } @@ -191,6 +193,7 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c, struct scheduler *s = &e->sched; const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt; + const int is_with_sourceterms = (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms; /* Is this the super-cell? */ if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) { @@ -225,6 +228,11 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c, /* Generate the ghost task. */ c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0, c, NULL, 0); + + /* add source terms */ + if (is_with_sourceterms) + c->sourceterms = scheduler_addtask(s, task_type_sourceterms, task_subtype_none, 0, 0, + c, NULL, 0); } } @@ -1610,6 +1618,10 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) { else if (t->type == task_type_grav_external) { scheduler_addunlock(sched, t->ci->init, t); scheduler_addunlock(sched, t, t->ci->kick); + } + /* source terms depend on kick (should eventually depend on cooling) */ + else if (t->type == task_type_sourceterms) { + scheduler_addunlock(sched, t->ci->kick, t); } } } @@ -2602,6 +2614,11 @@ void engine_step(struct engine *e) { mask |= 1 << task_type_grav_external; } + /* Add the tasks corresponding to sourceterms to the masks */ + if (e->policy & engine_policy_sourceterms) { + mask |= 1 << task_type_sourceterms; + } + /* Add MPI tasks if need be */ if (e->policy & engine_policy_mpi) { @@ -2927,7 +2944,8 @@ void engine_init(struct engine *e, struct space *s, const struct UnitSystem *internal_units, const struct phys_const *physical_constants, const struct hydro_props *hydro, - const struct external_potential *potential) { + const struct external_potential *potential, + const struct sourceterms *sourceterms) { /* Clean-up everything */ bzero(e, sizeof(struct engine)); @@ -2976,6 +2994,7 @@ void engine_init(struct engine *e, struct space *s, e->physical_constants = physical_constants; e->hydro_properties = hydro; e->external_potential = potential; + e->sourceterms = sourceterms; e->parameter_file = params; engine_rank = nodeID; diff --git a/src/engine.h b/src/engine.h index d708198c32b67c5118bbd7f4676f1ea0fe821c7d..164b0f1202f8760be92993fc45c280b8e6006d4b 100644 --- a/src/engine.h +++ b/src/engine.h @@ -41,6 +41,7 @@ #include "parser.h" #include "partition.h" #include "potentials.h" +#include "sourceterms.h" #include "runner.h" #include "scheduler.h" #include "space.h" @@ -61,7 +62,8 @@ enum engine_policy { engine_policy_hydro = (1 << 8), engine_policy_self_gravity = (1 << 9), engine_policy_external_gravity = (1 << 10), - engine_policy_cosmology = (1 << 11) + engine_policy_cosmology = (1 << 11), + engine_policy_sourceterms = (1 << 12) }; extern const char *engine_policy_names[]; @@ -204,6 +206,9 @@ struct engine { /* Properties of external gravitational potential */ const struct external_potential *external_potential; + /* Properties of source terms */ + const struct sourceterms *sourceterms; + /* The (parsed) parameter file */ const struct swift_params *parameter_file; }; @@ -218,7 +223,8 @@ void engine_init(struct engine *e, struct space *s, const struct UnitSystem *internal_units, const struct phys_const *physical_constants, const struct hydro_props *hydro, - const struct external_potential *potential); + const struct external_potential *potential, + const struct sourceterms *sourceterms); void engine_launch(struct engine *e, int nr_runners, unsigned int mask, unsigned int submask); void engine_prepare(struct engine *e); diff --git a/src/equation_of_state.h b/src/equation_of_state.h index b4a36e8a3ef0bcc281d1f939f89fde08ecf00be9..74cae84b31434a89cd08dae26df69bf85bd43f64 100644 --- a/src/equation_of_state.h +++ b/src/equation_of_state.h @@ -40,7 +40,9 @@ /* ------------------------------------------------------------------------- */ #if defined(EOS_IDEAL_GAS) - +#if defined(EOS_ISOTHERMAL_GAS) +#error can't have both! +#endif /** * @brief Returns the internal energy given density and entropy * diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index d2d4450fa12374a8a8dec624c5e54ba3d47b99aa..41acfc22e242c6fad5cba8cdc5d2aae9611be6f3 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -16,7 +16,8 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ - +#ifndef SWIFT_GADGET_HYDRO_H +#define SWIFT_GADGET_HYDRO_H #include "adiabatic_index.h" #include "dimension.h" #include "equation_of_state.h" @@ -361,3 +362,4 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( /* We read u in the entropy field. We now get S from u */ p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy); } +#endif /* SWIFT_GADGET_HYDRO_H */ diff --git a/src/runner.c b/src/runner.c index 6e99d8cece2cf086b1209ac2f8314301b9779e20..f8e03fba5c5ba674ca9c911c5fa4dde512dc00dc 100644 --- a/src/runner.c +++ b/src/runner.c @@ -47,6 +47,7 @@ #include "engine.h" #include "error.h" #include "gravity.h" +#include "sourceterms.h" #include "hydro.h" #include "hydro_properties.h" #include "kick.h" @@ -91,6 +92,70 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, #include "runner_doiact_fft.h" #include "runner_doiact_grav.h" +/** + * @brief Calculate gravity acceleration from external potential + * + * @param r runner task + * @param c cell + * @param timer 1 if the time is to be recorded. + */ +void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) { + + struct part *restrict parts = c->parts; + const int count = c->count; + const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]}; + const double cell_width[3] = {c->width[0], c->width[1], c->width[2]}; + const int ti_current = r->e->ti_current; + const struct sourceterms *sourceterms = r->e->sourceterms; + const double location[3] = {sourceterms->supernova.x, sourceterms->supernova.y, sourceterms->supernova.z}; + const int dimen=3; + const double timeBase = r->e->timeBase; + + TIMER_TIC; + + /* Recurse? */ + if (c->split) { + for (int k = 0; k < 8; k++) + if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0); + return; + } + + /* does cell contain explosion? */ + if(count > 0) { + const int incell = is_in_cell(cell_min, cell_width, location, dimen); + if (incell == 1) + { + + /* inject SN energy into particle with highest id, if it is active */ + int imax = 0; + struct part *restrict sn_part; + for (int i = 0; i < count; i++) { + + /* Get a direct pointer on the part. */ + struct part *restrict p = &parts[i]; + if(p->id > imax) + { + imax = p->id; + sn_part = p; + } + } + /* Is this part within the time step? */ + if (sn_part->ti_begin == ti_current) { + + /* does this time step sraddle the feedback injecton time? */ + const float t_begin = sn_part->ti_begin * timeBase; + const float t_end = sn_part->ti_end * timeBase; + if(t_begin <= sourceterms->supernova.time && t_end > sourceterms->supernova.time) + { + do_supernova_feedback(sourceterms, sn_part); + } + } + } + } + + if (timer) TIMER_TOC(timer_dosource); +} + /** * @brief Calculate gravity acceleration from external potential * @@ -1158,6 +1223,9 @@ void *runner_main(void *data) { case task_type_grav_external: runner_do_grav_external(r, t->ci, 1); break; + case task_type_sourceterms: + runner_do_sourceterms(r, t->ci, 1); + break; case task_type_part_sort: space_do_parts_sort(); break; diff --git a/src/sourceterms.c b/src/sourceterms.c new file mode 100644 index 0000000000000000000000000000000000000000..a3d32d829a02910ce9d98c6dd14ed905b21e2626 --- /dev/null +++ b/src/sourceterms.c @@ -0,0 +1,76 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk) + * Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +/* Config parameters. */ +#include "../config.h" + +/* Some standard headers. */ +#include <math.h> + +/* Local includes. */ +#include "const.h" +#include "error.h" +#include "parser.h" +#include "part.h" +#include "physical_constants.h" +#include "equation_of_state.h" +#include "units.h" + +/* This object's header. */ +#include "sourceterms.h" + +/** + * @brief Initialises the source terms + * + * @param parameter_file The parsed parameter file + * @param us The current internal system of units + * @param sourceterms: the structure that has all the source term properties + */ +void source_terms_init(const struct swift_params* parameter_file, + struct UnitSystem* us, + struct sourceterms* source) { + +#ifdef SN_FEEDBACK + source->supernova.time = parser_get_param_double(parameter_file, "SN:time"); + source->supernova.energy = parser_get_param_double(parameter_file, "SN:energy"); + source->supernova.x = parser_get_param_double(parameter_file, "SN:x"); + source->supernova.y = parser_get_param_double(parameter_file, "SN:y"); + source->supernova.z = parser_get_param_double(parameter_file, "SN:z"); +#endif /* SN_FEEDBCK */ +}; + +/** + * @brief Prints the properties of the external potential to stdout. + * + * @param potential The external potential properties. + */ +void source_terms_print(const struct sourceterms* source) { + +#ifdef SN_FEEDBACK + message(" Single SNe of energy= %e will explode at time= %e at location (%e,%e,%e)", + source->supernova.energy, source->supernova.time, source->supernova.x, source->supernova.y, source->supernova.z); +#endif /* SN_FEEDBACK */ + +}; + +__attribute__((always_inline)) INLINE float calculate_entropy(const float entropy_old, const float density, const float u_old, const float u_new) +{ + return entropy_old + hydro_gamma_minus_one * (u_new-u_old) * pow_minus_gamma_minus_one(density); +}; diff --git a/src/sourceterms.h b/src/sourceterms.h new file mode 100644 index 0000000000000000000000000000000000000000..506aadde82e66d9bc67c5bfe5bccb020d1dc7358 --- /dev/null +++ b/src/sourceterms.h @@ -0,0 +1,39 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_SOURCETERMS_H +#define SWIFT_SOURCETERMS_H + +#include "./const.h" +/* So far only one model here */ +struct sourceterms +{ +#ifdef SN_FEEDBACK +#include "sourceterms/sn_feedback/sn_feedback_struct.h" +#endif +}; + +void source_terms_init(const struct swift_params* parameter_file, + struct UnitSystem* us, + struct sourceterms* source); +void source_terms_print(const struct sourceterms* source); +float calculate_entropy(const float entropy_old, const float density, const float u_old, const float u_new); +#ifdef SN_FEEDBACK +#include "sourceterms/sn_feedback/sn_feedback.h" +#endif +#endif /* SWIFT_SOURCETERMS_H */ diff --git a/src/sourceterms/Default/sourceterms.h b/src/sourceterms/Default/sourceterms.h new file mode 100644 index 0000000000000000000000000000000000000000..a31a675fbe36629fe79849691f8d0377626cb5b5 --- /dev/null +++ b/src/sourceterms/Default/sourceterms.h @@ -0,0 +1,51 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * 2016 Tom Theuns (tom.theuns@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +#include <float.h> +#include "feedback.h" + +/** + * @brief Computes the feedback time-step of a given particle due to + * a source term + * + * This function only branches towards the potential chosen by the user. + * + * @param feedback The properties of the source terms. + * @param phys_const The physical constants in internal units. + * @param g Pointer to the particle data. + */ +__attribute__((always_inline)) INLINE static float +sourceterms_compute_timestep(const struct sourceterms *feedback, + const struct phys_const* const phys_const, + const struct part* const p){ + float dt = FLT_MAX; +#ifdef SN_FEEDBACK + dt = + fmin(dt, sn_feedback_timestep(feedback, phys_const, p)); +#endif +} +__attribute__((always_inline)) INLINE static void feedback(const struct sourceterms *feedback, + const struct phys_const* const phys_const, + struct part* p) +#ifdef SN_FEEDBACK + sn_feedback(feedback, phys_const, p); +#endif +} + diff --git a/src/sourceterms/Default/sourceterms_iact.h b/src/sourceterms/Default/sourceterms_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/src/sourceterms/sn_feedback/sn_feedback.c b/src/sourceterms/sn_feedback/sn_feedback.c new file mode 100644 index 0000000000000000000000000000000000000000..ebc9f69c55eb39bf300cef7b10ae157d97da97cb --- /dev/null +++ b/src/sourceterms/sn_feedback/sn_feedback.c @@ -0,0 +1,64 @@ +#error: this file is no longer in use! +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk) + * Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * Richard Bower (r.g.bower@durham.ac.uk) + * Stefan Arridge (stefan.arridge@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +/* Config parameters. */ +#include "../config.h" +#include "../sourceterms.h" + +void supernova_feedback_init(const struct swift_params* parameter_file, + struct UnitSystem* us, + struct supernova* sn){ + sn.time = parser_get_param_double(parameter_file, "SN:time"); + sn.energy = parser_get_param_double(parameter_file, "SN:energy"); + sn.x = parser_get_param_double(parameter_file, "SN:x"); + sn.y = parser_get_param_double(parameter_file, "SN:y"); + sn.z = parser_get_param_double(parameter_file, "SN:z"); +} + +void supernova_feedback_print(const struct supernova* sn){ + message(" Single SNe of energy= %e will explode at time= %e at location (%e,%e,%e)", + sn.energy, sn.time, sn.x, sn.y, sn.z); +}; + +__attribute__((always_inline)) INLINE static void do_supernova_feedback(const struct sourceterms* sourceterms, struct part* p) +{ +}; + +__attribute__((always_inline)) INLINE void update_entropy(const sourceterms *sourceterms, struct part* p) +{ + + /*updates the entropy of a particle due to feedback */ + float u_old; + float u_new; + float new_entropy; + float old_entropy = p->entropy; + float rho = p->rho; + + // u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1); + const float u_old = hydro_get_internal_energy(p,0); // dt = 0 because using current entropy + const float u_new = u_old + sourceterms->supernova.energy; + const float new_entropy = u_new*pow_minus_gamma_minus_one(-p>rho) * hydro_gamma_minus_one; + p->entropy = new_entropy; +} + + diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h new file mode 100644 index 0000000000000000000000000000000000000000..16a60c7b6f7bc10af00755281d660df0233e3696 --- /dev/null +++ b/src/sourceterms/sn_feedback/sn_feedback.h @@ -0,0 +1,59 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk) + * Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * Richard Bower (r.g.bower@durham.ac.uk) + * Stefan Arridge (stefan.arridge@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ +#ifndef SWIFT_SN_FEEDBACK_H +#define SWIFT_SN_FEEDBACK_H +#include <float.h> +#include "equation_of_state.h" +#include "hydro.h" +// void do_supernova_feedback(const struct source_terms* source_terms, const struct phys_const* const constants, const struct part* p); + + +/* determine whether location is in cell */ +__attribute__((always_inline)) INLINE static int is_in_cell(const double cell_min[], const double cell_width[], const double location[], const int dimen) +{ + for(int i=0; i<dimen; i++){ + if (cell_min[i] > location[i]) + return 0; + if ( (cell_min[i] + cell_width[i]) <= location[i]) + return 0; + } + return 1; +}; + +__attribute__((always_inline)) INLINE static void do_supernova_feedback(const struct sourceterms* sourceterms, struct part* p) +{ + const float density = p->rho; + const float u_old = hydro_get_internal_energy(p,0); + const float u_new = u_old + sourceterms->supernova.energy / p->mass; + const float entropy_old = gas_entropy_from_internal_energy(density, u_old); + const float entropy_new = calculate_entropy(entropy_old, density, u_old, u_new); + hydro_set_internal_energy(p, u_new); + hydro_set_entropy(p,entropy_new); + message(" u_old= %e u_new= %e S_old= %e S_new= %e rho= %e", u_old, u_new, entropy_old, entropy_new, density); + message(" p->s %e", p->entropy); +}; +#endif /* SWIFT_SN_FEEDBACK_H */ + + + + + diff --git a/src/sourceterms/sn_feedback/sn_feedback_struct.h b/src/sourceterms/sn_feedback/sn_feedback_struct.h new file mode 100644 index 0000000000000000000000000000000000000000..fe662e49c77f6870e20a9dc8b34815b4f2ac1f8d --- /dev/null +++ b/src/sourceterms/sn_feedback/sn_feedback_struct.h @@ -0,0 +1,7 @@ +/* sn feedback */ +struct { + double time; + double energy; + double x, y, z; + } supernova; + diff --git a/src/swift.h b/src/swift.h index 7e3116c1de8bc8e6cc2f89d0d5cbe9771ffbf33a..6d8d57968f1fd8283044f2dbb753d2943f0d43eb 100644 --- a/src/swift.h +++ b/src/swift.h @@ -43,6 +43,7 @@ #include "partition.h" #include "physical_constants.h" #include "potentials.h" +#include "sourceterms.h" #include "queue.h" #include "runner.h" #include "scheduler.h" diff --git a/src/task.c b/src/task.c index 13dd47e6cbf68de4ea6cd8ba6b898ee41a06618d..54b9ec5fae4e5d6040eb0f65194e1b4280a0f1a1 100644 --- a/src/task.c +++ b/src/task.c @@ -111,6 +111,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on( case task_type_sort: case task_type_ghost: + case task_type_sourceterms: return task_action_part; break; diff --git a/src/task.h b/src/task.h index a7cbf28c3d1c7bde45102e5ce85e18c5f736e343..a5cf89c9da62a20e3f0e356e81beee26d756255b 100644 --- a/src/task.h +++ b/src/task.h @@ -51,6 +51,7 @@ enum task_types { task_type_grav_mm, task_type_grav_up, task_type_grav_external, + task_type_sourceterms, task_type_part_sort, task_type_gpart_sort, task_type_split_cell, diff --git a/src/timers.h b/src/timers.h index aa8455397daf1e88b709a8332e3ae63694991e94..81a6df3786a0b2e2f423e869fe5b7ae2efd6c009 100644 --- a/src/timers.h +++ b/src/timers.h @@ -43,6 +43,7 @@ enum { timer_dopair_grav_pm, timer_dopair_grav_pp, timer_dograv_external, + timer_dosource, timer_dosub_self_density, timer_dosub_self_force, timer_dosub_self_grav,