From 0d354af02f058c52fd473666f4dd2901c75cff51 Mon Sep 17 00:00:00 2001 From: Sam Skillman Date: Tue, 21 Jul 2009 15:34:32 -0400 Subject: [PATCH] Adding necessary shock/cosmic ray code. --HG-- branch : week-of-code --- input/cosmic_ray.dat | 1024 +++++++++++++ src/enzo/CosmicRayData.h | 38 + src/enzo/Grid_FindShocks.C | 1854 +++++++++++++++++++++++ src/enzo/Grid_IdentifyCRSpeciesFields.C | 50 + src/enzo/Grid_ShocksHandler.C | 55 + src/enzo/InitializeCosmicRayData.C | 93 ++ 6 files changed, 3114 insertions(+) create mode 100644 input/cosmic_ray.dat create mode 100644 src/enzo/CosmicRayData.h create mode 100644 src/enzo/Grid_FindShocks.C create mode 100644 src/enzo/Grid_IdentifyCRSpeciesFields.C create mode 100644 src/enzo/Grid_ShocksHandler.C create mode 100644 src/enzo/InitializeCosmicRayData.C diff --git a/input/cosmic_ray.dat b/input/cosmic_ray.dat new file mode 100644 index 000000000..f8d678ed5 --- /dev/null +++ b/input/cosmic_ray.dat @@ -0,0 +1,1024 @@ +0.000000e+00 1.000000e+00 0.000000e+00 +0.000000e+00 1.013583e+00 5.350931e-05 +0.000000e+00 1.027351e+00 1.084816e-04 +0.000000e+00 1.041305e+00 1.649577e-04 +0.000000e+00 1.055450e+00 2.229785e-04 +0.000000e+00 1.069786e+00 2.825863e-04 +0.000000e+00 1.084317e+00 3.438243e-04 +0.000000e+00 1.099045e+00 4.067371e-04 +0.000000e+00 1.113974e+00 4.713705e-04 +0.000000e+00 1.129105e+00 5.377718e-04 +0.000000e+00 1.144442e+00 6.059892e-04 +0.000000e+00 1.159987e+00 6.760726e-04 +0.000000e+00 1.175743e+00 7.480729e-04 +0.000000e+00 1.191714e+00 8.220426e-04 +0.000000e+00 1.207901e+00 8.980352e-04 +0.000000e+00 1.224308e+00 9.761057e-04 +0.000000e+00 1.240938e+00 1.056312e-03 +0.000000e+00 1.257794e+00 1.138712e-03 +0.000000e+00 1.274878e+00 1.223366e-03 +0.000000e+00 1.292195e+00 1.310335e-03 +0.000000e+00 1.309747e+00 1.399682e-03 +0.000000e+00 1.327538e+00 1.491474e-03 +0.000000e+00 1.345570e+00 1.585776e-03 +0.000000e+00 1.363847e+00 1.682657e-03 +0.000000e+00 1.382372e+00 1.782189e-03 +0.000000e+00 1.401149e+00 1.884441e-03 +0.000000e+00 1.420181e+00 1.989492e-03 +0.000000e+00 1.439472e+00 2.097415e-03 +0.000000e+00 1.459024e+00 2.208289e-03 +0.000000e+00 1.478842e+00 2.322197e-03 +0.000000e+00 1.498930e+00 2.439220e-03 +0.000000e+00 1.519290e+00 2.559444e-03 +0.000000e+00 1.539927e+00 2.682955e-03 +0.000000e+00 1.560844e+00 2.809846e-03 +0.000000e+00 1.582045e+00 2.940206e-03 +0.000000e+00 1.603534e+00 3.074132e-03 +0.000000e+00 1.625315e+00 3.211721e-03 +0.000000e+00 1.647392e+00 3.353073e-03 +0.000000e+00 1.669768e+00 3.498291e-03 +0.000000e+00 1.692449e+00 3.647481e-03 +0.000000e+00 1.715438e+00 3.800751e-03 +0.000000e+00 1.738739e+00 3.958214e-03 +0.000000e+00 1.762356e+00 4.119983e-03 +0.000000e+00 1.786295e+00 4.286177e-03 +0.000000e+00 1.810558e+00 4.456916e-03 +0.000000e+00 1.835151e+00 4.632326e-03 +0.000000e+00 1.860078e+00 4.812532e-03 +0.000000e+00 1.885344e+00 4.997668e-03 +0.000000e+00 1.910953e+00 5.187866e-03 +0.000000e+00 1.936910e+00 5.383268e-03 +0.000000e+00 1.963219e+00 5.584014e-03 +0.000000e+00 1.989886e+00 5.790249e-03 +0.000000e+00 2.016915e+00 5.364985e-03 +0.000000e+00 2.044311e+00 5.489168e-03 +0.000000e+00 2.072079e+00 6.160623e-03 +0.000000e+00 2.100224e+00 7.327441e-03 +0.000000e+00 2.128752e+00 8.941695e-03 +0.000000e+00 2.157667e+00 1.095899e-02 +0.000000e+00 2.186975e+00 1.333848e-02 +0.000000e+00 2.216681e+00 1.604234e-02 +0.000000e+00 2.246790e+00 1.903586e-02 +0.000000e+00 2.277308e+00 2.228691e-02 +0.000000e+00 2.308242e+00 2.576607e-02 +0.000000e+00 2.339595e+00 2.944617e-02 +0.000000e+00 2.371374e+00 3.330239e-02 +0.000000e+00 2.403584e+00 3.731188e-02 +0.000000e+00 2.436233e+00 4.145370e-02 +0.000000e+00 2.469324e+00 4.570878e-02 +0.000000e+00 2.502865e+00 5.005958e-02 +0.000000e+00 2.536862e+00 5.449019e-02 +0.000000e+00 2.571321e+00 5.898609e-02 +0.000000e+00 2.606247e+00 6.353403e-02 +0.000000e+00 2.641648e+00 6.812201e-02 +0.000000e+00 2.677530e+00 7.273922e-02 +0.000000e+00 2.713899e+00 7.737575e-02 +0.000000e+00 2.750763e+00 8.202283e-02 +0.000000e+00 2.788127e+00 8.667243e-02 +0.000000e+00 2.825998e+00 9.131740e-02 +0.000000e+00 2.864384e+00 9.595138e-02 +0.000000e+00 2.903291e+00 1.005686e-01 +0.000000e+00 2.942727e+00 1.051641e-01 +0.000000e+00 2.982699e+00 1.097333e-01 +0.000000e+00 3.023213e+00 1.142723e-01 +0.000000e+00 3.064278e+00 1.187777e-01 +0.000000e+00 3.105900e+00 1.232465e-01 +0.000000e+00 3.148088e+00 1.276762e-01 +0.000000e+00 3.190849e+00 1.320645e-01 +0.000000e+00 3.234191e+00 1.364096e-01 +0.000000e+00 3.278121e+00 1.407099e-01 +0.000000e+00 3.322648e+00 1.449643e-01 +0.000000e+00 3.367780e+00 1.491717e-01 +0.000000e+00 3.413525e+00 1.533314e-01 +0.000000e+00 3.459892e+00 1.574428e-01 +0.000000e+00 3.506888e+00 1.615057e-01 +0.000000e+00 3.554522e+00 1.655198e-01 +0.000000e+00 3.602804e+00 1.694851e-01 +0.000000e+00 3.651741e+00 1.734016e-01 +0.000000e+00 3.701343e+00 1.772698e-01 +0.000000e+00 3.751619e+00 1.810898e-01 +0.000000e+00 3.802578e+00 1.848621e-01 +0.000000e+00 3.854229e+00 1.885873e-01 +0.000000e+00 3.906581e+00 1.922659e-01 +0.000000e+00 3.959645e+00 1.958985e-01 +0.000000e+00 4.013429e+00 1.994860e-01 +0.000000e+00 4.067945e+00 2.030290e-01 +0.000000e+00 4.123200e+00 2.065282e-01 +0.000000e+00 4.179206e+00 2.099846e-01 +0.000000e+00 4.235972e+00 2.133990e-01 +0.000000e+00 4.293510e+00 2.167722e-01 +0.000000e+00 4.351830e+00 2.201050e-01 +0.000000e+00 4.410941e+00 2.233984e-01 +0.000000e+00 4.470855e+00 2.266532e-01 +0.000000e+00 4.531584e+00 2.298704e-01 +0.000000e+00 4.593137e+00 2.330508e-01 +0.000000e+00 4.655526e+00 2.361953e-01 +0.000000e+00 4.718762e+00 2.393046e-01 +0.000000e+00 4.782858e+00 2.423798e-01 +0.000000e+00 4.847825e+00 2.454216e-01 +0.000000e+00 4.913673e+00 2.484308e-01 +0.000000e+00 4.980416e+00 2.514084e-01 +0.000000e+00 5.048066e+00 2.543549e-01 +0.000000e+00 5.116634e+00 2.572713e-01 +0.000000e+00 5.186134e+00 2.601582e-01 +0.000000e+00 5.256578e+00 2.630165e-01 +0.000000e+00 5.327979e+00 2.658468e-01 +0.000000e+00 5.400350e+00 2.686498e-01 +0.000000e+00 5.473703e+00 2.714261e-01 +0.000000e+00 5.548053e+00 2.741765e-01 +0.000000e+00 5.623413e+00 2.769015e-01 +0.000000e+00 5.699797e+00 2.796018e-01 +0.000000e+00 5.777218e+00 2.822779e-01 +0.000000e+00 5.855690e+00 2.849304e-01 +0.000000e+00 5.935229e+00 2.875599e-01 +0.000000e+00 6.015848e+00 2.901667e-01 +0.000000e+00 6.097562e+00 2.927515e-01 +0.000000e+00 6.180387e+00 2.953147e-01 +0.000000e+00 6.264335e+00 2.978567e-01 +0.000000e+00 6.349425e+00 3.003781e-01 +0.000000e+00 6.435670e+00 3.028791e-01 +0.000000e+00 6.523087e+00 3.053602e-01 +0.000000e+00 6.611690e+00 3.078218e-01 +0.000000e+00 6.701498e+00 3.102642e-01 +0.000000e+00 6.792525e+00 3.126878e-01 +0.000000e+00 6.884789e+00 3.150929e-01 +0.000000e+00 6.978306e+00 3.174797e-01 +0.000000e+00 7.073093e+00 3.198487e-01 +0.000000e+00 7.169168e+00 3.222000e-01 +0.000000e+00 7.266548e+00 3.245339e-01 +0.000000e+00 7.365250e+00 3.268507e-01 +0.000000e+00 7.465293e+00 3.291506e-01 +0.000000e+00 7.566695e+00 3.314338e-01 +0.000000e+00 7.669475e+00 3.337005e-01 +0.000000e+00 7.773650e+00 3.359510e-01 +0.000000e+00 7.879241e+00 3.381854e-01 +0.000000e+00 7.986266e+00 3.404039e-01 +0.000000e+00 8.094744e+00 3.426065e-01 +0.000000e+00 8.204696e+00 3.447936e-01 +0.000000e+00 8.316141e+00 3.469653e-01 +0.000000e+00 8.429101e+00 3.491216e-01 +0.000000e+00 8.543594e+00 3.512627e-01 +0.000000e+00 8.659643e+00 3.533887e-01 +0.000000e+00 8.777268e+00 3.554997e-01 +0.000000e+00 8.896491e+00 3.575959e-01 +0.000000e+00 9.017333e+00 3.596773e-01 +0.000000e+00 9.139817e+00 3.617441e-01 +0.000000e+00 9.263965e+00 3.637963e-01 +0.000000e+00 9.389798e+00 3.658339e-01 +0.000000e+00 9.517341e+00 3.678571e-01 +0.000000e+00 9.646616e+00 3.698659e-01 +0.000000e+00 9.777647e+00 3.718604e-01 +0.000000e+00 9.910459e+00 3.738407e-01 +0.000000e+00 1.004507e+01 3.758068e-01 +0.000000e+00 1.018152e+01 3.777588e-01 +0.000000e+00 1.031981e+01 3.796967e-01 +0.000000e+00 1.045999e+01 3.816206e-01 +0.000000e+00 1.060207e+01 3.835305e-01 +0.000000e+00 1.074608e+01 3.854264e-01 +0.000000e+00 1.089204e+01 3.873085e-01 +0.000000e+00 1.103999e+01 3.891766e-01 +0.000000e+00 1.118995e+01 3.910310e-01 +0.000000e+00 1.134194e+01 3.928715e-01 +0.000000e+00 1.149600e+01 3.946983e-01 +0.000000e+00 1.165215e+01 3.965114e-01 +0.000000e+00 1.181043e+01 3.983108e-01 +0.000000e+00 1.197085e+01 4.000965e-01 +0.000000e+00 1.213345e+01 4.018686e-01 +0.000000e+00 1.229826e+01 4.036270e-01 +0.000000e+00 1.246531e+01 4.053719e-01 +0.000000e+00 1.263463e+01 4.071033e-01 +0.000000e+00 1.280625e+01 4.088211e-01 +0.000000e+00 1.298020e+01 4.105255e-01 +0.000000e+00 1.315651e+01 4.122163e-01 +0.000000e+00 1.333521e+01 4.138938e-01 +0.000000e+00 1.351635e+01 4.155579e-01 +0.000000e+00 1.369994e+01 4.172085e-01 +0.000000e+00 1.388603e+01 4.188459e-01 +0.000000e+00 1.407465e+01 4.204699e-01 +0.000000e+00 1.426582e+01 4.220806e-01 +0.000000e+00 1.445960e+01 4.236781e-01 +0.000000e+00 1.465601e+01 4.252624e-01 +0.000000e+00 1.485508e+01 4.268335e-01 +0.000000e+00 1.505686e+01 4.283915e-01 +0.000000e+00 1.526138e+01 4.299363e-01 +0.000000e+00 1.546868e+01 4.314680e-01 +0.000000e+00 1.567879e+01 4.329867e-01 +0.000000e+00 1.589176e+01 4.344924e-01 +0.000000e+00 1.610761e+01 4.359852e-01 +0.000000e+00 1.632641e+01 4.374650e-01 +0.000000e+00 1.654817e+01 4.389319e-01 +0.000000e+00 1.677295e+01 4.403860e-01 +0.000000e+00 1.700078e+01 4.418273e-01 +0.000000e+00 1.723170e+01 4.432558e-01 +0.000000e+00 1.746576e+01 4.446716e-01 +0.000000e+00 1.770300e+01 4.460747e-01 +0.000000e+00 1.794346e+01 4.474652e-01 +0.000000e+00 1.818719e+01 4.488431e-01 +0.000000e+00 1.843423e+01 4.502085e-01 +0.000000e+00 1.868462e+01 4.515614e-01 +0.000000e+00 1.893842e+01 4.529019e-01 +0.000000e+00 1.919566e+01 4.542300e-01 +0.000000e+00 1.945640e+01 4.555458e-01 +0.000000e+00 1.972068e+01 4.568494e-01 +0.000000e+00 1.998855e+01 4.581407e-01 +0.000000e+00 2.026006e+01 4.594198e-01 +0.000000e+00 2.053525e+01 4.606869e-01 +0.000000e+00 2.081418e+01 4.619419e-01 +0.000000e+00 2.109690e+01 4.631850e-01 +0.000000e+00 2.138347e+01 4.644160e-01 +0.000000e+00 2.167392e+01 4.656353e-01 +0.000000e+00 2.196832e+01 4.668427e-01 +0.000000e+00 2.226672e+01 4.680384e-01 +0.000000e+00 2.256917e+01 4.692224e-01 +0.000000e+00 2.287573e+01 4.703948e-01 +0.000000e+00 2.318646e+01 4.715557e-01 +0.000000e+00 2.350140e+01 4.727050e-01 +0.000000e+00 2.382062e+01 4.738429e-01 +0.000000e+00 2.414418e+01 4.749695e-01 +0.000000e+00 2.447214e+01 4.760848e-01 +0.000000e+00 2.480454e+01 4.771889e-01 +0.000000e+00 2.514147e+01 4.782818e-01 +0.000000e+00 2.548297e+01 4.793636e-01 +0.000000e+00 2.582911e+01 4.804344e-01 +0.000000e+00 2.617995e+01 4.814944e-01 +0.000000e+00 2.653555e+01 4.825434e-01 +0.000000e+00 2.689599e+01 4.835816e-01 +0.000000e+00 2.726132e+01 4.846091e-01 +0.000000e+00 2.763161e+01 4.856260e-01 +0.000000e+00 2.800694e+01 4.866322e-01 +0.000000e+00 2.838736e+01 4.876280e-01 +0.000000e+00 2.877295e+01 4.886133e-01 +0.000000e+00 2.916378e+01 4.895883e-01 +0.000000e+00 2.955991e+01 4.905529e-01 +0.000000e+00 2.996143e+01 4.915074e-01 +0.000000e+00 3.036840e+01 4.924517e-01 +0.000000e+00 3.078090e+01 4.933860e-01 +0.000000e+00 3.119900e+01 4.943102e-01 +0.000000e+00 3.162278e+01 4.952246e-01 +0.000000e+00 3.205231e+01 4.961291e-01 +0.000000e+00 3.248768e+01 4.970239e-01 +0.000000e+00 3.292897e+01 4.979090e-01 +0.000000e+00 3.337625e+01 4.987845e-01 +0.000000e+00 3.382960e+01 4.996504e-01 +0.000000e+00 3.428911e+01 5.005069e-01 +0.000000e+00 3.475487e+01 5.013540e-01 +0.000000e+00 3.522695e+01 5.021918e-01 +0.000000e+00 3.570544e+01 5.030204e-01 +0.000000e+00 3.619043e+01 5.038398e-01 +0.000000e+00 3.668201e+01 5.046502e-01 +0.000000e+00 3.718027e+01 5.054516e-01 +0.000000e+00 3.768529e+01 5.062440e-01 +0.000000e+00 3.819717e+01 5.070276e-01 +0.000000e+00 3.871601e+01 5.078024e-01 +0.000000e+00 3.924190e+01 5.085685e-01 +0.000000e+00 3.977493e+01 5.093261e-01 +0.000000e+00 4.031519e+01 5.100751e-01 +0.000000e+00 4.086280e+01 5.108156e-01 +0.000000e+00 4.141785e+01 5.115477e-01 +0.000000e+00 4.198043e+01 5.122716e-01 +0.000000e+00 4.255066e+01 5.129871e-01 +0.000000e+00 4.312863e+01 5.136946e-01 +0.000000e+00 4.371445e+01 5.143939e-01 +0.000000e+00 4.430823e+01 5.150853e-01 +0.000000e+00 4.491007e+01 5.157686e-01 +0.000000e+00 4.552009e+01 5.164442e-01 +0.000000e+00 4.613840e+01 5.171119e-01 +0.000000e+00 4.676510e+01 5.177720e-01 +0.000000e+00 4.740032e+01 5.184243e-01 +0.000000e+00 4.804416e+01 5.190692e-01 +0.000000e+00 4.869675e+01 5.197065e-01 +0.000000e+00 4.935821e+01 5.203364e-01 +0.000000e+00 5.002864e+01 5.209590e-01 +0.000000e+00 5.070819e+01 5.215742e-01 +0.000000e+00 5.139697e+01 5.221822e-01 +0.000000e+00 5.209510e+01 5.227832e-01 +0.000000e+00 5.280272e+01 5.233771e-01 +0.000000e+00 5.351994e+01 5.239639e-01 +0.000000e+00 5.424691e+01 5.245438e-01 +0.000000e+00 5.498375e+01 5.251169e-01 +0.000000e+00 5.573060e+01 5.256832e-01 +0.000000e+00 5.648760e+01 5.262427e-01 +0.000000e+00 5.725488e+01 5.267957e-01 +0.000000e+00 5.803258e+01 5.273421e-01 +0.000000e+00 5.882084e+01 5.278819e-01 +0.000000e+00 5.961982e+01 5.284153e-01 +0.000000e+00 6.042964e+01 5.289423e-01 +0.000000e+00 6.125046e+01 5.294629e-01 +0.000000e+00 6.208244e+01 5.299774e-01 +0.000000e+00 6.292571e+01 5.304856e-01 +0.000000e+00 6.378044e+01 5.309878e-01 +0.000000e+00 6.464677e+01 5.314839e-01 +0.000000e+00 6.552488e+01 5.319740e-01 +0.000000e+00 6.641492e+01 5.324582e-01 +0.000000e+00 6.731704e+01 5.329365e-01 +0.000000e+00 6.823141e+01 5.334090e-01 +0.000000e+00 6.915821e+01 5.338758e-01 +0.000000e+00 7.009760e+01 5.343369e-01 +0.000000e+00 7.104974e+01 5.347924e-01 +0.000000e+00 7.201482e+01 5.352423e-01 +0.000000e+00 7.299300e+01 5.356868e-01 +0.000000e+00 7.398448e+01 5.361258e-01 +0.000000e+00 7.498942e+01 5.365595e-01 +0.000000e+00 7.600801e+01 5.369878e-01 +0.000000e+00 7.704044e+01 5.374109e-01 +0.000000e+00 7.808689e+01 5.378288e-01 +0.000000e+00 7.914755e+01 5.382416e-01 +0.000000e+00 8.022263e+01 5.386492e-01 +0.000000e+00 8.131230e+01 5.390519e-01 +0.000000e+00 8.241678e+01 5.394495e-01 +0.000000e+00 8.353625e+01 5.398423e-01 +0.000000e+00 8.467094e+01 5.402302e-01 +0.000000e+00 8.582104e+01 5.406133e-01 +0.000000e+00 8.698676e+01 5.409917e-01 +0.000000e+00 8.816830e+01 5.413653e-01 +0.000000e+00 8.936591e+01 5.417343e-01 +0.000000e+00 9.057978e+01 5.420988e-01 +0.000000e+00 9.181013e+01 5.424587e-01 +0.000000e+00 9.305721e+01 5.428141e-01 +0.000000e+00 9.432121e+01 5.431651e-01 +0.000000e+00 9.560239e+01 5.435117e-01 +0.000000e+00 9.690097e+01 5.438540e-01 +0.000000e+00 9.821719e+01 5.441920e-01 +0.000000e+00 9.955128e+01 5.445257e-01 +0.000000e+00 1.009035e+02 5.446360e-01 +0.000000e+00 1.022741e+02 5.446360e-01 +0.000000e+00 1.036633e+02 5.446360e-01 +0.000000e+00 1.050714e+02 5.446360e-01 +0.000000e+00 1.064986e+02 5.446360e-01 +0.000000e+00 1.079451e+02 5.446360e-01 +0.000000e+00 1.094114e+02 5.446360e-01 +0.000000e+00 1.108975e+02 5.446360e-01 +0.000000e+00 1.124039e+02 5.446360e-01 +0.000000e+00 1.139307e+02 5.446360e-01 +0.000000e+00 1.154782e+02 5.446360e-01 +0.000000e+00 1.170468e+02 5.446360e-01 +0.000000e+00 1.186366e+02 5.446360e-01 +0.000000e+00 1.202481e+02 5.446360e-01 +0.000000e+00 1.218814e+02 5.446360e-01 +0.000000e+00 1.235369e+02 5.446360e-01 +0.000000e+00 1.252150e+02 5.446360e-01 +0.000000e+00 1.269158e+02 5.446360e-01 +0.000000e+00 1.286397e+02 5.446360e-01 +0.000000e+00 1.303870e+02 5.446360e-01 +0.000000e+00 1.321581e+02 5.446360e-01 +0.000000e+00 1.339532e+02 5.446360e-01 +0.000000e+00 1.357727e+02 5.446360e-01 +0.000000e+00 1.376169e+02 5.446360e-01 +0.000000e+00 1.394862e+02 5.446360e-01 +0.000000e+00 1.413809e+02 5.446360e-01 +0.000000e+00 1.433013e+02 5.446360e-01 +0.000000e+00 1.452477e+02 5.446360e-01 +0.000000e+00 1.472207e+02 5.446360e-01 +0.000000e+00 1.492204e+02 5.446360e-01 +0.000000e+00 1.512473e+02 5.446360e-01 +0.000000e+00 1.533017e+02 5.446360e-01 +0.000000e+00 1.553840e+02 5.446360e-01 +0.000000e+00 1.574946e+02 5.446360e-01 +0.000000e+00 1.596339e+02 5.446360e-01 +0.000000e+00 1.618022e+02 5.446360e-01 +0.000000e+00 1.640000e+02 5.446360e-01 +0.000000e+00 1.662276e+02 5.446360e-01 +0.000000e+00 1.684855e+02 5.446360e-01 +0.000000e+00 1.707740e+02 5.446360e-01 +0.000000e+00 1.730937e+02 5.446360e-01 +0.000000e+00 1.754449e+02 5.446360e-01 +0.000000e+00 1.778279e+02 5.446360e-01 +0.000000e+00 1.802434e+02 5.446360e-01 +0.000000e+00 1.826917e+02 5.446360e-01 +0.000000e+00 1.851732e+02 5.446360e-01 +0.000000e+00 1.876884e+02 5.446360e-01 +0.000000e+00 1.902378e+02 5.446360e-01 +0.000000e+00 1.928219e+02 5.446360e-01 +0.000000e+00 1.954410e+02 5.446360e-01 +0.000000e+00 1.980957e+02 5.446360e-01 +0.000000e+00 2.007864e+02 5.446360e-01 +0.000000e+00 2.035137e+02 5.446360e-01 +0.000000e+00 2.062781e+02 5.446360e-01 +0.000000e+00 2.090800e+02 5.446360e-01 +0.000000e+00 2.119200e+02 5.446360e-01 +0.000000e+00 2.147985e+02 5.446360e-01 +0.000000e+00 2.177161e+02 5.446360e-01 +0.000000e+00 2.206734e+02 5.446360e-01 +0.000000e+00 2.236708e+02 5.446360e-01 +0.000000e+00 2.267090e+02 5.446360e-01 +0.000000e+00 2.297884e+02 5.446360e-01 +0.000000e+00 2.329097e+02 5.446360e-01 +0.000000e+00 2.360733e+02 5.446360e-01 +0.000000e+00 2.392799e+02 5.446360e-01 +0.000000e+00 2.425301e+02 5.446360e-01 +0.000000e+00 2.458244e+02 5.446360e-01 +0.000000e+00 2.491635e+02 5.446360e-01 +0.000000e+00 2.525479e+02 5.446360e-01 +0.000000e+00 2.559783e+02 5.446360e-01 +0.000000e+00 2.594553e+02 5.446360e-01 +0.000000e+00 2.629795e+02 5.446360e-01 +0.000000e+00 2.665516e+02 5.446360e-01 +0.000000e+00 2.701722e+02 5.446360e-01 +0.000000e+00 2.738419e+02 5.446360e-01 +0.000000e+00 2.775616e+02 5.446360e-01 +0.000000e+00 2.813318e+02 5.446360e-01 +0.000000e+00 2.851531e+02 5.446360e-01 +0.000000e+00 2.890264e+02 5.446360e-01 +0.000000e+00 2.929523e+02 5.446360e-01 +0.000000e+00 2.969315e+02 5.446360e-01 +0.000000e+00 3.009648e+02 5.446360e-01 +0.000000e+00 3.050528e+02 5.446360e-01 +0.000000e+00 3.091964e+02 5.446360e-01 +0.000000e+00 3.133962e+02 5.446360e-01 +0.000000e+00 3.176531e+02 5.446360e-01 +0.000000e+00 3.219678e+02 5.446360e-01 +0.000000e+00 3.263412e+02 5.446360e-01 +0.000000e+00 3.307739e+02 5.446360e-01 +0.000000e+00 3.352668e+02 5.446360e-01 +0.000000e+00 3.398208e+02 5.446360e-01 +0.000000e+00 3.444367e+02 5.446360e-01 +0.000000e+00 3.491152e+02 5.446360e-01 +0.000000e+00 3.538573e+02 5.446360e-01 +0.000000e+00 3.586638e+02 5.446360e-01 +0.000000e+00 3.635356e+02 5.446360e-01 +0.000000e+00 3.684735e+02 5.446360e-01 +0.000000e+00 3.734785e+02 5.446360e-01 +0.000000e+00 3.785515e+02 5.446360e-01 +0.000000e+00 3.836935e+02 5.446360e-01 +0.000000e+00 3.889052e+02 5.446360e-01 +0.000000e+00 3.941877e+02 5.446360e-01 +0.000000e+00 3.995421e+02 5.446360e-01 +0.000000e+00 4.049691e+02 5.446360e-01 +0.000000e+00 4.104698e+02 5.446360e-01 +0.000000e+00 4.160453e+02 5.446360e-01 +0.000000e+00 4.216965e+02 5.446360e-01 +0.000000e+00 4.274245e+02 5.446360e-01 +0.000000e+00 4.332302e+02 5.446360e-01 +0.000000e+00 4.391148e+02 5.446360e-01 +0.000000e+00 4.450794e+02 5.446360e-01 +0.000000e+00 4.511250e+02 5.446360e-01 +0.000000e+00 4.572527e+02 5.446360e-01 +0.000000e+00 4.634636e+02 5.446360e-01 +0.000000e+00 4.697589e+02 5.446360e-01 +0.000000e+00 4.761397e+02 5.446360e-01 +0.000000e+00 4.826071e+02 5.446360e-01 +0.000000e+00 4.891625e+02 5.446360e-01 +0.000000e+00 4.958068e+02 5.446360e-01 +0.000000e+00 5.025414e+02 5.446360e-01 +0.000000e+00 5.093675e+02 5.446360e-01 +0.000000e+00 5.162863e+02 5.446360e-01 +0.000000e+00 5.232991e+02 5.446360e-01 +0.000000e+00 5.304072e+02 5.446360e-01 +0.000000e+00 5.376118e+02 5.446360e-01 +0.000000e+00 5.449142e+02 5.446360e-01 +0.000000e+00 5.523159e+02 5.446360e-01 +0.000000e+00 5.598180e+02 5.446360e-01 +0.000000e+00 5.674221e+02 5.446360e-01 +0.000000e+00 5.751295e+02 5.446360e-01 +0.000000e+00 5.829415e+02 5.446360e-01 +0.000000e+00 5.908597e+02 5.446360e-01 +0.000000e+00 5.988854e+02 5.446360e-01 +0.000000e+00 6.070202e+02 5.446360e-01 +0.000000e+00 6.152654e+02 5.446360e-01 +0.000000e+00 6.236226e+02 5.446360e-01 +0.000000e+00 6.320934e+02 5.446360e-01 +0.000000e+00 6.406792e+02 5.446360e-01 +0.000000e+00 6.493817e+02 5.446360e-01 +0.000000e+00 6.582023e+02 5.446360e-01 +0.000000e+00 6.671427e+02 5.446360e-01 +0.000000e+00 6.762046e+02 5.446360e-01 +0.000000e+00 6.853896e+02 5.446360e-01 +0.000000e+00 6.946993e+02 5.446360e-01 +0.000000e+00 7.041355e+02 5.446360e-01 +0.000000e+00 7.136999e+02 5.446360e-01 +0.000000e+00 7.233942e+02 5.446360e-01 +0.000000e+00 7.332201e+02 5.446360e-01 +0.000000e+00 7.431796e+02 5.446360e-01 +0.000000e+00 7.532742e+02 5.446360e-01 +0.000000e+00 7.635061e+02 5.446360e-01 +0.000000e+00 7.738769e+02 5.446360e-01 +0.000000e+00 7.843885e+02 5.446360e-01 +0.000000e+00 7.950430e+02 5.446360e-01 +0.000000e+00 8.058422e+02 5.446360e-01 +0.000000e+00 8.167880e+02 5.446360e-01 +0.000000e+00 8.278826e+02 5.446360e-01 +0.000000e+00 8.391278e+02 5.446360e-01 +0.000000e+00 8.505258e+02 5.446360e-01 +0.000000e+00 8.620786e+02 5.446360e-01 +0.000000e+00 8.737883e+02 5.446360e-01 +0.000000e+00 8.856571e+02 5.446360e-01 +0.000000e+00 8.976871e+02 5.446360e-01 +0.000000e+00 9.098806e+02 5.446360e-01 +0.000000e+00 9.222396e+02 5.446360e-01 +0.000000e+00 9.347665e+02 5.446360e-01 +0.000000e+00 9.474635e+02 5.446360e-01 +0.000000e+00 9.603331e+02 5.446360e-01 +0.000000e+00 9.733774e+02 5.446360e-01 +0.000000e+00 9.865989e+02 5.446360e-01 +3.000000e-01 1.000000e+00 0.000000e+00 +3.000000e-01 1.013583e+00 3.750115e-06 +3.000000e-01 1.027351e+00 2.883254e-05 +3.000000e-01 1.041305e+00 9.515434e-05 +3.000000e-01 1.055450e+00 2.200341e-04 +3.000000e-01 1.069786e+00 4.196625e-04 +3.000000e-01 1.084317e+00 7.085034e-04 +3.000000e-01 1.099045e+00 1.099257e-03 +3.000000e-01 1.113974e+00 1.603124e-03 +3.000000e-01 1.129105e+00 2.229906e-03 +3.000000e-01 1.144442e+00 2.988489e-03 +3.000000e-01 1.159987e+00 3.886489e-03 +3.000000e-01 1.175743e+00 4.930094e-03 +3.000000e-01 1.191714e+00 6.125025e-03 +3.000000e-01 1.207901e+00 7.475133e-03 +3.000000e-01 1.224308e+00 8.984170e-03 +3.000000e-01 1.240938e+00 1.065498e-02 +3.000000e-01 1.257794e+00 1.248916e-02 +3.000000e-01 1.274878e+00 1.448819e-02 +3.000000e-01 1.292195e+00 1.665242e-02 +3.000000e-01 1.309747e+00 1.898190e-02 +3.000000e-01 1.327538e+00 2.147606e-02 +3.000000e-01 1.345570e+00 2.413343e-02 +3.000000e-01 1.363847e+00 2.695237e-02 +3.000000e-01 1.382372e+00 2.993114e-02 +3.000000e-01 1.401149e+00 3.306693e-02 +3.000000e-01 1.420181e+00 3.635694e-02 +3.000000e-01 1.439472e+00 3.979807e-02 +3.000000e-01 1.459024e+00 4.338667e-02 +3.000000e-01 1.478842e+00 4.711908e-02 +3.000000e-01 1.498930e+00 5.099142e-02 +3.000000e-01 1.519290e+00 5.584033e-02 +3.000000e-01 1.539927e+00 6.085568e-02 +3.000000e-01 1.560844e+00 6.596592e-02 +3.000000e-01 1.582045e+00 7.114946e-02 +3.000000e-01 1.603534e+00 7.638680e-02 +3.000000e-01 1.625315e+00 8.166029e-02 +3.000000e-01 1.647392e+00 8.695409e-02 +3.000000e-01 1.669768e+00 9.225400e-02 +3.000000e-01 1.692449e+00 9.754725e-02 +3.000000e-01 1.715438e+00 1.028225e-01 +3.000000e-01 1.738739e+00 1.080698e-01 +3.000000e-01 1.762356e+00 1.132801e-01 +3.000000e-01 1.786295e+00 1.184456e-01 +3.000000e-01 1.810558e+00 1.235595e-01 +3.000000e-01 1.835151e+00 1.286159e-01 +3.000000e-01 1.860078e+00 1.336097e-01 +3.000000e-01 1.885344e+00 1.385366e-01 +3.000000e-01 1.910953e+00 1.433927e-01 +3.000000e-01 1.936910e+00 1.481753e-01 +3.000000e-01 1.963219e+00 1.528819e-01 +3.000000e-01 1.989886e+00 1.575104e-01 +3.000000e-01 2.016915e+00 1.620595e-01 +3.000000e-01 2.044311e+00 1.665282e-01 +3.000000e-01 2.072079e+00 1.709157e-01 +3.000000e-01 2.100224e+00 1.752219e-01 +3.000000e-01 2.128752e+00 1.794467e-01 +3.000000e-01 2.157667e+00 1.835904e-01 +3.000000e-01 2.186975e+00 1.876535e-01 +3.000000e-01 2.216681e+00 1.916369e-01 +3.000000e-01 2.246790e+00 1.955415e-01 +3.000000e-01 2.277308e+00 1.993683e-01 +3.000000e-01 2.308242e+00 2.031187e-01 +3.000000e-01 2.339595e+00 2.067941e-01 +3.000000e-01 2.371374e+00 2.103959e-01 +3.000000e-01 2.403584e+00 2.139257e-01 +3.000000e-01 2.436233e+00 2.173853e-01 +3.000000e-01 2.469324e+00 2.207763e-01 +3.000000e-01 2.502865e+00 2.241005e-01 +3.000000e-01 2.536862e+00 2.273598e-01 +3.000000e-01 2.571321e+00 2.305560e-01 +3.000000e-01 2.606247e+00 2.336910e-01 +3.000000e-01 2.641648e+00 2.367665e-01 +3.000000e-01 2.677530e+00 2.397846e-01 +3.000000e-01 2.713899e+00 2.427470e-01 +3.000000e-01 2.750763e+00 2.456557e-01 +3.000000e-01 2.788127e+00 2.485124e-01 +3.000000e-01 2.825998e+00 2.513190e-01 +3.000000e-01 2.864384e+00 2.540773e-01 +3.000000e-01 2.903291e+00 2.567890e-01 +3.000000e-01 2.942727e+00 2.594558e-01 +3.000000e-01 2.982699e+00 2.620796e-01 +3.000000e-01 3.023213e+00 2.646618e-01 +3.000000e-01 3.064278e+00 2.672042e-01 +3.000000e-01 3.105900e+00 2.697082e-01 +3.000000e-01 3.148088e+00 2.721755e-01 +3.000000e-01 3.190849e+00 2.746074e-01 +3.000000e-01 3.234191e+00 2.770055e-01 +3.000000e-01 3.278121e+00 2.793711e-01 +3.000000e-01 3.322648e+00 2.817056e-01 +3.000000e-01 3.367780e+00 2.840101e-01 +3.000000e-01 3.413525e+00 2.862862e-01 +3.000000e-01 3.459892e+00 2.885348e-01 +3.000000e-01 3.506888e+00 2.907572e-01 +3.000000e-01 3.554522e+00 2.929544e-01 +3.000000e-01 3.602804e+00 2.951275e-01 +3.000000e-01 3.651741e+00 2.972776e-01 +3.000000e-01 3.701343e+00 2.994056e-01 +3.000000e-01 3.751619e+00 3.015125e-01 +3.000000e-01 3.802578e+00 3.035991e-01 +3.000000e-01 3.854229e+00 3.056662e-01 +3.000000e-01 3.906581e+00 3.077148e-01 +3.000000e-01 3.959645e+00 3.097455e-01 +3.000000e-01 4.013429e+00 3.117590e-01 +3.000000e-01 4.067945e+00 3.137561e-01 +3.000000e-01 4.123200e+00 3.157374e-01 +3.000000e-01 4.179206e+00 3.177036e-01 +3.000000e-01 4.235972e+00 3.196552e-01 +3.000000e-01 4.293510e+00 3.215926e-01 +3.000000e-01 4.351830e+00 3.235167e-01 +3.000000e-01 4.410941e+00 3.254276e-01 +3.000000e-01 4.470855e+00 3.273259e-01 +3.000000e-01 4.531584e+00 3.292121e-01 +3.000000e-01 4.593137e+00 3.310865e-01 +3.000000e-01 4.655526e+00 3.329495e-01 +3.000000e-01 4.718762e+00 3.348015e-01 +3.000000e-01 4.782858e+00 3.366428e-01 +3.000000e-01 4.847825e+00 3.384736e-01 +3.000000e-01 4.913673e+00 3.402943e-01 +3.000000e-01 4.980416e+00 3.421051e-01 +3.000000e-01 5.048066e+00 3.439062e-01 +3.000000e-01 5.116634e+00 3.456979e-01 +3.000000e-01 5.186134e+00 3.474804e-01 +3.000000e-01 5.256578e+00 3.492538e-01 +3.000000e-01 5.327979e+00 3.510182e-01 +3.000000e-01 5.400350e+00 3.527740e-01 +3.000000e-01 5.473703e+00 3.545210e-01 +3.000000e-01 5.548053e+00 3.562596e-01 +3.000000e-01 5.623413e+00 3.579898e-01 +3.000000e-01 5.699797e+00 3.597116e-01 +3.000000e-01 5.777218e+00 3.614253e-01 +3.000000e-01 5.855690e+00 3.631307e-01 +3.000000e-01 5.935229e+00 3.648280e-01 +3.000000e-01 6.015848e+00 3.665173e-01 +3.000000e-01 6.097562e+00 3.681986e-01 +3.000000e-01 6.180387e+00 3.698718e-01 +3.000000e-01 6.264335e+00 3.715371e-01 +3.000000e-01 6.349425e+00 3.731944e-01 +3.000000e-01 6.435670e+00 3.748437e-01 +3.000000e-01 6.523087e+00 3.764851e-01 +3.000000e-01 6.611690e+00 3.781185e-01 +3.000000e-01 6.701498e+00 3.797440e-01 +3.000000e-01 6.792525e+00 3.813615e-01 +3.000000e-01 6.884789e+00 3.829710e-01 +3.000000e-01 6.978306e+00 3.845724e-01 +3.000000e-01 7.073093e+00 3.861657e-01 +3.000000e-01 7.169168e+00 3.877510e-01 +3.000000e-01 7.266548e+00 3.893282e-01 +3.000000e-01 7.365250e+00 3.908971e-01 +3.000000e-01 7.465293e+00 3.924579e-01 +3.000000e-01 7.566695e+00 3.940104e-01 +3.000000e-01 7.669475e+00 3.955546e-01 +3.000000e-01 7.773650e+00 3.970905e-01 +3.000000e-01 7.879241e+00 3.986180e-01 +3.000000e-01 7.986266e+00 4.001371e-01 +3.000000e-01 8.094744e+00 4.016476e-01 +3.000000e-01 8.204696e+00 4.031497e-01 +3.000000e-01 8.316141e+00 4.046431e-01 +3.000000e-01 8.429101e+00 4.061280e-01 +3.000000e-01 8.543594e+00 4.076042e-01 +3.000000e-01 8.659643e+00 4.090717e-01 +3.000000e-01 8.777268e+00 4.105303e-01 +3.000000e-01 8.896491e+00 4.119802e-01 +3.000000e-01 9.017333e+00 4.134211e-01 +3.000000e-01 9.139817e+00 4.148532e-01 +3.000000e-01 9.263965e+00 4.162763e-01 +3.000000e-01 9.389798e+00 4.176904e-01 +3.000000e-01 9.517341e+00 4.190955e-01 +3.000000e-01 9.646616e+00 4.204914e-01 +3.000000e-01 9.777647e+00 4.218783e-01 +3.000000e-01 9.910459e+00 4.232559e-01 +3.000000e-01 1.004507e+01 4.246244e-01 +3.000000e-01 1.018152e+01 4.259836e-01 +3.000000e-01 1.031981e+01 4.273335e-01 +3.000000e-01 1.045999e+01 4.286741e-01 +3.000000e-01 1.060207e+01 4.300054e-01 +3.000000e-01 1.074608e+01 4.313273e-01 +3.000000e-01 1.089204e+01 4.326399e-01 +3.000000e-01 1.103999e+01 4.339429e-01 +3.000000e-01 1.118995e+01 4.352366e-01 +3.000000e-01 1.134194e+01 4.365207e-01 +3.000000e-01 1.149600e+01 4.377955e-01 +3.000000e-01 1.165215e+01 4.390606e-01 +3.000000e-01 1.181043e+01 4.403163e-01 +3.000000e-01 1.197085e+01 4.415624e-01 +3.000000e-01 1.213345e+01 4.427990e-01 +3.000000e-01 1.229826e+01 4.440260e-01 +3.000000e-01 1.246531e+01 4.452434e-01 +3.000000e-01 1.263463e+01 4.464513e-01 +3.000000e-01 1.280625e+01 4.476496e-01 +3.000000e-01 1.298020e+01 4.488383e-01 +3.000000e-01 1.315651e+01 4.500174e-01 +3.000000e-01 1.333521e+01 4.511870e-01 +3.000000e-01 1.351635e+01 4.523469e-01 +3.000000e-01 1.369994e+01 4.534972e-01 +3.000000e-01 1.388603e+01 4.546381e-01 +3.000000e-01 1.407465e+01 4.557694e-01 +3.000000e-01 1.426582e+01 4.568911e-01 +3.000000e-01 1.445960e+01 4.580033e-01 +3.000000e-01 1.465601e+01 4.591059e-01 +3.000000e-01 1.485508e+01 4.601990e-01 +3.000000e-01 1.505686e+01 4.612827e-01 +3.000000e-01 1.526138e+01 4.623569e-01 +3.000000e-01 1.546868e+01 4.634216e-01 +3.000000e-01 1.567879e+01 4.644768e-01 +3.000000e-01 1.589176e+01 4.655227e-01 +3.000000e-01 1.610761e+01 4.665592e-01 +3.000000e-01 1.632641e+01 4.675863e-01 +3.000000e-01 1.654817e+01 4.686042e-01 +3.000000e-01 1.677295e+01 4.696127e-01 +3.000000e-01 1.700078e+01 4.706119e-01 +3.000000e-01 1.723170e+01 4.716018e-01 +3.000000e-01 1.746576e+01 4.725826e-01 +3.000000e-01 1.770300e+01 4.735542e-01 +3.000000e-01 1.794346e+01 4.745167e-01 +3.000000e-01 1.818719e+01 4.754700e-01 +3.000000e-01 1.843423e+01 4.764144e-01 +3.000000e-01 1.868462e+01 4.773496e-01 +3.000000e-01 1.893842e+01 4.782760e-01 +3.000000e-01 1.919566e+01 4.791933e-01 +3.000000e-01 1.945640e+01 4.801017e-01 +3.000000e-01 1.972068e+01 4.810013e-01 +3.000000e-01 1.998855e+01 4.818920e-01 +3.000000e-01 2.026006e+01 4.827741e-01 +3.000000e-01 2.053525e+01 4.836473e-01 +3.000000e-01 2.081418e+01 4.845119e-01 +3.000000e-01 2.109690e+01 4.853679e-01 +3.000000e-01 2.138347e+01 4.862152e-01 +3.000000e-01 2.167392e+01 4.870540e-01 +3.000000e-01 2.196832e+01 4.878844e-01 +3.000000e-01 2.226672e+01 4.887063e-01 +3.000000e-01 2.256917e+01 4.895198e-01 +3.000000e-01 2.287573e+01 4.903250e-01 +3.000000e-01 2.318646e+01 4.911219e-01 +3.000000e-01 2.350140e+01 4.919106e-01 +3.000000e-01 2.382062e+01 4.926911e-01 +3.000000e-01 2.414418e+01 4.934635e-01 +3.000000e-01 2.447214e+01 4.942278e-01 +3.000000e-01 2.480454e+01 4.949840e-01 +3.000000e-01 2.514147e+01 4.957324e-01 +3.000000e-01 2.548297e+01 4.964728e-01 +3.000000e-01 2.582911e+01 4.972054e-01 +3.000000e-01 2.617995e+01 4.979301e-01 +3.000000e-01 2.653555e+01 4.986472e-01 +3.000000e-01 2.689599e+01 4.993567e-01 +3.000000e-01 2.726132e+01 5.000584e-01 +3.000000e-01 2.763161e+01 5.007526e-01 +3.000000e-01 2.800694e+01 5.014392e-01 +3.000000e-01 2.838736e+01 5.021185e-01 +3.000000e-01 2.877295e+01 5.027903e-01 +3.000000e-01 2.916378e+01 5.034547e-01 +3.000000e-01 2.955991e+01 5.041120e-01 +3.000000e-01 2.996143e+01 5.047621e-01 +3.000000e-01 3.036840e+01 5.054048e-01 +3.000000e-01 3.078090e+01 5.060406e-01 +3.000000e-01 3.119900e+01 5.066693e-01 +3.000000e-01 3.162278e+01 5.072910e-01 +3.000000e-01 3.205231e+01 5.079058e-01 +3.000000e-01 3.248768e+01 5.085137e-01 +3.000000e-01 3.292897e+01 5.091149e-01 +3.000000e-01 3.337625e+01 5.097092e-01 +3.000000e-01 3.382960e+01 5.102969e-01 +3.000000e-01 3.428911e+01 5.108780e-01 +3.000000e-01 3.475487e+01 5.114524e-01 +3.000000e-01 3.522695e+01 5.120203e-01 +3.000000e-01 3.570544e+01 5.125818e-01 +3.000000e-01 3.619043e+01 5.131370e-01 +3.000000e-01 3.668201e+01 5.136857e-01 +3.000000e-01 3.718027e+01 5.142282e-01 +3.000000e-01 3.768529e+01 5.147644e-01 +3.000000e-01 3.819717e+01 5.152945e-01 +3.000000e-01 3.871601e+01 5.158184e-01 +3.000000e-01 3.924190e+01 5.163364e-01 +3.000000e-01 3.977493e+01 5.168483e-01 +3.000000e-01 4.031519e+01 5.173542e-01 +3.000000e-01 4.086280e+01 5.178544e-01 +3.000000e-01 4.141785e+01 5.183486e-01 +3.000000e-01 4.198043e+01 5.188372e-01 +3.000000e-01 4.255066e+01 5.193200e-01 +3.000000e-01 4.312863e+01 5.197971e-01 +3.000000e-01 4.371445e+01 5.202686e-01 +3.000000e-01 4.430823e+01 5.207346e-01 +3.000000e-01 4.491007e+01 5.211951e-01 +3.000000e-01 4.552009e+01 5.216501e-01 +3.000000e-01 4.613840e+01 5.220998e-01 +3.000000e-01 4.676510e+01 5.225441e-01 +3.000000e-01 4.740032e+01 5.229833e-01 +3.000000e-01 4.804416e+01 5.234171e-01 +3.000000e-01 4.869675e+01 5.238457e-01 +3.000000e-01 4.935821e+01 5.242693e-01 +3.000000e-01 5.002864e+01 5.246878e-01 +3.000000e-01 5.070819e+01 5.251014e-01 +3.000000e-01 5.139697e+01 5.255099e-01 +3.000000e-01 5.209510e+01 5.259135e-01 +3.000000e-01 5.280272e+01 5.263123e-01 +3.000000e-01 5.351994e+01 5.267062e-01 +3.000000e-01 5.424691e+01 5.270954e-01 +3.000000e-01 5.498375e+01 5.274799e-01 +3.000000e-01 5.573060e+01 5.278599e-01 +3.000000e-01 5.648760e+01 5.282351e-01 +3.000000e-01 5.725488e+01 5.286059e-01 +3.000000e-01 5.803258e+01 5.289720e-01 +3.000000e-01 5.882084e+01 5.293338e-01 +3.000000e-01 5.961982e+01 5.296912e-01 +3.000000e-01 6.042964e+01 5.300441e-01 +3.000000e-01 6.125046e+01 5.303928e-01 +3.000000e-01 6.208244e+01 5.307372e-01 +3.000000e-01 6.292571e+01 5.310774e-01 +3.000000e-01 6.378044e+01 5.314134e-01 +3.000000e-01 6.464677e+01 5.317453e-01 +3.000000e-01 6.552488e+01 5.320731e-01 +3.000000e-01 6.641492e+01 5.323969e-01 +3.000000e-01 6.731704e+01 5.327167e-01 +3.000000e-01 6.823141e+01 5.330325e-01 +3.000000e-01 6.915821e+01 5.333444e-01 +3.000000e-01 7.009760e+01 5.336525e-01 +3.000000e-01 7.104974e+01 5.339568e-01 +3.000000e-01 7.201482e+01 5.342573e-01 +3.000000e-01 7.299300e+01 5.345541e-01 +3.000000e-01 7.398448e+01 5.348471e-01 +3.000000e-01 7.498942e+01 5.351366e-01 +3.000000e-01 7.600801e+01 5.354224e-01 +3.000000e-01 7.704044e+01 5.357047e-01 +3.000000e-01 7.808689e+01 5.359835e-01 +3.000000e-01 7.914755e+01 5.362588e-01 +3.000000e-01 8.022263e+01 5.365305e-01 +3.000000e-01 8.131230e+01 5.367990e-01 +3.000000e-01 8.241678e+01 5.370641e-01 +3.000000e-01 8.353625e+01 5.373259e-01 +3.000000e-01 8.467094e+01 5.375843e-01 +3.000000e-01 8.582104e+01 5.378395e-01 +3.000000e-01 8.698676e+01 5.380915e-01 +3.000000e-01 8.816830e+01 5.383404e-01 +3.000000e-01 8.936591e+01 5.385861e-01 +3.000000e-01 9.057978e+01 5.388288e-01 +3.000000e-01 9.181013e+01 5.390682e-01 +3.000000e-01 9.305721e+01 5.393048e-01 +3.000000e-01 9.432121e+01 5.395384e-01 +3.000000e-01 9.560239e+01 5.397689e-01 +3.000000e-01 9.690097e+01 5.399966e-01 +3.000000e-01 9.821719e+01 5.402215e-01 +3.000000e-01 9.955128e+01 5.404434e-01 +3.000000e-01 1.009035e+02 5.406625e-01 +3.000000e-01 1.022741e+02 5.408788e-01 +3.000000e-01 1.036633e+02 5.410924e-01 +3.000000e-01 1.050714e+02 5.413033e-01 +3.000000e-01 1.064986e+02 5.415115e-01 +3.000000e-01 1.079451e+02 5.417171e-01 +3.000000e-01 1.094114e+02 5.419199e-01 +3.000000e-01 1.108975e+02 5.421203e-01 +3.000000e-01 1.124039e+02 5.423181e-01 +3.000000e-01 1.139307e+02 5.425134e-01 +3.000000e-01 1.154782e+02 5.427061e-01 +3.000000e-01 1.170468e+02 5.428964e-01 +3.000000e-01 1.186366e+02 5.430843e-01 +3.000000e-01 1.202481e+02 5.432698e-01 +3.000000e-01 1.218814e+02 5.434529e-01 +3.000000e-01 1.235369e+02 5.436336e-01 +3.000000e-01 1.252150e+02 5.438121e-01 +3.000000e-01 1.269158e+02 5.439882e-01 +3.000000e-01 1.286397e+02 5.441620e-01 +3.000000e-01 1.303870e+02 5.443338e-01 +3.000000e-01 1.321581e+02 5.445032e-01 +3.000000e-01 1.339532e+02 5.446705e-01 +3.000000e-01 1.357727e+02 5.448355e-01 +3.000000e-01 1.376169e+02 5.449985e-01 +3.000000e-01 1.394862e+02 5.451594e-01 +3.000000e-01 1.413809e+02 5.453182e-01 +3.000000e-01 1.433013e+02 5.454749e-01 +3.000000e-01 1.452477e+02 5.456297e-01 +3.000000e-01 1.472207e+02 5.457825e-01 +3.000000e-01 1.492204e+02 5.459332e-01 +3.000000e-01 1.512473e+02 5.460820e-01 +3.000000e-01 1.533017e+02 5.462289e-01 +3.000000e-01 1.553840e+02 5.463740e-01 +3.000000e-01 1.574946e+02 5.465171e-01 +3.000000e-01 1.596339e+02 5.466583e-01 +3.000000e-01 1.618022e+02 5.467978e-01 +3.000000e-01 1.640000e+02 5.469354e-01 +3.000000e-01 1.662276e+02 5.470713e-01 +3.000000e-01 1.684855e+02 5.472053e-01 +3.000000e-01 1.707740e+02 5.473377e-01 +3.000000e-01 1.730937e+02 5.474683e-01 +3.000000e-01 1.754449e+02 5.475972e-01 +3.000000e-01 1.778279e+02 5.477245e-01 +3.000000e-01 1.802434e+02 5.478501e-01 +3.000000e-01 1.826917e+02 5.479741e-01 +3.000000e-01 1.851732e+02 5.480964e-01 +3.000000e-01 1.876884e+02 5.482172e-01 +3.000000e-01 1.902378e+02 5.483363e-01 +3.000000e-01 1.928219e+02 5.484540e-01 +3.000000e-01 1.954410e+02 5.485701e-01 +3.000000e-01 1.980957e+02 5.486847e-01 +3.000000e-01 2.007864e+02 5.487977e-01 +3.000000e-01 2.035137e+02 5.489094e-01 +3.000000e-01 2.062781e+02 5.490196e-01 +3.000000e-01 2.090800e+02 5.491284e-01 +3.000000e-01 2.119200e+02 5.492356e-01 +3.000000e-01 2.147985e+02 5.493416e-01 +3.000000e-01 2.177161e+02 5.494460e-01 +3.000000e-01 2.206734e+02 5.495492e-01 +3.000000e-01 2.236708e+02 5.496510e-01 +3.000000e-01 2.267090e+02 5.497515e-01 +3.000000e-01 2.297884e+02 5.498506e-01 +3.000000e-01 2.329097e+02 5.499485e-01 +3.000000e-01 2.360733e+02 5.500451e-01 +3.000000e-01 2.392799e+02 5.501404e-01 +3.000000e-01 2.425301e+02 5.502344e-01 +3.000000e-01 2.458244e+02 5.503273e-01 +3.000000e-01 2.491635e+02 5.504189e-01 +3.000000e-01 2.525479e+02 5.505093e-01 +3.000000e-01 2.559783e+02 5.505986e-01 +3.000000e-01 2.594553e+02 5.506866e-01 +3.000000e-01 2.629795e+02 5.507734e-01 +3.000000e-01 2.665516e+02 5.508593e-01 +3.000000e-01 2.701722e+02 5.509439e-01 +3.000000e-01 2.738419e+02 5.510274e-01 +3.000000e-01 2.775616e+02 5.511099e-01 +3.000000e-01 2.813318e+02 5.511912e-01 +3.000000e-01 2.851531e+02 5.512714e-01 +3.000000e-01 2.890264e+02 5.513507e-01 +3.000000e-01 2.929523e+02 5.514289e-01 +3.000000e-01 2.969315e+02 5.515060e-01 +3.000000e-01 3.009648e+02 5.515821e-01 +3.000000e-01 3.050528e+02 5.516572e-01 +3.000000e-01 3.091964e+02 5.517314e-01 +3.000000e-01 3.133962e+02 5.518045e-01 +3.000000e-01 3.176531e+02 5.518767e-01 +3.000000e-01 3.219678e+02 5.519480e-01 +3.000000e-01 3.263412e+02 5.520183e-01 +3.000000e-01 3.307739e+02 5.520877e-01 +3.000000e-01 3.352668e+02 5.521561e-01 +3.000000e-01 3.398208e+02 5.522237e-01 +3.000000e-01 3.444367e+02 5.522904e-01 +3.000000e-01 3.491152e+02 5.523562e-01 +3.000000e-01 3.538573e+02 5.524211e-01 +3.000000e-01 3.586638e+02 5.524852e-01 +3.000000e-01 3.635356e+02 5.525484e-01 +3.000000e-01 3.684735e+02 5.526108e-01 +3.000000e-01 3.734785e+02 5.526723e-01 +3.000000e-01 3.785515e+02 5.527331e-01 +3.000000e-01 3.836935e+02 5.527930e-01 +3.000000e-01 3.889052e+02 5.528522e-01 +3.000000e-01 3.941877e+02 5.529105e-01 +3.000000e-01 3.995421e+02 5.529681e-01 +3.000000e-01 4.049691e+02 5.530250e-01 +3.000000e-01 4.104698e+02 5.530811e-01 +3.000000e-01 4.160453e+02 5.531364e-01 +3.000000e-01 4.216965e+02 5.531911e-01 +3.000000e-01 4.274245e+02 5.532450e-01 +3.000000e-01 4.332302e+02 5.532982e-01 +3.000000e-01 4.391148e+02 5.533506e-01 +3.000000e-01 4.450794e+02 5.534024e-01 +3.000000e-01 4.511250e+02 5.534534e-01 +3.000000e-01 4.572527e+02 5.535039e-01 +3.000000e-01 4.634636e+02 5.535536e-01 +3.000000e-01 4.697589e+02 5.536027e-01 +3.000000e-01 4.761397e+02 5.536512e-01 +3.000000e-01 4.826071e+02 5.536990e-01 +3.000000e-01 4.891625e+02 5.537461e-01 +3.000000e-01 4.958068e+02 5.537927e-01 +3.000000e-01 5.025414e+02 5.538386e-01 +3.000000e-01 5.093675e+02 5.538839e-01 +3.000000e-01 5.162863e+02 5.539286e-01 +3.000000e-01 5.232991e+02 5.539728e-01 +3.000000e-01 5.304072e+02 5.540163e-01 +3.000000e-01 5.376118e+02 5.540592e-01 +3.000000e-01 5.449142e+02 5.541016e-01 +3.000000e-01 5.523159e+02 5.541434e-01 +3.000000e-01 5.598180e+02 5.541847e-01 +3.000000e-01 5.674221e+02 5.542254e-01 +3.000000e-01 5.751295e+02 5.542656e-01 +3.000000e-01 5.829415e+02 5.543053e-01 +3.000000e-01 5.908597e+02 5.543444e-01 +3.000000e-01 5.988854e+02 5.543830e-01 +3.000000e-01 6.070202e+02 5.544211e-01 +3.000000e-01 6.152654e+02 5.544586e-01 +3.000000e-01 6.236226e+02 5.544958e-01 +3.000000e-01 6.320934e+02 5.545323e-01 +3.000000e-01 6.406792e+02 5.545684e-01 +3.000000e-01 6.493817e+02 5.546041e-01 +3.000000e-01 6.582023e+02 5.546392e-01 +3.000000e-01 6.671427e+02 5.546739e-01 +3.000000e-01 6.762046e+02 5.547081e-01 +3.000000e-01 6.853896e+02 5.547419e-01 +3.000000e-01 6.946993e+02 5.547752e-01 +3.000000e-01 7.041355e+02 5.548081e-01 +3.000000e-01 7.136999e+02 5.548405e-01 +3.000000e-01 7.233942e+02 5.548725e-01 +3.000000e-01 7.332201e+02 5.549041e-01 +3.000000e-01 7.431796e+02 5.549353e-01 +3.000000e-01 7.532742e+02 5.549660e-01 +3.000000e-01 7.635061e+02 5.549964e-01 +3.000000e-01 7.738769e+02 5.550263e-01 +3.000000e-01 7.843885e+02 5.550559e-01 +3.000000e-01 7.950430e+02 5.550849e-01 +3.000000e-01 8.058422e+02 5.551137e-01 +3.000000e-01 8.167880e+02 5.551420e-01 +3.000000e-01 8.278826e+02 5.551701e-01 +3.000000e-01 8.391278e+02 5.551977e-01 +3.000000e-01 8.505258e+02 5.552250e-01 +3.000000e-01 8.620786e+02 5.552518e-01 +3.000000e-01 8.737883e+02 5.552784e-01 +3.000000e-01 8.856571e+02 5.553045e-01 +3.000000e-01 8.976871e+02 5.553303e-01 +3.000000e-01 9.098806e+02 5.553558e-01 +3.000000e-01 9.222396e+02 5.553810e-01 +3.000000e-01 9.347665e+02 5.554058e-01 +3.000000e-01 9.474635e+02 5.554304e-01 +3.000000e-01 9.603331e+02 5.554545e-01 +3.000000e-01 9.733774e+02 5.554783e-01 +3.000000e-01 9.865989e+02 5.555018e-01 diff --git a/src/enzo/CosmicRayData.h b/src/enzo/CosmicRayData.h new file mode 100644 index 000000000..860c356b4 --- /dev/null +++ b/src/enzo/CosmicRayData.h @@ -0,0 +1,38 @@ +/*********************************************************************** +/ +/ Cosmic Ray Acceleration Efficiency Data +/ +***********************************************************************/ + +struct CosmicRayDataType +{ + // Number of dimension to interpolate over + int CosmicRayGridRank; + + // Maximum number of values for Pre-CR population + int CRMaxNumberPrePopValues; + // Maximum number of values for Mach + int CRMaxNumberMachValues; + + // Cooling grid run file + char *CosmicRayGridRunFile; + + // Values for first parameter space (Pre-Shock CR Percentage) + float *CRPrePopValues; + // Values for second parameter space (Mach Number) + float *CRMachValues; + + //Min/max of the paramters + float CRMinPrePop; + float CRMaxPrePop; + float CRMinMach; + float CRMaxMach; + + // Number of values for first parameter space + int CRNumberPrePopValues; + // Number of values for second parameter space + int CRNumberMachValues; + + // Efficiency Values + float **CREfficiency; +}; diff --git a/src/enzo/Grid_FindShocks.C b/src/enzo/Grid_FindShocks.C new file mode 100644 index 000000000..032c9fa4a --- /dev/null +++ b/src/enzo/Grid_FindShocks.C @@ -0,0 +1,1854 @@ +/*********************************************************************** +/ +/ GRID CLASS (Find Shocks and Accelerate Cosmic Rays) +/ +/ written by: Sam Skillman +/ date: May, 2008 +/ modified1: +/ +/ PURPOSE:Finds all shock mach numbers and injects energy into CR color +/ fields +/ +/ RETURNS: +/ SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include +#include "ErrorExceptions.h" +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "Fluxes.h" +#include "GridList.h" +#include "ExternalBoundary.h" +#include "Grid.h" +#include "fortran.def" +#include "phys_constants.h" +#include "CosmologyParameters.h" + +int GetUnits(float *DensityUnits, float *LengthUnits, + float *TemperatureUnits, float *TimeUnits, + float *VelocityUnits, float *MassUnits, FLOAT Time); + +//Use Unsplit Temperature Jump Identification +int grid::FindShocks() +{ + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + + if (NumberOfBaryonFields == 0) + return SUCCESS; + + this->DebugCheck("FindShocks"); + + long shockcounter, stopcounter; + + float invdx = 1./(CellWidth[0][0]); + float inv2dx = 1./(2.*CellWidth[0][0]); + float inv2dx2 = invdx*invdx; + + int is=GridStartIndex[0]; + int js=GridStartIndex[1]; + int ks=GridStartIndex[2]; + int ie=GridEndIndex[0]; + int je=GridEndIndex[1]; + int ke=GridEndIndex[2]; + + int i, j, k, index, + tempi, posti, prei; + float mu, + preT, postT, tempjumpmag, + gradtx, gradty, gradtz, + maxdiv, temprat, tempmach; + float sarea,vol, + Csound; + float num; + double energyin; + + float prepop,CRfraction; + int thisprepopbin,thismachbin; + + int DensNum, TENum, GENum, + Vel1Num, Vel2Num, Vel3Num; + + int MachNum, CRNum, PSTempNum, PSDenNum; + + /* Compute size (in floats) of the current grid. */ + int size = 1; + for (int dim = 0; dim < GridRank; dim++) + size *= GridDimension[dim]; + + if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, + Vel3Num, TENum) == FAIL) { + ENZO_FAIL("Error in IdentifyPhysicalQuantities."); + } + + // Get CR species fields. + + if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) { + ENZO_FAIL("Error in IdentifyCRSpeciesFields."); + } + + /* Get easy to handle pointers for each variable. */ + + float *density = BaryonField[DensNum]; + float *totalenergy = BaryonField[TENum]; + float *gasenergy = BaryonField[GENum]; + float *velocity1 = BaryonField[Vel1Num]; + float *velocity2 = BaryonField[Vel2Num]; + float *velocity3 = BaryonField[Vel3Num]; + float *mach = BaryonField[MachNum]; + float *cr = BaryonField[CRNum]; + float *pstemp = BaryonField[PSTempNum]; + float *psden = BaryonField[PSDenNum]; + + /* Create temperature, entropy fields */ + float *entropy = new float[size]; + float *tempgrad_dot_entropygrad = new float[size]; + double *flowdivergence = new double[size]; + /* If using cosmology, compute the expansion factor and get units. */ + + float *temperature = new float[size]; + if (this->ComputeTemperatureField(temperature) == FAIL){ + ENZO_FAIL("Error in grid->ComputeTemperatureField."); + } + + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1, MassUnits = 1, aUnits = 1; + + if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, + &TimeUnits, &VelocityUnits, &MassUnits, Time) == FAIL) { + ENZO_FAIL("Error in GetUnits."); + } + + // calculate cell temperature, set default values for mach + float energy; + for (i=0;i 1) + tempgrad_dot_entropygrad[index] += inv2dx2* + ((max(ShockTemperatureFloor,temperature[index+GridDimension[0]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]]))* + (entropy[index+GridDimension[0]]- + entropy[index-GridDimension[0]])); + if (GridRank > 2) + tempgrad_dot_entropygrad[index] += inv2dx2* + ((max(ShockTemperatureFloor,temperature[index+GridDimension[0]*GridDimension[1]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]*GridDimension[1]]))* + (entropy[index+GridDimension[0]*GridDimension[1]]- + entropy[index-GridDimension[0]*GridDimension[1]])); + + flowdivergence[index] = (double)(inv2dx)* + ((double)(velocity1[index+1]) - (double)(velocity1[index-1])); + if (GridRank > 1) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity2[index+GridDimension[0]])- + (double)(velocity2[index-GridDimension[0]])); + if (GridRank > 2) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity3[index+GridDimension[0]*GridDimension[1]])- + (double)(velocity3[index-GridDimension[0]*GridDimension[1]])); + + + + } + } + } + +//Pseudo-Code +/* -----------------Pseudo Code------------------- / + Check a cell if it is is a shock + If not, continue + If yes, then determine temperature gradient + Search for pre/postshock cell + If we come to a higher divergence, break + If we come to a local minimum in divergence, make it the pre/post-shock + cell, and stop looking + If we come to a non-negative divergence, make it the pre/post-shock + cell, and stop looking + Make sure that the temperature/density keep increasing or decreasing + After the ends are found, compute the mach number. + Then compute the CR energy injected, multiplied by the timestep! + + In order to get statistics in terms of pre-shock quantities later, put + the shock in the pre-shock cell(i.e. mach number, cr...) + Done! +/ ----------------------------------------------- */ + if(GridRank < 3){ + ks=0; + ke=0; + } + if(GridRank < 2){ + js=0; + je=0; + } + for(k=ks; k<=ke;k++){ + for(j=js; j<=je;j++){ + for(i=is; i<=ie;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + + if(tempgrad_dot_entropygrad[index] <= 0.0 || + flowdivergence[index] >= 0.0) + continue; + + gradtx = gradty = gradtz = 0.0; + + preT = temperature[index]; + postT = temperature[index]; + + tempjumpmag = + (max(ShockTemperatureFloor,temperature[index+1])-max(ShockTemperatureFloor,temperature[index-1]))* + (max(ShockTemperatureFloor,temperature[index+1])-max(ShockTemperatureFloor,temperature[index-1])); + if (GridRank > 1) + tempjumpmag += + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]]))* + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]])); + if (GridRank > 2) + tempjumpmag += + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]*GridDimension[1]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]*GridDimension[1]]))* + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]*GridDimension[1]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]*GridDimension[1]])); + + tempjumpmag = sqrt(tempjumpmag); + + gradtx = + (max(ShockTemperatureFloor,temperature[index+1])-max(ShockTemperatureFloor,temperature[index-1]))/ + tempjumpmag; + if (GridRank > 1) + gradty = + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]]))/tempjumpmag; + if (GridRank > 2) + gradtz = + (max(ShockTemperatureFloor,temperature[index+GridDimension[0]*GridDimension[1]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]*GridDimension[1]]))/ + tempjumpmag; + + + num=0.0; + maxdiv = flowdivergence[index]; + tempi = index; + while(true){ + //Find next post-cell along temperature gradient + //Make sure you are still in the grid + if( ((i+(int)(num*gradtx)) > (GridDimension[0]-1)) || + ((i+(int)(num*gradtx)) < 0) ) + break; + posti = index + (int)(num*gradtx); + + if (GridRank > 1){ + if( ((j+(int)(num*gradty)) > (GridDimension[1]-1)) || + ((j+(int)(num*gradty)) < 0) ) + break; + posti += ((int)(num*gradty))*GridDimension[0]; + } + if (GridRank > 2){ + if( ((k+(int)(num*gradtz)) > (GridDimension[2]-1)) || + ((k+(int)(num*gradtz)) < 0) ) + break; + posti += ((int)(num*gradtz))*GridDimension[0]*GridDimension[1]; + } + + //If we haven't gone anywhere, increment num. + if(posti == tempi){ + num++; + continue; + } + //Make sure temperature keeps increasing + if(temperature[posti] < postT){ + posti = tempi; + break; + } + //Check for a shock in the current cell. If not, set postT + //and break out. + if(tempgrad_dot_entropygrad[posti] <= 0.0 || + flowdivergence[posti] >= 0.0){ + postT = temperature[posti]; + break; + } + //Check for better center of the shock. If so, get out. + if(flowdivergence[posti] < flowdivergence[index]){ + num=-1; + break; + } + //Check for local maximum in divergence. If so, set + //postT and break out + if(flowdivergence[posti] < maxdiv){ + // postT = temperature[tempi]; //Debatable + postT = temperature[posti]; + break; + } + //Update temporary i, maximum divergence, and increment num. + tempi=posti; + maxdiv = flowdivergence[posti]; + postT = temperature[posti]; + num++; + } + //If a center was found, continue to next cell. + if(num == -1) + continue; + + + //Now find pre-shock cell + num=0.0; + maxdiv = flowdivergence[index]; + tempi = index; + while(true){ + //Find next pre-cell along max(ShockTemperatureFloor,temperature gradient + //Make sure you are still in the grid + if( ((i-(int)(num*gradtx)) > (GridDimension[0]-1)) || + ((i-(int)(num*gradtx)) < 0) ) + break; + prei = index - (int)(num*gradtx); + + if (GridRank > 1){ + if( ((j-(int)(num*gradty)) > (GridDimension[1]-1)) || + ((j-(int)(num*gradty)) < 0) ) + break; + prei -= ((int)(num*gradty))*GridDimension[0]; + } + if (GridRank > 2){ + if( ((k-(int)(num*gradtz)) > (GridDimension[2]-1)) || + ((k-(int)(num*gradtz)) < 0) ) + break; + prei -= ((int)(num*gradtz))*GridDimension[0]*GridDimension[1]; + } + + //If we haven't gone anywhere, increment num. + if(prei == tempi){ + num++; + continue; + } + //Make sure temperature keeps decreasing + if(temperature[prei] > preT){ + prei = tempi; + break; + } + //Check for a shock in the current cell. If not, set preT + //and break out. + if(tempgrad_dot_entropygrad[prei] <= 0.0 || + flowdivergence[prei] >= 0.0){ + preT = temperature[prei]; + break; + } + //Check for better center of the shock. If so, get out. + if(flowdivergence[prei] < flowdivergence[index]){ + num=-1; + break; + } + //Check for local maximum in divergence. If so, set + //preT and break out + if(flowdivergence[prei] < maxdiv){ + // preT = temperature[tempi]; //Debatable + preT = temperature[prei]; + break; + } + //Update temporary i, maximum divergence, and increment num. + tempi=prei; + maxdiv = flowdivergence[prei]; + preT = temperature[prei]; + num++; + } + //If a center was found, continue to next cell. + if(num == -1) + continue; + + temprat = max(ShockTemperatureFloor,postT)/(max(ShockTemperatureFloor,preT)); + //temprat = max(postT,ShockTemperatureFloor)/(max(preT,ShockTemperatureFloor)); + + if(temprat < 1.0) + continue; + + if(density[posti] < density[prei]) + continue; + + tempmach = + sqrt(( 8.0*temprat - 7.0e0 + + sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat) + + 15.0e0) )/5.0e0); + if(tempmach <= 1.0) + continue; + + mach[index] = tempmach; + + + //Calculate mass flux self-consistently + temprat = postT/preT; + tempmach = + sqrt(( 8.0*temprat - 7.0e0 + + sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat) + + 15.0e0) )/5.0e0); + //Speed of sound in code units of velocity + Csound = sqrt(Gamma*kboltz*preT/(DEFAULT_MU*mh))/VelocityUnits; + + /*-------------Lots of Extra Crap------------/ + sarea = CellWidth[0][0] * CellWidth[0][0]; + vol = CellWidth[0][0] * sarea; + //For now calculate the kinetic energy flux through the shock. + + //This is in code units of energy. + //Should multiply by vol*DensityUnits*VelUnits^2 + //to get units of ergs + energyin = 0.5*1.19*sarea*density[prei]* + pow(Csound*mach[index],3.0)*dtFixed; + + //Now in honest to god ergs per cell + energyin *= LengthUnits*LengthUnits*DensityUnits* + VelocityUnits*VelocityUnits*VelocityUnits*TimeUnits; + + //Now convert to the same as gas_energy: + energyin /= vol*LengthUnits*LengthUnits*LengthUnits; + energyin /= density[index]*DensityUnits; + energyin /= VelocityUnits*VelocityUnits; + /--------------------------------------------*/ + //*1.19e0 + //Put energy in postshock cell + energyin = 0.5e0*pow(Csound*tempmach,3.0)* + (density[prei]/density[posti])*dtFixed* + TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + //Put energy in middle cell +// energyin = 0.5e0*pow(Csound*tempmach,3.0)* +// (density[prei]/density[index])*dtFixed* +// TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + prepop = cr[prei]/gasenergy[prei]; + + thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)* + (log10(tempmach) - CosmicRayData.CRMinMach)/ + (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach)); + + thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)* + (prepop - CosmicRayData.CRMinPrePop)/ + (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop)); + + thismachbin = max(0, (int)(thismachbin)); + thismachbin = min((CosmicRayData.CRNumberMachValues -1), + (int)(thismachbin)); + + thisprepopbin = max(0, (int)(thisprepopbin)); + thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1), + (int)(thisprepopbin)); + + CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin]; + + //Want this in Code Units + if(CRModel == 1) + cr[posti]+= CRfraction*energyin; + // cr[index]+= CRfraction*energyin; + + if(CRModel == 2) + cr[index] = CRfraction*energyin/(dtFixed*TimeUnits); + + if(CRModel == 3){ + if(CRfraction*energyin > gasenergy[index] || + CRfraction*energyin > totalenergy[index]){ + gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + cr[index] += 0.9*min(gasenergy[index],totalenergy[index]); + } else { + cr[index] += CRfraction*energyin; + gasenergy[index] -= CRfraction*energyin; + totalenergy[index] -= CRfraction*energyin; + } + } + if(StorePreShockFields){ + pstemp[index] = max(temperature[prei],ShockTemperatureFloor); + psden[index] = density[prei]; + } + } + } + } + + /* deallocate temporary space for solver */ + + delete [] temperature; + delete [] flowdivergence; + delete [] tempgrad_dot_entropygrad; + delete [] entropy; + + return SUCCESS; +} + +//Use Unsplit Velocity Jump Identification +int grid::FindVelShocks() +{ + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + + if (NumberOfBaryonFields == 0) + return SUCCESS; + + this->DebugCheck("FindShocks"); + + long shockcounter, stopcounter; + + float invdx = 1./(CellWidth[0][0]); + float inv2dx = 1./(2.*CellWidth[0][0]); + float inv2dx2 = invdx*invdx; + + int is=GridStartIndex[0]; + int js=GridStartIndex[1]; + int ks=GridStartIndex[2]; + int ie=GridEndIndex[0]; + int je=GridEndIndex[1]; + int ke=GridEndIndex[2]; + + int i, j, k, index, + tempi, posti, prei; + float mu,centervelx,centervely,centervelz, + preV, postV,veljumpmag, + v1jump,v2jump,v3jump, + gradvx, gradvy, gradvz, + maxdiv, thisjump, oldjump, velmach; + float sarea,vol, + Csound; + float num; + double energyin; + + float prepop,CRfraction; + int thisprepopbin,thismachbin; + + int DensNum, TENum, GENum, + Vel1Num, Vel2Num, Vel3Num; + + int MachNum, CRNum, PSTempNum, PSDenNum; + + /* Compute size (in floats) of the current grid. */ + int size = 1; + for (int dim = 0; dim < GridRank; dim++) + size *= GridDimension[dim]; + + if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, + Vel3Num, TENum) == FAIL) { + ENZO_FAIL("Error in IdentifyPhysicalQuantities."); + } + + // Get CR species fields. + + if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) { + ENZO_FAIL("Error in IdentifyCRSpeciesFields."); + } + + /* Get easy to handle pointers for each variable. */ + + float *density = BaryonField[DensNum]; + float *totalenergy = BaryonField[TENum]; + float *gasenergy = BaryonField[GENum]; + float *velocity1 = BaryonField[Vel1Num]; + float *velocity2 = BaryonField[Vel2Num]; + float *velocity3 = BaryonField[Vel3Num]; + float *mach = BaryonField[MachNum]; + float *cr = BaryonField[CRNum]; + float *pstemp = BaryonField[PSTempNum]; + float *psden = BaryonField[PSDenNum]; + + /* Create temperature, entropy fields */ + float *entropy = new float[size]; + float *tempgrad_dot_entropygrad = new float[size]; + double *flowdivergence = new double[size]; + /* If using cosmology, compute the expansion factor and get units. */ + + float *temperature = new float[size]; + + if (this->ComputeTemperatureField(temperature) == FAIL){ + ENZO_FAIL("Error in grid->ComputeTemperatureField."); + } + + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1, MassUnits = 1, aUnits = 1; + + if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, + &TimeUnits, &VelocityUnits, &MassUnits, Time) == FAIL) { + ENZO_FAIL("Error in GetUnits."); + } + + // calculate cell temperature, set default values for mach + float energy; + for (i=0;i 1) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity2[index+GridDimension[0]])- + (double)(velocity2[index-GridDimension[0]])); + if (GridRank > 2) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity3[index+GridDimension[0]*GridDimension[1]])- + (double)(velocity3[index-GridDimension[0]*GridDimension[1]])); + } + } + } + +//Pseudo-Code +/* ----------------------------Pseudo Code---------------------------- / + * Check a cell if it is is a shock If not, continue If yes, then + * determine temperature * gradient Search for pre/postshock cell If we + * come to a higher * divergence, break If we come to a local minimum in + * divergence, make * it the pre/post-shock cell, and stop looking If we + * come to a * non-negative divergence, make it the pre/post-shock cell, + * and stop * looking Make sure that the temperature/density keep + * increasing or * decreasing After the ends are found, compute the mach + * number. Then * compute the CR energy injected, multiplied by the + * timestep! In order to get statistics in terms of pre-shock + * quantities later, put the shock in the pre-shock cell(i.e. mach + * number, cr...) Done! / + * ------------------------------------------------------------------ */ + if(GridRank < 3){ + ks=0; + ke=0; + } + if(GridRank < 2){ + js=0; + je=0; + } + for(k=ks; k<=ke;k++){ + for(j=js; j<=je;j++){ + for(i=is; i<=ie;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + + if(flowdivergence[index] >= 0.0) + continue; + + gradvx = gradvy = gradvz = 0.0; + + preV = index; + postV = index; + + veljumpmag = + (velocity1[index+1]-velocity1[index-1])* + (velocity1[index+1]-velocity1[index-1]); + if (GridRank > 1) + veljumpmag += + (velocity2[index+GridDimension[0]]- + velocity2[index-GridDimension[0]])* + (velocity2[index+GridDimension[0]]- + velocity2[index-GridDimension[0]]); + if (GridRank > 2) + veljumpmag += + (velocity3[index+GridDimension[0]*GridDimension[1]]- + velocity3[index-GridDimension[0]*GridDimension[1]])* + (velocity3[index+GridDimension[0]*GridDimension[1]]- + velocity3[index-GridDimension[0]*GridDimension[1]]); + + veljumpmag = sqrt(veljumpmag); + oldjump = 0.0; + thisjump = 0.0; + + if(veljumpmag == 0.0) + continue; + + gradvx = + (velocity1[index+1]-velocity1[index-1])/veljumpmag; + if (GridRank > 1) + gradvy = + (velocity2[index+GridDimension[0]]- + velocity2[index-GridDimension[0]])/veljumpmag; + if (GridRank > 2) + gradvz = + (velocity3[index+GridDimension[0]*GridDimension[1]]- + velocity3[index-GridDimension[0]*GridDimension[1]])/veljumpmag; + + //Correct gradvx to point in the same direction as temperature gradient + + gradvx = ( (temperature[index+1] > temperature[index-1]) + ? (sqrt(gradvx*gradvx)) : (-sqrt(gradvx*gradvx)) ); + if(GridRank > 1) + gradvy = ( (temperature[index+GridDimension[0]] > + temperature[index-GridDimension[0]]) + ? (sqrt(gradvy*gradvy)) : (-sqrt(gradvy*gradvy)) ); + if(GridRank > 2) + gradvz = ( (temperature[index+GridDimension[0]*GridDimension[1]] > + temperature[index-GridDimension[0]*GridDimension[1]]) + ? (sqrt(gradvz*gradvz)) : (-sqrt(gradvz*gradvz)) ); + + + num=0.0; + centervelx = velocity1[index]; + if(GridRank > 1) + centervely = velocity2[index]; + if(GridRank > 2) + centervelz = velocity3[index]; + + maxdiv = flowdivergence[index]; + tempi = index; + while(true){ + //Find next post-cell along velocity gradient + //Make sure you are still in the grid + if( ((i+(int)(num*gradvx)) > (GridDimension[0]-1)) || + ((i+(int)(num*gradvx)) < 0) ) + break; + posti = index + (int)(num*gradvx); + + if (GridRank > 1){ + if( ((j+(int)(num*gradvy)) > (GridDimension[1]-1)) || + ((j+(int)(num*gradvy)) < 0) ) + break; + posti += ((int)(num*gradvy))*GridDimension[0]; + } + if (GridRank > 2){ + if( ((k+(int)(num*gradvz)) > (GridDimension[2]-1)) || + ((k+(int)(num*gradvz)) < 0) ) + break; + posti += ((int)(num*gradvz))*GridDimension[0]*GridDimension[1]; + } + + //If we haven't gone anywhere, increment num. + if(posti == tempi){ + num++; + continue; + } + v1jump = gradvx*(velocity1[posti] - centervelx); + if(GridRank > 1) + v2jump = gradvy*(velocity2[posti] - centervely); + if(GridRank > 2) + v3jump = gradvz*(velocity3[posti] - centervelz); + + thisjump = v1jump*v1jump; + if(GridRank > 1) + thisjump += v2jump*v2jump; + if(GridRank > 2) + thisjump += v3jump*v3jump; + thisjump = sqrt(thisjump); + + //Make sure velocity (dot) velgrad keeps increasing + if(thisjump < oldjump){ + posti = tempi; + break; + } + //Make sure density/temperature keeps increasing +// if(temperature[posti] < temperature[tempi] || +// density[posti] < density[tempi]){ +// posti = tempi; +// break; +// } + oldjump = thisjump; + //Check for a shock in the current cell. If not, set postV + //and break out. Use last actual shocked cell as cell, not + //the next one. + if(flowdivergence[posti] >= 0.0){ + posti = tempi; + break; + } + //Check for better center of the shock. If so, get out. + if(flowdivergence[posti] < flowdivergence[index]){ + num=-1; + break; + } + //Check for local maximum in divergence. If so, set + //postV and break out + if(flowdivergence[posti] < maxdiv){ + // postV = temperature[tempi]; //Debatable + postV = posti; + break; + } + //Update temporary i, maximum divergence, and increment num. + tempi=posti; + maxdiv = flowdivergence[posti]; + num++; + } + //If a center was found, continue to next cell. + if(num == -1) + continue; + + thisjump = 0.0; + + //Now find pre-shock cell + num=0.0; + maxdiv = flowdivergence[index]; + tempi = index; + while(true){ + //Find next pre-cell along max(ShockTemperatureFloor,temperature gradient + //Make sure you are still in the grid + if( ((i-(int)(num*gradvx)) > (GridDimension[0]-1)) || + ((i-(int)(num*gradvx)) < 0) ) + break; + prei = index - (int)(num*gradvx); + + if (GridRank > 1){ + if( ((j-(int)(num*gradvy)) > (GridDimension[1]-1)) || + ((j-(int)(num*gradvy)) < 0) ) + break; + prei -= ((int)(num*gradvy))*GridDimension[0]; + } + if (GridRank > 2){ + if( ((k-(int)(num*gradvz)) > (GridDimension[2]-1)) || + ((k-(int)(num*gradvz)) < 0) ) + break; + prei -= ((int)(num*gradvz))*GridDimension[0]*GridDimension[1]; + } + + //If we haven't gone anywhere, increment num. + if(prei == tempi){ + num++; + continue; + } + v1jump = gradvx*(velocity1[posti] - velocity1[prei]); + if(GridRank > 1) + v2jump = gradvy*(velocity2[posti] - velocity2[prei]); + if(GridRank > 2) + v3jump = gradvz*(velocity3[posti] - velocity3[prei]); + + thisjump = v1jump*v1jump; + if(GridRank > 1) + thisjump += v2jump*v2jump; + if(GridRank > 2) + thisjump += v3jump*v3jump; + thisjump = sqrt(thisjump); + + //Make sure velocity jump keeps increasing + if(thisjump < oldjump){ + prei = tempi; + break; + } + //Make sure density/temperature keeps decreasing +// if(temperature[prei] > temperature[tempi] || +// density[prei] > density[tempi]){ +// posti = tempi; +// break; +// } + //Check for a shock in the current cell. If not, set preV + //and break out. + if(flowdivergence[prei] >= 0.0){ + preV = prei; + break; + } + //Check for better center of the shock. If so, get out. + if(flowdivergence[prei] < flowdivergence[index]){ + num=-1; + break; + } + //Check for local maximum in divergence. If so, set + //preV and break out + if(flowdivergence[prei] < maxdiv){ + // preV = temperature[tempi]; //Debatable + preV = prei; + break; + } + //Update temporary i, maximum divergence, and increment num. + tempi=prei; + maxdiv = flowdivergence[prei]; + oldjump=thisjump; + num++; + } + //If a center was found, continue to next cell. + if(num == -1) + continue; + +// temprat = max(ShockTemperatureFloor,postV)/(max(ShockTemperatureFloor,preV)); +// //temprat = max(postV,ShockTemperatureFloor)/(max(preV,ShockTemperatureFloor)); + +// if(max(temperature[posti],ShockTemperatureFloor) <= max(temperature[prei],ShockTemperatureFloor)) +// continue; + +// if(density[posti] <= density[prei]) +// continue; + + v1jump = gradvx*(velocity1[posti] - velocity1[prei]); + if(GridRank > 1) + v2jump = gradvy*(velocity2[posti] - velocity2[prei]); + if(GridRank > 2) + v3jump = gradvz*(velocity3[posti] - velocity3[prei]); + + thisjump = v1jump*v1jump; + if(GridRank > 1) + thisjump += v2jump*v2jump; + if(GridRank > 2) + thisjump += v3jump*v3jump; + thisjump = sqrt(thisjump); + + //Speed of sound in code units of velocity + + //Figure out which one is the pre-shock cell: + prei = ( (temperature[prei] < temperature[posti]) ? prei : posti); + + Csound = sqrt(Gamma*kboltz*temperature[prei]/(DEFAULT_MU*mh))/VelocityUnits; + + velmach = (1.0/3.0)*(2.0*thisjump/Csound + + sqrt(9.0+4.0*thisjump*thisjump/Csound/Csound)); + + if(velmach < 1.0) + continue; + + mach[index] = velmach; + + energyin = 0.5e0*pow(thisjump,3.0)* + (density[prei]/density[index])*dtFixed* + TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + //Energy Injection + energyin = 0.5e0*pow(Csound*mach[index],3.0)* + (density[prei]/density[index])*dtFixed* + TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + prepop = cr[prei]/gasenergy[prei]; + + thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)* + (log10(velmach) - CosmicRayData.CRMinMach)/ + (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach)); + + thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)* + (prepop - CosmicRayData.CRMinPrePop)/ + (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop)); + + thismachbin = max(0, (int)(thismachbin)); + thismachbin = min((CosmicRayData.CRNumberMachValues -1), + (int)(thismachbin)); + + thisprepopbin = max(0, (int)(thisprepopbin)); + thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1), + (int)(thisprepopbin)); + + CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin]; + +// fprintf(stdout,"thisprepopbin: %i thismachbin: %i energyin: %e s +// CRfraction: +// %e\n",thisprepopbin,thismachbin,energyin,CRfraction); +// fflush(stdout); + + //Want this in Code Units + if(CRModel == 1) + cr[index]+= CRfraction*energyin; + + if(CRModel == 2) + cr[index] = CRfraction*energyin/(dtFixed*TimeUnits); + + if(CRModel == 3){ + if(CRfraction*energyin > gasenergy[index] || + CRfraction*energyin > totalenergy[index]){ + gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + cr[index] += 0.9*min(gasenergy[index],totalenergy[index]); + } else { + cr[index] += CRfraction*energyin; + gasenergy[index] -= CRfraction*energyin; + totalenergy[index] -= CRfraction*energyin; + } + } + if(StorePreShockFields){ + pstemp[index] = max(temperature[prei],ShockTemperatureFloor); + psden[index] = density[prei]; + } + } + } + } + + /* deallocate temporary space for solver */ + + delete [] temperature; + delete [] flowdivergence; + delete [] tempgrad_dot_entropygrad; + delete [] entropy; + + return SUCCESS; +} + +//Use Split Velocity Jump Identification +int grid::FindVelSplitShocks() +{ + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + + if (NumberOfBaryonFields == 0) + return SUCCESS; + + this->DebugCheck("FindShocks"); + + long shockcounter, stopcounter; + + float invdx = 1./(CellWidth[0][0]); + float inv2dx = 1./(2.*CellWidth[0][0]); + float inv2dx2 = invdx*invdx; + + int is=GridStartIndex[0]; + int js=GridStartIndex[1]; + int ks=GridStartIndex[2]; + int ie=GridEndIndex[0]; + int je=GridEndIndex[1]; + int ke=GridEndIndex[2]; + + int i, j, k, index, + tempi, posti, prei; + float mu,centervelx,centervely,centervelz, + preV, postV,veljumpmag, + v1jump,v2jump,v3jump, + gradvx, gradvy, gradvz, + maxdiv, thisjump, oldjump, velmach; + float sarea,vol, + Csound; + float num; + double energyin; + + float prepop,CRfraction; + int thisprepopbin,thismachbin; + + int DensNum, TENum, GENum, + Vel1Num, Vel2Num, Vel3Num; + + int MachNum, CRNum, PSTempNum, PSDenNum; + + /* Compute size (in floats) of the current grid. */ + int size = 1; + for (int dim = 0; dim < GridRank; dim++) + size *= GridDimension[dim]; + + if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, + Vel3Num, TENum) == FAIL) { + ENZO_FAIL("Error in IdentifyPhysicalQuantities."); + } + + // Get CR species fields. + + if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) { + ENZO_FAIL("Error in IdentifyCRSpeciesFields."); + } + + /* Get easy to handle pointers for each variable. */ + + float *density = BaryonField[DensNum]; + float *totalenergy = BaryonField[TENum]; + float *gasenergy = BaryonField[GENum]; + float *velocity1 = BaryonField[Vel1Num]; + float *velocity2 = BaryonField[Vel2Num]; + float *velocity3 = BaryonField[Vel3Num]; + float *mach = BaryonField[MachNum]; + float *cr = BaryonField[CRNum]; + float *pstemp = BaryonField[PSTempNum]; + float *psden = BaryonField[PSDenNum]; + + /* Create temperature, entropy fields */ + float *entropy = new float[size]; + float *tempgrad_dot_entropygrad = new float[size]; + double *flowdivergence = new double[size]; + /* If using cosmology, compute the expansion factor and get units. */ + + float *temperature = new float[size]; + + if (this->ComputeTemperatureField(temperature) == FAIL){ + ENZO_FAIL("Error in grid->ComputeTemperatureField."); + } + + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1, MassUnits = 1, aUnits = 1; + + if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, + &TimeUnits, &VelocityUnits, &MassUnits, Time) == FAIL) { + ENZO_FAIL("Error in GetUnits."); + } + + // calculate cell temperature, set default values for mach + float energy; + for (i=0;i 1) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity2[index+GridDimension[0]])- + (double)(velocity2[index-GridDimension[0]])); + if (GridRank > 2) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity3[index+GridDimension[0]*GridDimension[1]])- + (double)(velocity3[index-GridDimension[0]*GridDimension[1]])); + } + } + } + +//Pseudo-Code +/* -----------------Pseudo Code------------------- / + Check a cell if it is is a shock + If not, continue + If yes, then determine temperature gradient + Search for pre/postshock cell + If we come to a higher divergence, break + If we come to a local minimum in divergence, make it the pre/post-shock + cell, and stop looking + If we come to a non-negative divergence, make it the pre/post-shock + cell, and stop looking + Make sure that the temperature/density keep increasing or decreasing + After the ends are found, compute the mach number. + Then compute the CR energy injected, multiplied by the timestep! + + In order to get statistics in terms of pre-shock quantities later, put + the shock in the pre-shock cell(i.e. mach number, cr...) + Done! +/ ----------------------------------------------- */ + +//Different starting/ending to make sure we examine all pairs. For 2 or 1D calculations, we only want to analyze k=j=0 + if(GridRank < 3){ + ks=1; + ke=-1; + } + if(GridRank < 2){ + js=1; + je=-1; + } + for(k=ks-1; k<=ke+1;k++){ + for(j=js-1; j<=je+1;j++){ + for(i=is-1; i<=ie+1;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + + if(flowdivergence[index] >= 0.0) + continue; + + //x-direction + + v1jump = fabs(velocity1[index]-velocity1[index-1]); + if(velocity1[index] > velocity1[index-1]) + prei = index-1; + else + prei = index; + + Csound = sqrt(Gamma*kboltz*temperature[prei]/ + (DEFAULT_MU*mh))/VelocityUnits; + mach[prei] = (1.0/3.0)*(2.0*v1jump/Csound + + sqrt(9.0+4.0*v1jump*v1jump/Csound/Csound)); + if(flowdivergence[prei] >= 0.0) + mach[prei]=0.0; + + //y-direction + + if (GridRank > 1){ + v2jump = fabs(velocity2[index]-velocity2[index-GridDimension[0]]); + if(velocity2[index] > velocity2[index-GridDimension[0]]) + prei = index-GridDimension[0]; + else + prei = index; + Csound = sqrt(Gamma*kboltz*temperature[prei]/ + (DEFAULT_MU*mh))/VelocityUnits; + mach[prei] = sqrt(mach[prei]*mach[prei] + + (1.0/3.0)*(2.0*v2jump/Csound + + sqrt(9.0+4.0*v2jump*v2jump/Csound/Csound))* + (1.0/3.0)*(2.0*v2jump/Csound + + sqrt(9.0+4.0*v2jump*v2jump/Csound/Csound))); + if(flowdivergence[prei] >= 0.0) + mach[prei]=0.0; + } + + //z-direction + + if (GridRank > 2){ + v3jump = fabs(velocity3[index]-velocity3[index-GridDimension[0]]); + if(velocity3[index] > velocity3[index- + GridDimension[0]*GridDimension[1]]) + prei = index-GridDimension[0]*GridDimension[1]; + else + prei = index; + Csound = sqrt(Gamma*kboltz*temperature[prei]/ + (DEFAULT_MU*mh))/VelocityUnits; + mach[prei] = sqrt(mach[prei]*mach[prei] + + (1.0/3.0)*(2.0*v3jump/Csound + + sqrt(9.0+4.0*v3jump*v3jump/Csound/Csound))* + (1.0/3.0)*(2.0*v3jump/Csound + + sqrt(9.0+4.0*v3jump*v3jump/Csound/Csound))); + if(flowdivergence[prei] >= 0.0) + mach[prei]=0.0; + } + } + } + } + for(k=ks; k<=ke;k++){ + for(j=js; j<=je;j++){ + for(i=is; i<=ie;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + Csound = sqrt(Gamma*kboltz*temperature[index]/ + (DEFAULT_MU*mh))/VelocityUnits; + //Energy Injection + energyin = 0.5e0*pow(Csound*mach[index],3.0)* + (density[index]/density[index])*dtFixed* + TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + prepop = cr[index]/gasenergy[index]; + + thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)* + (log10(mach[index]) - CosmicRayData.CRMinMach)/ + (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach)); + + thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)* + (prepop - CosmicRayData.CRMinPrePop)/ + (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop)); + + thismachbin = max(0, (int)(thismachbin)); + thismachbin = min((CosmicRayData.CRNumberMachValues -1), + (int)(thismachbin)); + + thisprepopbin = max(0, (int)(thisprepopbin)); + thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1), + (int)(thisprepopbin)); + + CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin]; + + //Want this in Code Units + if(CRModel == 1) + cr[index]+= CRfraction*energyin; + + if(CRModel == 2) + cr[index] = CRfraction*energyin/(dtFixed*TimeUnits); + + if(CRModel == 3){ + if(CRfraction*energyin > gasenergy[index] || + CRfraction*energyin > totalenergy[index]){ + gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + cr[index] += 0.9*min(gasenergy[index],totalenergy[index]); + } else { + cr[index] += CRfraction*energyin; + gasenergy[index] -= CRfraction*energyin; + totalenergy[index] -= CRfraction*energyin; + } + } + if(StorePreShockFields){ + pstemp[index] = max(temperature[prei],ShockTemperatureFloor); + psden[index] = density[prei]; + } + } + } + } + + /* deallocate temporary space for solver */ + + delete [] temperature; + delete [] flowdivergence; + delete [] tempgrad_dot_entropygrad; + delete [] entropy; + + return SUCCESS; +} + +//Use Split Temperature Jump Identification +int grid::FindTempSplitShocks() +{ + /* Return if this doesn't concern us. */ + + if (ProcessorNumber != MyProcessorNumber) + return SUCCESS; + + if (NumberOfBaryonFields == 0) + return SUCCESS; + + this->DebugCheck("FindShocks"); + + long shockcounter, stopcounter; + + float invdx = 1./(CellWidth[0][0]); + float inv2dx = 1./(2.*CellWidth[0][0]); + float inv2dx2 = invdx*invdx; + + int is=GridStartIndex[0]; + int js=GridStartIndex[1]; + int ks=GridStartIndex[2]; + int ie=GridEndIndex[0]; + int je=GridEndIndex[1]; + int ke=GridEndIndex[2]; + + int i, j, k, index, + tempi, posti, prei; + float mu,centervelx,centervely,centervelz, + preV, postV,veljumpmag, + v1jump,v2jump,v3jump, + gradvx, gradvy, gradvz, + maxdiv, thisjump, oldjump, velmach; + float sarea,vol, + Csound; + float num; + double energyin; + + float prepop,CRfraction; + int thisprepopbin,thismachbin; + + //Specific for Split Temperature: + float tempden, temptemp, mach1, + postden, preden, postT, preT, + temprat; + int centerfound; + + int DensNum, TENum, GENum, + Vel1Num, Vel2Num, Vel3Num; + + int MachNum, CRNum, PSTempNum, PSDenNum; + + /* Compute size (in floats) of the current grid. */ + int size = 1; + for (int dim = 0; dim < GridRank; dim++) + size *= GridDimension[dim]; + + if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, + Vel3Num, TENum) == FAIL) { + ENZO_FAIL("Error in IdentifyPhysicalQuantities."); + } + + // Get CR species fields. + + if (IdentifyCRSpeciesFields(MachNum,CRNum,PSTempNum,PSDenNum) == FAIL) { + ENZO_FAIL("Error in IdentifyCRSpeciesFields."); + } + + /* Get easy to handle pointers for each variable. */ + + float *density = BaryonField[DensNum]; + float *totalenergy = BaryonField[TENum]; + float *gasenergy = BaryonField[GENum]; + float *velocity1 = BaryonField[Vel1Num]; + float *velocity2 = BaryonField[Vel2Num]; + float *velocity3 = BaryonField[Vel3Num]; + float *mach = BaryonField[MachNum]; + float *cr = BaryonField[CRNum]; + float *pstemp = BaryonField[PSTempNum]; + float *psden = BaryonField[PSDenNum]; + + /* Create temperature, entropy fields */ + float *entropy = new float[size]; + float *tempgrad_dot_entropygrad = new float[size]; + double *flowdivergence = new double[size]; + /* If using cosmology, compute the expansion factor and get units. */ + + float *temperature = new float[size]; + + if (this->ComputeTemperatureField(temperature) == FAIL){ + ENZO_FAIL("Error in grid->ComputeTemperatureField."); + } + + float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, + VelocityUnits = 1, TimeUnits = 1, MassUnits = 1, aUnits = 1; + + if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, + &TimeUnits, &VelocityUnits, &MassUnits, Time) == FAIL) { + ENZO_FAIL("Error in GetUnits."); + } + + // calculate cell temperature, set default values for mach + float energy; + for (i=0;i 1) + tempgrad_dot_entropygrad[index] += inv2dx2* + ((max(ShockTemperatureFloor,temperature[index+GridDimension[0]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]]))* + (entropy[index+GridDimension[0]]- + entropy[index-GridDimension[0]])); + if (GridRank > 2) + tempgrad_dot_entropygrad[index] += inv2dx2* + ((max(ShockTemperatureFloor,temperature[index+GridDimension[0]*GridDimension[1]])- + max(ShockTemperatureFloor,temperature[index-GridDimension[0]*GridDimension[1]]))* + (entropy[index+GridDimension[0]*GridDimension[1]]- + entropy[index-GridDimension[0]*GridDimension[1]])); + + flowdivergence[index] = (double)(inv2dx)* + ((double)(velocity1[index+1]) - (double)(velocity1[index-1])); + if (GridRank > 1) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity2[index+GridDimension[0]])- + (double)(velocity2[index-GridDimension[0]])); + if (GridRank > 2) + flowdivergence[index] += (double)(inv2dx)* + ((double)(velocity3[index+GridDimension[0]*GridDimension[1]])- + (double)(velocity3[index-GridDimension[0]*GridDimension[1]])); + } + } + } + +//Pseudo-Code +/* -----------------Pseudo Code------------------- / + Check a cell if it is is a shock + If not, continue + If yes, then determine temperature gradient + Search for pre/postshock cell + If we come to a higher divergence, break + If we come to a local minimum in divergence, make it the pre/post-shock + cell, and stop looking + If we come to a non-negative divergence, make it the pre/post-shock + cell, and stop looking + Make sure that the temperature/density keep increasing or decreasing + After the ends are found, compute the mach number. + Then compute the CR energy injected, multiplied by the timestep! + + In order to get statistics in terms of pre-shock quantities later, put + the shock in the pre-shock cell(i.e. mach number, cr...) + Done! +/ ----------------------------------------------- */ + +//Different starting/ending to make sure we examine all pairs. For 2 or 1D calculations, we only want to analyze k=j=0 + if(GridRank < 3){ + ks=1; + ke=-1; + } + if(GridRank < 2){ + js=1; + je=-1; + } + for(k=ks-1; k<=ke+1;k++){ + for(j=js-1; j<=je+1;j++){ + for(i=is-1; i<=ie+1;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + + if(tempgrad_dot_entropygrad[index] <= 0.0 || + flowdivergence[index] >= 0.0) + continue; + + //x-direction + + while(true){ + posti = index; + tempden = density[index]; + maxdiv = flowdivergence[index]; + temptemp = temperature[index]; + if(temperature[index+1] >= temperature[index-1]){ + posti++; + }else{ + posti--; + } + //Make sure temperature/density keeps increasing + if(temperature[posti] < temptemp || + density[posti] < tempden){ + break; + }else{ + temptemp = temperature[posti]; + tempden = density[posti]; + } + + if(flowdivergence[posti] < maxdiv || + flowdivergence[posti] >= 0.0){ + postT = temperature[posti]; + break; + } + if(flowdivergence[posti] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[posti]; + } + while(true){ + prei = index; + tempden = density[index]; + temptemp = temperature[index]; + if(temperature[index+1] < temperature[index-1]){ + prei++; + }else{ + prei--; + } + //Make sure temperature/density keeps decreasing + if(temperature[prei] > temptemp || + density[prei] > tempden){ + break; + }else{ + temptemp = temperature[prei]; + tempden = density[prei]; + } + if(flowdivergence[prei] < maxdiv || + flowdivergence[prei] >= 0.0){ + preT = temperature[prei]; + break; + } + if(flowdivergence[prei] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[prei]; + } + if(centerfound !=1 && + density[posti]>density[prei] && + temperature[posti]>temperature[prei]){ + temprat = max(postT,ShockTemperatureFloor)/(max(preT,ShockTemperatureFloor)); + + mach[index] = + sqrt(( 8.0*temprat - 7.0e0 + + sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat) + + 15.0e0) )/5.0e0); + + } + + + //y-direction + if(GridRank > 1){ + while(true){ + posti = index; + tempden = density[index]; + temptemp = temperature[index]; + maxdiv = flowdivergence[index]; + if(temperature[index+GridDimension[0]] >= + temperature[index-GridDimension[0]]){ + posti+=GridDimension[0]; + }else{ + posti-=GridDimension[0]; + } + //Make sure temperature/density keeps increasing + if(temperature[posti] < temptemp || + density[posti] < tempden){ + break; + }else{ + temptemp = temperature[posti]; + tempden = density[posti]; + } + + if(flowdivergence[posti] < maxdiv || + flowdivergence[posti] >= 0.0){ + postT = temperature[posti]; + break; + } + if(flowdivergence[posti] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[posti]; + } + while(true){ + prei = index; + tempden = density[index]; + temptemp = temperature[index]; + if(temperature[index+GridDimension[0]] < + temperature[index-GridDimension[0]]){ + prei+=GridDimension[0]; + }else{ + prei-=GridDimension[0]; + } + //Make sure temperature/density keeps decreasing + if(temperature[prei] > temptemp || + density[prei] > tempden){ + break; + }else{ + temptemp = temperature[prei]; + tempden = density[prei]; + } + + if(flowdivergence[prei] < maxdiv || + flowdivergence[prei] >= 0.0){ + preT = temperature[prei]; + break; + } + if(flowdivergence[prei] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[prei]; + } + if(centerfound !=1 && + density[posti]>density[prei] && + temperature[posti]>temperature[prei]){ + temprat = max(postT,ShockTemperatureFloor)/(max(preT,ShockTemperatureFloor)); + + mach1 = + sqrt(( 8.0*temprat - 7.0e0 + + sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat) + + 15.0e0) )/5.0e0); + if(mach1 > mach[index]) + mach[index]=mach1; + } + } + + //z-direction + if(GridRank > 2){ + while(true){ + posti = index; + tempden = density[index]; + temptemp = temperature[index]; + maxdiv = flowdivergence[index]; + if(temperature[index+GridDimension[0]*GridDimension[1]] >= + temperature[index-GridDimension[0]*GridDimension[1]]){ + posti+=GridDimension[0]*GridDimension[1]; + }else{ + posti-=GridDimension[0]*GridDimension[1]; + } + //Make sure temperature/density keeps increasing + if(temperature[posti] < temptemp || + density[posti] < tempden){ + break; + }else{ + temptemp = temperature[posti]; + tempden = density[posti]; + } + + if(flowdivergence[posti] < maxdiv || + flowdivergence[posti] >= 0.0){ + postT = temperature[posti]; + break; + } + if(flowdivergence[posti] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[posti]; + } + while(true){ + prei = index; + tempden = density[index]; + temptemp = temperature[index]; + if(temperature[index+GridDimension[0]*GridDimension[1]] < + temperature[index-GridDimension[0]*GridDimension[1]]){ + prei+=GridDimension[0]*GridDimension[1]; + }else{ + prei-=GridDimension[0]*GridDimension[1]; + } + //Make sure temperature/density keeps decreasing + if(temperature[prei] > temptemp || + density[prei] > tempden){ + break; + }else{ + temptemp = temperature[prei]; + tempden = density[prei]; + } + + if(flowdivergence[prei] < maxdiv || + flowdivergence[prei] >= 0.0){ + preT = temperature[prei]; + break; + } + if(flowdivergence[prei] < flowdivergence[index]){ + centerfound = 1; + break; + } + maxdiv = flowdivergence[prei]; + } + if(centerfound !=1 && + density[posti]>density[prei] && + temperature[posti]>temperature[prei]){ + temprat = max(postT,ShockTemperatureFloor)/(max(preT,ShockTemperatureFloor)); + + mach1 = + sqrt(( 8.0*temprat - 7.0e0 + + sqrt( (7.0e0 - 8.0e0*temprat)*(7.0e0 - 8.0e0*temprat) + + 15.0e0) )/5.0e0); + if(mach1 > mach[index]) + mach[index]=mach1; + } + } // GridRank > 2 + } // for i + } // for j + } // for k + for(k=ks; k<=ke;k++){ + for(j=js; j<=je;j++){ + for(i=is; i<=ie;i++){ + + index = i + GridDimension[0]*(j + GridDimension[1]*k); + Csound = sqrt(Gamma*kboltz*temperature[index]/ + (DEFAULT_MU*mh))/VelocityUnits; + //Energy Injection + energyin = 0.5e0*pow(Csound*mach[index],3.0)* + (density[index]/density[index])*dtFixed* + TimeUnits*VelocityUnits/LengthUnits/CellWidth[0][0]; + + prepop = cr[index]/gasenergy[index]; + + thismachbin = (int)((float)(CosmicRayData.CRNumberMachValues)* + (log10(mach[index]) - CosmicRayData.CRMinMach)/ + (CosmicRayData.CRMaxMach - CosmicRayData.CRMinMach)); + + thisprepopbin = (int)((float)(CosmicRayData.CRNumberPrePopValues)* + (prepop - CosmicRayData.CRMinPrePop)/ + (CosmicRayData.CRMaxPrePop - CosmicRayData.CRMinPrePop)); + + thismachbin = max(0, (int)(thismachbin)); + thismachbin = min((CosmicRayData.CRNumberMachValues -1), + (int)(thismachbin)); + + thisprepopbin = max(0, (int)(thisprepopbin)); + thisprepopbin = min((CosmicRayData.CRNumberPrePopValues -1), + (int)(thisprepopbin)); + + CRfraction = CosmicRayData.CREfficiency[thisprepopbin][thismachbin]; + + //Want this in Code Units + if(CRModel == 1) + cr[index]+= CRfraction*energyin; + + if(CRModel == 2) + cr[index] = CRfraction*energyin/(dtFixed*TimeUnits); + + if(CRModel == 3){ + if(CRfraction*energyin > gasenergy[index] || + CRfraction*energyin > totalenergy[index]){ + gasenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + totalenergy[index] -= 0.9*min(gasenergy[index],totalenergy[index]); + cr[index] += 0.9*min(gasenergy[index],totalenergy[index]); + } else { + cr[index] += CRfraction*energyin; + gasenergy[index] -= CRfraction*energyin; + totalenergy[index] -= CRfraction*energyin; + } + } + if(StorePreShockFields){ + pstemp[index] = max(temperature[prei],ShockTemperatureFloor); + psden[index] = density[prei]; + } + } + } + } + + /* deallocate temporary space for solver */ + + delete [] temperature; + delete [] flowdivergence; + delete [] tempgrad_dot_entropygrad; + delete [] entropy; + + return SUCCESS; +} diff --git a/src/enzo/Grid_IdentifyCRSpeciesFields.C b/src/enzo/Grid_IdentifyCRSpeciesFields.C new file mode 100644 index 000000000..9dc0a87b7 --- /dev/null +++ b/src/enzo/Grid_IdentifyCRSpeciesFields.C @@ -0,0 +1,50 @@ +/*********************************************************************** +/ +/ GRID CLASS (IDENTIFY THE SPECIES FIELDS FOR SAM SKILLMAN'S COSMIC RAYS) +/ +/ written by: Samuel Skillman +/ date: May, 2008 +/ modified1: +/ +/ PURPOSE: +/ +/ NOTE: +/ +************************************************************************/ + +#include +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "Fluxes.h" +#include "GridList.h" +#include "ExternalBoundary.h" +#include "Grid.h" +#include "fortran.def" + +/* function prototypes */ + +int FindField(int f, int farray[], int n); + +int grid::IdentifyCRSpeciesFields(int &MachNum,int &CRNum, + int &PSTempNum, int &PSDenNum) +{ + MachNum = CRNum = PSTempNum = PSDenNum = 0; + + // Basic: Mach, CR Protons + if (CRModel) { + MachNum = FindField(Mach, FieldType, NumberOfBaryonFields); + CRNum = FindField(CRDensity, FieldType, NumberOfBaryonFields); + if (StorePreShockFields){ + PSTempNum= FindField(PreShockTemperature, FieldType, NumberOfBaryonFields); + PSDenNum = FindField(PreShockDensity, FieldType, NumberOfBaryonFields); + } + } + + if ((MachNum < 0) || (CRNum < 0) || (PSTempNum < 0) || (PSDenNum < 0)) { + fprintf(stderr,"Error identifying species for CRModel = %"ISYM" MachNum= %"ISYM" CRNum = %"ISYM" PSTempNum = %"ISYM" PSDenNum = %"ISYM" NBaryonFs = %"ISYM".\n", + CRModel,MachNum,CRNum,PSTempNum,PSDenNum,NumberOfBaryonFields); + return FAIL; + } + return SUCCESS; +} diff --git a/src/enzo/Grid_ShocksHandler.C b/src/enzo/Grid_ShocksHandler.C new file mode 100644 index 000000000..6c109f7f1 --- /dev/null +++ b/src/enzo/Grid_ShocksHandler.C @@ -0,0 +1,55 @@ +/*********************************************************************** +/ +/ GRID CLASS (HANDLE CALLING AND SOLVING SHOCK ANALYSIS) +/ +/ written by: Samuel Skillman +/ date: July, 2009 +/ modified1: +/ +/ PURPOSE: Move logic for shock module selection here +/ +/ RETURNS: +/ SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include "ErrorExceptions.h" +#include "performance.h" +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "Fluxes.h" +#include "GridList.h" +#include "ExternalBoundary.h" +#include "Grid.h" + +int grid::ShocksHandler() +{ + if (!CRModel) return SUCCESS; + int shock_status; + JBPERF_START("grid_ShocksHandler"); + + switch(ShockMethod){ + case 0: + shock_status = this->FindShocks(); + break; + case 1: + shock_status = this->FindTempSplitShocks(); + break; + case 2: + shock_status = this->FindVelShocks(); + break; + case 3: + shock_status = this->FindVelSplitShocks(); + break; + default: + shock_status = this->FindShocks(); + } + JBPERF_STOP("grid_ShocksHandler"); + + if(shock_status == FAIL){ + ENZO_FAIL("Error in grid->ShocksHandler."); + } + return SUCCESS; +} diff --git a/src/enzo/InitializeCosmicRayData.C b/src/enzo/InitializeCosmicRayData.C new file mode 100644 index 000000000..d9a1dde02 --- /dev/null +++ b/src/enzo/InitializeCosmicRayData.C @@ -0,0 +1,93 @@ +/*********************************************************************** +/ +/ INITIALIZE THE COSMIC RAY ACCELERATION EFFICIENCIES +/ +/ written by: Sam Skillman +/ date: June, 2008 +/ modified1: +/ +/ PURPOSE: Initializes the cosmic ray efficiency table that is used to +/ interpolate the pre-shock CR component and Mach number +/ +/ RETURNS: SUCCESS or FAIL +/ +************************************************************************/ + +#include +#include +#include +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "CosmologyParameters.h" + +/* function prototypes */ +int InitializeCosmicRayData(void) +{ + /* Open input file for data. */ + + FILE *fptr = fopen("cosmic_ray.dat", "r"); + if (fptr == NULL) { + fprintf(stderr, "Error opening cosmic_ray.dat\n"); + return FAIL; + } + + /* Read CR efficiency data, skipping over comments. */ + //Pre-CR is in normal, Mach in log. + int index = 0; + int i, j; + float precr, mach, creff; + char line[MAX_LINE_LENGTH]; + CosmicRayData.CRNumberPrePopValues = 2; //Pre-CR + CosmicRayData.CRNumberMachValues = 512; //Mach + + CosmicRayData.CREfficiency = + new float*[CosmicRayData.CRNumberPrePopValues]; + for(i=0; i < CosmicRayData.CRNumberPrePopValues; i++) + CosmicRayData.CREfficiency[i] = new float[CosmicRayData.CRNumberMachValues]; + + while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) + if (line[0] != '#') + if (sscanf(line, "%"FSYM" %"FSYM" %"FSYM, &precr, &mach, + &creff) == 3) { + + i = index / CosmicRayData.CRNumberMachValues; + j = index - (i * CosmicRayData.CRNumberMachValues); + + CosmicRayData.CREfficiency[i][j] = creff; + + if (index == 0){ + CosmicRayData.CRMinPrePop = precr; + CosmicRayData.CRMinMach = log10(mach); + } + if (index == (CosmicRayData.CRNumberMachValues -1)) + CosmicRayData.CRMaxMach = log10(mach); + if (index == (CosmicRayData.CRNumberMachValues* + CosmicRayData.CRNumberPrePopValues -1)) + CosmicRayData.CRMaxPrePop = precr; + + index++; + } + fclose(fptr); + + if (index < (CosmicRayData.CRNumberPrePopValues * + CosmicRayData.CRNumberMachValues - 1)){ + printf("Number of entries smaller than expected: %i\n",index); + return FAIL; + } + + if (debug) { + printf("InitializeCosmicRayData: MinPrePop = %"GSYM"\n", + CosmicRayData.CRMinPrePop); + printf("InitializeCosmicRayData: MaxPrePop = %"GSYM"\n", + CosmicRayData.CRMaxPrePop); + printf("InitializeCosmicRayData: MinMach = %"GSYM"\n", + CosmicRayData.CRMinMach); + printf("InitializeCosmicRayData: MaxMach = %"GSYM"\n", + CosmicRayData.CRMaxMach); + printf("InitializeCosmicRayData: Lines Read = %"ISYM"\n", + index); + } + + return SUCCESS; +}