diff --git a/input/metal_cool.dat b/input/metal_cool.dat index dd4be6d22..fc7446103 100644 --- a/input/metal_cool.dat +++ b/input/metal_cool.dat @@ -168,7 +168,7 @@ 6.9499E-29 6.9503E-29 6.9509E-29 6.9516E-29 6.9525E-29 6.9536E-29 6.9550E-29 6.9567E-29 6.9590E-29 6.9618E-29 6.9654E-29 6.9699E-29 6.9756E-29 6.9828E-29 6.9918E-29 7.0033E-29 7.0179E-29 7.0362E-29 7.0594E-29 7.0887E-29 7.1257E-29 7.1725E-29 7.2317E-29 7.3064E-29 7.4009E-29 7.5202E-29 7.6711E-29 7.8618E-29 8.1028E-29 8.4073E-29 8.7922E-29 9.2786E-29 9.8934E-29 1.0670E-28 1.1652E-28 1.2893E-28 1.4461E-28 1.6444E-28 1.8949E-28 2.2115E-28 2.6116E-28 3.1172E-28 3.7563E-28 4.5639E-28 5.5846E-28 6.8745E-28 8.5045E-28 1.0564E-27 1.3167E-27 1.6457E-27 2.0613E-27 2.5864E-27 3.2498E-27 4.0879E-27 5.1465E-27 6.4834E-27 8.1716E-27 1.0303E-26 1.2992E-26 1.6384E-26 7.3342E-29 7.3347E-29 7.3353E-29 7.3360E-29 7.3370E-29 7.3382E-29 7.3397E-29 7.3416E-29 7.3441E-29 7.3471E-29 7.3510E-29 7.3559E-29 7.3621E-29 7.3699E-29 7.3798E-29 7.3923E-29 7.4081E-29 7.4280E-29 7.4532E-29 7.4851E-29 7.5253E-29 7.5762E-29 7.6405E-29 7.7217E-29 7.8244E-29 7.9541E-29 8.1181E-29 8.3253E-29 8.5873E-29 8.9183E-29 9.3366E-29 9.8653E-29 1.0534E-28 1.1378E-28 1.2445E-28 1.3794E-28 1.5499E-28 1.7653E-28 2.0376E-28 2.3817E-28 2.8166E-28 3.3662E-28 4.0608E-28 4.9387E-28 6.0481E-28 7.4501E-28 9.2218E-28 1.1461E-27 1.4290E-27 1.7865E-27 2.2383E-27 2.8090E-27 3.5302E-27 4.4411E-27 5.5918E-27 7.0449E-27 8.8799E-27 1.1196E-26 1.4119E-26 1.7807E-26 7.7359E-29 7.7364E-29 7.7371E-29 7.7379E-29 7.7389E-29 7.7402E-29 7.7419E-29 7.7439E-29 7.7466E-29 7.7499E-29 7.7541E-29 7.7594E-29 7.7661E-29 7.7746E-29 7.7853E-29 7.7989E-29 7.8159E-29 7.8376E-29 7.8649E-29 7.8994E-29 7.9430E-29 7.9981E-29 8.0678E-29 8.1558E-29 8.2671E-29 8.4077E-29 8.5854E-29 8.8100E-29 9.0939E-29 9.4526E-29 9.9060E-29 1.0479E-28 1.1203E-28 1.2118E-28 1.3275E-28 1.4737E-28 1.6584E-28 1.8919E-28 2.1869E-28 2.5599E-28 3.0312E-28 3.6268E-28 4.3796E-28 5.3309E-28 6.5332E-28 8.0526E-28 9.9726E-28 1.2399E-27 1.5465E-27 1.9340E-27 2.4235E-27 3.0421E-27 3.8236E-27 4.8108E-27 6.0578E-27 7.6327E-27 9.6212E-27 1.2132E-26 1.5300E-26 1.9296E-26 - 8.1555E-29 8.1560E-29 8.1567E-29 8.1576E-29 8.1588E-29 8.1602E-29 8.1619E-29 8.1642E-29 8.1670E-29 8.1706E-29 8.1752E-29 8.1809E-29 8.1882E-29 8.1973E-29 8.2089E-29 8.2235E-29 8.2420E-29 8.2654E-29 8.2949E-29 8.3322E-29 8.3793E-29 8.4388E-29 8.5141E-29 8.6092E-29 8.7295E-29 8.8814E-29 9.0735E-29 9.3161E-29 9.6229E-29 1.0010E-28 1.0500E-28 1.1119E-28 1.1902E-28 1.2891E-28 1.4141E-28 1.5720E-28 1.7716E-28 2.0239E-28 2.3428E-28 2.7458E-28 3.2550E-28 3.8987E-28 4.7121E-28 5.7401E-28 7.0393E-28 8.6811E-28 1.0756E-27 1.3378E-27 1.6691E-27 2.0878E-27 2.6168E-27 3.2852E-27 4.1297E-27 5.1965E-27 6.5440E-27 8.2458E-27 1.0395E-26 1.3107E-26 1.6531E-26 2.0850E-26 + 8.1555E-29 8.1560E-29 8.1568E-29 8.1576E-29 8.1588E-29 8.1602E-29 8.1619E-29 8.1642E-29 8.1670E-29 8.1706E-29 8.1752E-29 8.1809E-29 8.1882E-29 8.1973E-29 8.2089E-29 8.2235E-29 8.2420E-29 8.2654E-29 8.2949E-29 8.3322E-29 8.3793E-29 8.4388E-29 8.5141E-29 8.6092E-29 8.7295E-29 8.8814E-29 9.0735E-29 9.3161E-29 9.6229E-29 1.0010E-28 1.0500E-28 1.1119E-28 1.1902E-28 1.2891E-28 1.4141E-28 1.5720E-28 1.7716E-28 2.0239E-28 2.3428E-28 2.7458E-28 3.2550E-28 3.8987E-28 4.7121E-28 5.7401E-28 7.0393E-28 8.6811E-28 1.0756E-27 1.3378E-27 1.6691E-27 2.0878E-27 2.6168E-27 3.2852E-27 4.1297E-27 5.1965E-27 6.5440E-27 8.2458E-27 1.0395E-26 1.3107E-26 1.6531E-26 2.0850E-26 8.5935E-29 8.5941E-29 8.5949E-29 8.5958E-29 8.5970E-29 8.5986E-29 8.6005E-29 8.6029E-29 8.6060E-29 8.6098E-29 8.6147E-29 8.6209E-29 8.6287E-29 8.6386E-29 8.6511E-29 8.6668E-29 8.6867E-29 8.7119E-29 8.7437E-29 8.7839E-29 8.8347E-29 8.8988E-29 8.9800E-29 9.0825E-29 9.2121E-29 9.3757E-29 9.5827E-29 9.8442E-29 1.0175E-28 1.0592E-28 1.1120E-28 1.1787E-28 1.2631E-28 1.3696E-28 1.5043E-28 1.6745E-28 1.8896E-28 2.1615E-28 2.5051E-28 2.9393E-28 3.4881E-28 4.1817E-28 5.0582E-28 6.1659E-28 7.5658E-28 9.3351E-28 1.1571E-27 1.4396E-27 1.7966E-27 2.2478E-27 2.8179E-27 3.5381E-27 4.4481E-27 5.5977E-27 7.0497E-27 8.8836E-27 1.1199E-26 1.4122E-26 1.7811E-26 2.2465E-26 9.0505E-29 9.0511E-29 9.0519E-29 9.0529E-29 9.0542E-29 9.0559E-29 9.0579E-29 9.0605E-29 9.0638E-29 9.0680E-29 9.0732E-29 9.0799E-29 9.0883E-29 9.0989E-29 9.1123E-29 9.1292E-29 9.1506E-29 9.1777E-29 9.2118E-29 9.2550E-29 9.3096E-29 9.3786E-29 9.4657E-29 9.5759E-29 9.7151E-29 9.8910E-29 1.0113E-28 1.0394E-28 1.0750E-28 1.1198E-28 1.1766E-28 1.2483E-28 1.3389E-28 1.4534E-28 1.5981E-28 1.7810E-28 2.0122E-28 2.3044E-28 2.6736E-28 3.1402E-28 3.7300E-28 4.4753E-28 5.4172E-28 6.6077E-28 8.1121E-28 1.0013E-27 1.2416E-27 1.5452E-27 1.9289E-27 2.4137E-27 3.0263E-27 3.8004E-27 4.7783E-27 6.0136E-27 7.5741E-27 9.5448E-27 1.2033E-26 1.5175E-26 1.9139E-26 2.4141E-26 9.5268E-29 9.5275E-29 9.5283E-29 9.5294E-29 9.5308E-29 9.5326E-29 9.5348E-29 9.5376E-29 9.5411E-29 9.5456E-29 9.5512E-29 9.5583E-29 9.5673E-29 9.5787E-29 9.5930E-29 9.6112E-29 9.6341E-29 9.6631E-29 9.6997E-29 9.7460E-29 9.8045E-29 9.8784E-29 9.9719E-29 1.0090E-28 1.0239E-28 1.0428E-28 1.0666E-28 1.0967E-28 1.1348E-28 1.1829E-28 1.2437E-28 1.3205E-28 1.4177E-28 1.5404E-28 1.6955E-28 1.8916E-28 2.1393E-28 2.4524E-28 2.8482E-28 3.3483E-28 3.9804E-28 4.7793E-28 5.7888E-28 7.0647E-28 8.6771E-28 1.0715E-27 1.3290E-27 1.6544E-27 2.0657E-27 2.5853E-27 3.2419E-27 4.0715E-27 5.1196E-27 6.4437E-27 8.1162E-27 1.0228E-26 1.2896E-26 1.6263E-26 2.0512E-26 2.5873E-26 diff --git a/input/metal_cool_ratios.dat b/input/metal_cool_ratios.dat index 91739c750..43f800093 100644 --- a/input/metal_cool_ratios.dat +++ b/input/metal_cool_ratios.dat @@ -15229,7 +15229,7 @@ 269.60 0.00026478 0.26249 0.056946 3.3441E-08 0.14765 0.49034 1.1984E-05 0.023418 0.013208 6.2553E-08 0.0056664 278.02 0.00026417 0.26498 0.055478 3.2602E-08 0.14403 0.49358 1.1526E-05 0.023874 0.012383 5.9447E-08 0.0054052 286.70 0.00026364 0.26727 0.054067 3.1795E-08 0.14054 0.49676 1.1094E-05 0.024328 0.011609 5.6483E-08 0.0051542 - 295.65 0.00026317 0.26936 0.052711 3.1020E-08 0.13719 0.49988 1.0686E-05 0.024780 0.010884 5.3655E-08 0.0049133 + 295.65 0.00026317 0.26936 0.052711 3.1020E-08 0.13719 0.49988 1.0686E-05 0.024780 0.010884 5.3655E-08 0.0049132 304.89 0.00026277 0.27127 0.051408 3.0274E-08 0.13397 0.50297 1.0301E-05 0.025230 0.010205 5.0959E-08 0.0046821 314.41 0.00026245 0.27299 0.050155 2.9557E-08 0.13087 0.50601 9.9368E-06 0.025679 0.0095678 4.8389E-08 0.0044605 324.23 0.00026219 0.27453 0.048951 2.8867E-08 0.12788 0.50902 9.5927E-06 0.026127 0.0089710 4.5941E-08 0.0042483 @@ -22456,7 +22456,7 @@ 427.61 0.00011962 0.67268 0.017967 1.0661E-08 0.047460 0.24324 3.2624E-06 0.013828 0.0027875 1.9419E-08 0.0019153 440.97 0.00012088 0.66985 0.017754 1.0541E-08 0.046954 0.24669 3.1993E-06 0.014150 0.0026424 1.8585E-08 0.0018359 454.74 0.00012221 0.66688 0.017553 1.0429E-08 0.046477 0.25022 3.1401E-06 0.014479 0.0025056 1.7787E-08 0.0017599 - 468.94 0.00012361 0.66375 0.017362 1.0322E-08 0.046027 0.25385 3.0844E-06 0.014816 0.0023765 1.7025E-08 0.0016871 + 468.94 0.00012361 0.66375 0.017363 1.0322E-08 0.046027 0.25385 3.0844E-06 0.014816 0.0023765 1.7025E-08 0.0016871 483.59 0.00012507 0.66048 0.017182 1.0222E-08 0.045603 0.25758 3.0320E-06 0.015160 0.0022545 1.6296E-08 0.0016172 498.69 0.00012660 0.65706 0.017012 1.0127E-08 0.045203 0.26139 2.9827E-06 0.015513 0.0021393 1.5599E-08 0.0015503 514.26 0.00012820 0.65351 0.016850 1.0037E-08 0.044827 0.26529 2.9363E-06 0.015873 0.0020304 1.4932E-08 0.0014861 @@ -24214,7 +24214,7 @@ 107.16 0.00010882 0.64946 0.049681 2.8485E-08 0.12330 0.12351 2.0322E-05 0.0036492 0.037913 1.4578E-07 0.012358 110.51 0.00010536 0.66364 0.046751 2.6830E-08 0.11623 0.12276 1.8486E-05 0.0037098 0.034990 1.3792E-07 0.011797 113.96 0.00010209 0.67686 0.044033 2.5293E-08 0.10966 0.12200 1.6852E-05 0.0037696 0.032301 1.3048E-07 0.011255 - 117.52 9.8990E-05 0.68917 0.041515 2.3869E-08 0.10356 0.12125 1.5397E-05 0.0038291 0.029828 1.2343E-07 0.010733 + 117.52 9.8990E-05 0.68917 0.041515 2.3869E-08 0.10356 0.12125 1.5397E-05 0.0038290 0.029828 1.2343E-07 0.010733 121.19 9.6072E-05 0.70060 0.039184 2.2549E-08 0.097906 0.12053 1.4098E-05 0.0038884 0.027558 1.1676E-07 0.010232 124.98 9.3329E-05 0.71119 0.037028 2.1327E-08 0.092669 0.11983 1.2938E-05 0.0039481 0.025474 1.1047E-07 0.0097521 128.88 9.0760E-05 0.72100 0.035034 2.0196E-08 0.087821 0.11917 1.1900E-05 0.0040084 0.023563 1.0454E-07 0.0092934 @@ -24890,7 +24890,7 @@ 1075.8 0.00011369 0.71539 0.0088105 5.3297E-09 0.024090 0.23408 1.4432E-06 0.016526 0.00048825 4.8598E-09 0.00050420 1109.4 0.00011631 0.71000 0.0088266 5.3428E-09 0.024160 0.23898 1.4454E-06 0.016961 0.00046715 4.6699E-09 0.00048474 1144.0 0.00011900 0.70451 0.0088435 5.3562E-09 0.024233 0.24398 1.4481E-06 0.017405 0.00044693 4.4869E-09 0.00046596 - 1179.7 0.00012176 0.69892 0.0088608 5.3701E-09 0.024307 0.24906 1.4511E-06 0.017859 0.00042757 4.3104E-09 0.00044785 + 1179.7 0.00012176 0.69892 0.0088608 5.3701E-09 0.024307 0.24906 1.4511E-06 0.017859 0.00042757 4.3105E-09 0.00044785 1216.6 0.00012459 0.69323 0.0088787 5.3842E-09 0.024382 0.25422 1.4545E-06 0.018322 0.00040902 4.1404E-09 0.00043039 1254.6 0.00012748 0.68744 0.0088969 5.3985E-09 0.024459 0.25947 1.4582E-06 0.018794 0.00039124 3.9765E-09 0.00041354 1293.8 0.00013044 0.68156 0.0089153 5.4130E-09 0.024536 0.26481 1.4623E-06 0.019275 0.00037421 3.8185E-09 0.00039729 @@ -26089,7 +26089,7 @@ 980.96 7.6485E-05 0.80639 0.0062468 3.7715E-09 0.017021 0.15814 1.0359E-06 0.011069 0.00049967 5.2892E-09 0.00055252 1011.6 7.8440E-05 0.80227 0.0062708 3.7889E-09 0.017110 0.16187 1.0393E-06 0.011390 0.00047957 5.0964E-09 0.00053260 1043.2 8.0456E-05 0.79803 0.0062982 3.8078E-09 0.017203 0.16569 1.0431E-06 0.011720 0.00046030 4.9103E-09 0.00051337 - 1075.8 8.2534E-05 0.79368 0.0063268 3.8274E-09 0.017300 0.16961 1.0472E-06 0.012058 0.00044180 4.7307E-09 0.00049479 + 1075.8 8.2534E-05 0.79369 0.0063268 3.8274E-09 0.017300 0.16961 1.0472E-06 0.012058 0.00044180 4.7307E-09 0.00049479 1109.4 8.4675E-05 0.78923 0.0063563 3.8476E-09 0.017400 0.17362 1.0518E-06 0.012406 0.00042405 4.5574E-09 0.00047684 1144.0 8.6879E-05 0.78468 0.0063867 3.8685E-09 0.017502 0.17772 1.0567E-06 0.012762 0.00040701 4.3900E-09 0.00045951 1179.7 8.9149E-05 0.78001 0.0064180 3.8898E-09 0.017607 0.18191 1.0619E-06 0.013128 0.00039065 4.2285E-09 0.00044276 @@ -29031,7 +29031,7 @@ 141.34 1.7226E-05 0.93573 0.0055510 3.2086E-09 0.013984 0.028297 1.8709E-06 0.0012453 0.0088311 6.3863E-08 0.0063383 145.75 1.6717E-05 0.93815 0.0052483 3.0362E-09 0.013241 0.027893 1.7262E-06 0.0012471 0.0081889 6.0409E-08 0.0060152 150.30 1.6253E-05 0.94035 0.0049724 2.8789E-09 0.012564 0.027530 1.5973E-06 0.0012500 0.0076046 5.7191E-08 0.0057123 - 154.99 1.5832E-05 0.94236 0.0047207 2.7354E-09 0.011946 0.027206 1.4822E-06 0.0012540 0.0070722 5.4190E-08 0.0054284 + 154.99 1.5832E-05 0.94236 0.0047207 2.7354E-09 0.011946 0.027206 1.4822E-06 0.0012540 0.0070722 5.4189E-08 0.0054284 159.84 1.5451E-05 0.94418 0.0044907 2.6043E-09 0.011381 0.026920 1.3791E-06 0.0012590 0.0065863 5.1387E-08 0.0051619 164.83 1.5106E-05 0.94585 0.0042803 2.4843E-09 0.010863 0.026670 1.2867E-06 0.0012651 0.0061424 4.8768E-08 0.0049117 169.97 1.4794E-05 0.94737 0.0040877 2.3743E-09 0.010390 0.026454 1.2036E-06 0.0012723 0.0057360 4.6318E-08 0.0046766 @@ -29679,7 +29679,7 @@ 599.74 1.5098E-05 0.95760 0.0015931 9.5223E-10 0.0042641 0.032648 3.0429E-07 0.0022661 0.00071453 8.5010E-09 0.00089856 618.47 1.5460E-05 0.95685 0.0015967 9.5502E-10 0.0042787 0.033373 3.0356E-07 0.0023274 0.00068759 8.2121E-09 0.00086834 637.79 1.5840E-05 0.95607 0.0016013 9.5834E-10 0.0042958 0.034127 3.0307E-07 0.0023910 0.00066185 7.9343E-09 0.00083926 - 657.70 1.6236E-05 0.95525 0.0016067 9.6219E-10 0.0043152 0.034909 3.0281E-07 0.0024570 0.00063724 7.6670E-09 0.00081126 + 657.70 1.6236E-05 0.95525 0.0016067 9.6219E-10 0.0043152 0.034909 3.0281E-07 0.0024571 0.00063724 7.6670E-09 0.00081126 678.24 1.6651E-05 0.95439 0.0016129 9.6653E-10 0.0043368 0.035721 3.0277E-07 0.0025255 0.00061369 7.4097E-09 0.00078429 699.43 1.7083E-05 0.95349 0.0016199 9.7137E-10 0.0043607 0.036564 3.0295E-07 0.0025966 0.00059116 7.1620E-09 0.00075832 721.27 1.7535E-05 0.95256 0.0016278 9.7670E-10 0.0043868 0.037438 3.0333E-07 0.0026703 0.00056959 6.9235E-09 0.00073328 @@ -31435,7 +31435,7 @@ 141.34 7.0775E-06 0.96363 0.0019639 1.1353E-09 0.0049487 0.015167 7.6498E-07 0.00081142 0.0073771 6.0230E-08 0.0060910 145.75 6.8604E-06 0.96516 0.0018551 1.0734E-09 0.0046818 0.014859 7.0560E-07 0.00080706 0.0068488 5.6977E-08 0.0057774 150.30 6.6632E-06 0.96656 0.0017562 1.0170E-09 0.0044389 0.014577 6.5279E-07 0.00080345 0.0063683 5.3947E-08 0.0054840 - 154.99 6.4845E-06 0.96785 0.0016661 9.6557E-10 0.0042174 0.014321 6.0569E-07 0.00080056 0.0059304 5.1124E-08 0.0052094 + 154.99 6.4845E-06 0.96785 0.0016661 9.6557E-10 0.0042174 0.014321 6.0569E-07 0.00080056 0.0059303 5.1124E-08 0.0052094 159.84 6.3227E-06 0.96902 0.0015839 9.1869E-10 0.0040154 0.014088 5.6359E-07 0.00079836 0.0055307 4.8489E-08 0.0049521 164.83 6.1766E-06 0.97010 0.0015089 8.7587E-10 0.0038308 0.013877 5.2588E-07 0.00079686 0.0051653 4.6028E-08 0.0047109 169.97 6.0451E-06 0.97109 0.0014403 8.3672E-10 0.0036619 0.013687 4.9201E-07 0.00079604 0.0048307 4.3726E-08 0.0044844 @@ -35631,7 +35631,7 @@ 100.77 2.4882E-06 0.95884 2.8558E-06 1.9735E-12 1.0182E-05 0.012437 4.7718E-07 0.00070138 0.016630 1.1540E-07 0.011373 103.92 2.3486E-06 0.96162 2.6597E-06 1.8421E-12 9.4969E-06 0.011861 4.2399E-07 0.00068617 0.015143 1.0785E-07 0.010677 107.16 2.2221E-06 0.96413 2.4838E-06 1.7238E-12 8.8813E-06 0.011330 3.7832E-07 0.00067183 0.013817 1.0092E-07 0.010033 - 110.51 2.1074E-06 0.96642 2.3258E-06 1.6173E-12 8.3272E-06 0.010838 3.3895E-07 0.00065832 0.012633 9.4550E-08 0.0094364 + 110.51 2.1074E-06 0.96642 2.3258E-06 1.6173E-12 8.3272E-06 0.010838 3.3895E-07 0.00065832 0.012633 9.4551E-08 0.0094364 113.96 2.0034E-06 0.96850 2.1836E-06 1.5212E-12 7.8278E-06 0.010384 3.0489E-07 0.00064562 0.011572 8.8688E-08 0.0088843 117.52 1.9090E-06 0.97040 2.0557E-06 1.4344E-12 7.3771E-06 0.0099639 2.7532E-07 0.00063368 0.010622 8.3287E-08 0.0083724 121.19 1.8232E-06 0.97213 1.9403E-06 1.3559E-12 6.9699E-06 0.0095752 2.4957E-07 0.00062248 0.0097676 7.8304E-08 0.0078975 diff --git a/src/enzo/ComputePotentialFieldLevelZero.C b/src/enzo/ComputePotentialFieldLevelZero.C index 652c68a9a..9bc217ae4 100644 --- a/src/enzo/ComputePotentialFieldLevelZero.C +++ b/src/enzo/ComputePotentialFieldLevelZero.C @@ -277,7 +277,7 @@ int ComputePotentialFieldLevelZeroPer(TopGridData *MetaData, float coef = GravitationalConstant/a; // for (int dim = 0; dim < MetaData->TopGridRank; dim++) - // coef *= (DomainRightEdge[dim] - DomainLeftEdge[dim])/float(DomainDim[dim]); + // coef *= (DomainRightEdge[dim] - DomainLeftEdge[dim])/float(DomainDim[dim]); /* Multiply density by Green's function to get potential. */ diff --git a/src/enzo/DepositBaryons.C b/src/enzo/DepositBaryons.C index 93c643cf6..d44fcee01 100644 --- a/src/enzo/DepositBaryons.C +++ b/src/enzo/DepositBaryons.C @@ -36,8 +36,7 @@ int DepositBaryonsChildren(HierarchyEntry *DepositGrid, int DepositBaryons(HierarchyEntry *Grid, FLOAT When) { - - /* Get the time and dt for this grid. Compute time+1/2 dt. */ + /* Get the time and dt for this grid. Compute time+1/2 dt. */ FLOAT TimeMidStep = Grid->GridData->ReturnTime() + When*Grid->GridData->ReturnTimeStep(); @@ -47,7 +46,7 @@ int DepositBaryons(HierarchyEntry *Grid, FLOAT When) /* Set the under_subgrid field (indicating if a cell is refined or not) on this grid. */ - + // printf("DepositBaryons:\n"); HierarchyEntry *Temp = Grid->NextGridNextLevel; Grid->GridData->ZeroSolutionUnderSubgrid(NULL, ZERO_UNDER_SUBGRID_FIELD); while (Temp != NULL) { @@ -78,8 +77,6 @@ int DepositBaryons(HierarchyEntry *Grid, FLOAT When) } -/* this doesn't quite work yet (subgrid zeroing needs to be worked out). */ - int DepositBaryonsChildren(HierarchyEntry *DepositGrid, HierarchyEntry *Grid, FLOAT DepositTime) { diff --git a/src/enzo/EvolveLevel.C b/src/enzo/EvolveLevel.C index bee4e4d18..f3c9c113e 100644 --- a/src/enzo/EvolveLevel.C +++ b/src/enzo/EvolveLevel.C @@ -219,6 +219,8 @@ static int StaticLevelZero = 1; static int StaticLevelZero = 0; #endif +extern int RK2SecondStepBaryonDeposit; + int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], int level, float dtLevelAbove, ExternalBoundary *Exterior) { @@ -253,7 +255,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], LevelHierarchyEntry **SUBlingList; #endif - /* Initialize the chaining mesh used in the FastSiblingLocator. */ if (dbx) fprintf(stderr, "EL: Initialize FSL \n"); @@ -280,8 +281,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], ENZO_FAIL(""); #endif - - /* Clear the boundary fluxes for all Grids (this will be accumulated over the subcycles below (i.e. during one current grid step) and used to by the current grid to correct the zones surrounding this subgrid (step #18). */ @@ -461,7 +460,7 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], //dcc cut second potential cut: Duplicate? - if (ComovingCoordinates && SelfGravity && WritePotential) { + if (SelfGravity && WritePotential) { CopyGravPotential = TRUE; When = 0.0; @@ -471,7 +470,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], PrepareDensityField(LevelArray, level, MetaData, When); #endif // end FAST_SIB - CopyGravPotential = FALSE; for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { if (level <= MaximumGravityRefinementLevel) { @@ -483,6 +481,8 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], Grids[grid1]->GridData->CopyPotentialToBaryonField(); } } // end loop over grids + CopyGravPotential = FALSE; + } // if WritePotential @@ -492,13 +492,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], Grids[grid1]->GridData->DeleteGravitatingMassFieldParticles(); - /* Run the Divergence Cleaing */ - - for (grid1 = 0; grid1 < NumberOfGrids; grid1++) - Grids[grid1]->GridData->PoissonSolver(level); - - - /* ----------------------------------------- */ /* Evolve the next level down (recursively). */ @@ -564,12 +557,7 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], /* Recompute radiation field, if requested. */ - if (RadiationFieldType >= 10 && RadiationFieldType <= 11 && - level <= RadiationFieldLevelRecompute) - if (RadiationFieldUpdate(LevelArray, level, MetaData) == FAIL) { - fprintf(stderr, "Error in RecomputeRadiationField.\n"); - ENZO_FAIL(""); - } + RadiationFieldUpdate(LevelArray, level, MetaData); /* Rebuild the Grids on the next level down. Don't bother on the last cycle, as we'll rebuild this grid soon. */ @@ -589,8 +577,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], LevelZoneCycleCountPerProc[level] += NumberOfCells; } - - cycle++; LevelCycleCount[level]++; @@ -610,12 +596,6 @@ int EvolveLevel(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], /* Clean up. */ -#ifdef UNUSED - if (level > MaximumGravityRefinementLevel && - level == MaximumRefinementLevel) - ZEUSQuadraticArtificialViscosity /= 1; -#endif - delete [] NumberOfSubgrids; delete [] Grids; delete [] SubgridFluxesEstimate; diff --git a/src/enzo/GenerateGridArray.C b/src/enzo/GenerateGridArray.C new file mode 100644 index 000000000..12574e123 --- /dev/null +++ b/src/enzo/GenerateGridArray.C @@ -0,0 +1,49 @@ +/*********************************************************************** +/ +/ EVOLVE LEVEL ROUTINES (CALLED BY EVOLVE LEVEL) +/ +/ written by: Greg Bryan +/ date: June, 1999 +/ +/ ======================================================================= +/ This routine simply converts a linked list of grids into an array of +/ pointers. +/ +************************************************************************/ + +#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 "Hierarchy.h" +#include "LevelHierarchy.h" + +int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, + HierarchyEntry **Grids[]) +{ + /* Count the number of grids on this level. */ + + int NumberOfGrids = 0, counter = 0; + LevelHierarchyEntry *Temp = LevelArray[level]; + while (Temp != NULL) { + NumberOfGrids++; + Temp = Temp->NextGridThisLevel; + } + /* Create a list of pointers and number of subgrids (and fill it out). */ + + typedef HierarchyEntry* HierarchyEntryPointer; + *Grids = new HierarchyEntryPointer[NumberOfGrids]; + Temp = LevelArray[level]; + while (Temp != NULL) { + (*Grids)[counter++] = Temp->GridHierarchyEntry; + Temp = Temp->NextGridThisLevel; + } + + return NumberOfGrids; +} + + diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index 67e481722..1cc69184d 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -1955,8 +1955,8 @@ int CollapseTestInitializeGrid(int NumberOfSpheres, float fluxcoef, int fallback); int TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FLOAT CloudRadius, float CloudMachNumber, float CloudAngularVelocity, float InitialBField, - int SetTurbulence, int CloudType, int TurbulenceSeed, int level, - int SetBaryonFields); + int SetTurbulence, int CloudType, int TurbulenceSeed, int PutSink, + int level, int SetBaryonFields); int Collapse3DInitializeGrid(int n_sphere, FLOAT r_sphere[MAX_SPHERES], FLOAT rc_sphere[MAX_SPHERES], diff --git a/src/enzo/Grid_AccelerationBoundaryRoutines.C b/src/enzo/Grid_AccelerationBoundaryRoutines.C index 2bd337e27..434a8fa17 100644 --- a/src/enzo/Grid_AccelerationBoundaryRoutines.C +++ b/src/enzo/Grid_AccelerationBoundaryRoutines.C @@ -122,29 +122,28 @@ int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids, int CycleNumber) { - if ( ! (SelfGravity || UniformGravity || PointSourceGravity) && level > 0 ) + if ( ! ( (SelfGravity || UniformGravity || PointSourceGravity) && level > 0 )) return SUCCESS; //Set the boundary on the Acceleration field. Reuse SetBoundaryConditions. //Juggle pointers around. - int grid, ConservativeTruth; + int grid1, ConservativeTruth; char basename[30]; //We don't want conservative interpolation actually being done for the acceleration field. ConservativeTruth = ConservativeInterpolation; ConservativeInterpolation = FALSE; - for (grid = 0; grid < NumberOfGrids; grid++) { - if( Grids[grid]->GridData->AttachAcceleration() == FAIL ) { + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { + if( Grids[grid1]->GridData->AttachAcceleration() == FAIL ) { fprintf(stderr,"Error in AttachAcceleration \n"); ENZO_FAIL(""); } - if( Grids[grid]->ParentGrid->GridData->AttachAcceleration() ==FAIL ){ + if( Grids[grid1]->ParentGrid->GridData->AttachAcceleration() ==FAIL ){ fprintf(stderr,"Error in AttachAcceleration, Parent \n"); ENZO_FAIL(""); } - } #ifdef FAST_SIB @@ -159,13 +158,13 @@ int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids, - for (grid = 0; grid < NumberOfGrids; grid++) { + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { - if( Grids[grid]->GridData->DetachAcceleration() == FAIL ) { + if( Grids[grid1]->GridData->DetachAcceleration() == FAIL ) { fprintf(stderr,"Error in DetachAcceleration\n"); ENZO_FAIL(""); } - if( Grids[grid]->ParentGrid->GridData->DetachAcceleration() == FAIL ) { + if( Grids[grid1]->ParentGrid->GridData->DetachAcceleration() == FAIL ) { fprintf(stderr,"Error in DetachAcceleration, parent\n"); ENZO_FAIL(""); } diff --git a/src/enzo/Grid_ComputeTemperatureField.C b/src/enzo/Grid_ComputeTemperatureField.C index e7de3a4ef..e3658c9e0 100644 --- a/src/enzo/Grid_ComputeTemperatureField.C +++ b/src/enzo/Grid_ComputeTemperatureField.C @@ -30,7 +30,7 @@ /* Set the mean molecular mass. */ -#define DEFAULT_MU 0.6 +//#define DEFAULT_MU 0.6 #define MU_METAL 16.0 /* This is minimum returned temperature. (K) */ @@ -122,7 +122,9 @@ int grid::ComputeTemperatureField(float *temperature) /* For Sedov Explosion compute temperature without floor */ - float mol_weight = DEFAULT_MU, min_temperature = 1.0; + float mol_weight = Mu, min_temperature = 1.0; // Mu is now read from parameter file + + // this : would not be needed as we can specify it in parameter file: if (ProblemType == 7) {//AK for Sedov explosion test mol_weight = 1.0; min_temperature = tiny_number; @@ -131,7 +133,7 @@ int grid::ComputeTemperatureField(float *temperature) if (MultiSpecies == FALSE) /* If the multi-species flag is not set, - Compute temperature T = p/d and assume mu = DEFAULT_MU. */ + Compute temperature T = p/d and assume mu = Mu (global data). */ for (i = 0; i < size; i++) temperature[i] = max((TemperatureUnits*temperature[i]*mol_weight diff --git a/src/enzo/Grid_ComputeTimeStep.C b/src/enzo/Grid_ComputeTimeStep.C index 1d6f018c9..1a1dcc817 100644 --- a/src/enzo/Grid_ComputeTimeStep.C +++ b/src/enzo/Grid_ComputeTimeStep.C @@ -420,7 +420,7 @@ float grid::ComputeTimeStep() if (HydroMethod != MHD_RK && NumberOfBaryonFields > 0) printf("Bar = %"FSYM" ", dtBaryons); if (HydroMethod == MHD_RK) - printf("dtMHD = %"FSYM" ", dtMHD); + printf("dtMHD = %e ", dtMHD); if (HydroMethod == Zeus_Hydro) printf("Vis = %"FSYM" ", dtViscous); if (ComovingCoordinates) diff --git a/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C b/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C index 4d167a9c1..5d71ac161 100644 --- a/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C +++ b/src/enzo/Grid_CopyBaryonFieldToOldBaryonField.C @@ -48,11 +48,6 @@ int grid::CopyBaryonFieldToOldBaryonField() int size = 1; -/* BUG reported 4th June 2006 - for (int dim = 0; dim < GridDimension[dim]; dim++) - size *= GridDimension[dim]; -*/ - for (int dim = 0; dim < GridRank; dim++) { size *= GridDimension[dim]; } diff --git a/src/enzo/Grid_CopyPotentialToBaryonField.C b/src/enzo/Grid_CopyPotentialToBaryonField.C index c63c6c769..57e2ea132 100644 --- a/src/enzo/Grid_CopyPotentialToBaryonField.C +++ b/src/enzo/Grid_CopyPotentialToBaryonField.C @@ -101,7 +101,7 @@ int grid::CopyPotentialToBaryonField() for (i = 0; i < GridDimension[0]; i++, index++) { - // BaryonField[field+1][jj] = GravitatingMassField[index]; // use this for debugging + BaryonField[field+1][jj] = GravitatingMassField[index]; // use this for debugging BaryonField[field][jj++] = PotentialField[index]; // debuggin: maxPot = max(maxPot,PotentialField[index]); diff --git a/src/enzo/Grid_DepositBaryons.C b/src/enzo/Grid_DepositBaryons.C index 059325a04..07afc1c0c 100644 --- a/src/enzo/Grid_DepositBaryons.C +++ b/src/enzo/Grid_DepositBaryons.C @@ -51,6 +51,8 @@ extern "C" void FORTRAN_NAME(dep_grid_cic)( int *ddim1, int *ddim2, int *ddim3, int *refine1, int *refine2, int *refine3); +int RK2SecondStepBaryonDeposit = 0; + /* InterpolateBoundaryFromParent function */ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime) @@ -145,6 +147,8 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime) if (RegionDim[dim] < 2) { fprintf(stderr, "RegionDim[%"ISYM"] = %"ISYM" < 2\n", dim, RegionDim[dim]); + fprintf(stderr, "GridOffsetEnd[%"ISYM"] = %"ISYM" < 2\n", dim, GridOffsetEnd[dim]); + fprintf(stderr, "GridOffset[%"ISYM"] = %"ISYM" < 2\n", dim, GridOffset[dim]); ENZO_FAIL(""); } @@ -172,8 +176,9 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime) ENZO_FAIL(""); } float dt = (DepositTime - Time)/a; + // if (HydroMethod < 3) dt = 0; dt = 0; - + /* Set up a float version of cell size to pass to fortran. */ float dxfloat[MAX_DIMENSION] = {0,0,0}; @@ -188,21 +193,50 @@ int grid::DepositBaryons(grid *TargetGrid, FLOAT DepositTime) velocity field. */ // fprintf(stderr, "Grid_DepositBaryons - call dep_grid_cic\n"); - - FORTRAN_NAME(dep_grid_cic)(BaryonField[DensNum], dens_field, vel_field, - BaryonField[Vel1Num], BaryonField[Vel2Num], - BaryonField[Vel3Num], &dt, - BaryonField[NumberOfBaryonFields], &GridRank, - &HydroMethod, - dxfloat, dxfloat+1, dxfloat+2, - GridDimension, GridDimension+1, GridDimension+2, - GridStartIndex, GridStartIndex+1, GridStartIndex+2, - GridEndIndex, GridEndIndex+1, GridEndIndex+2, - GridStart, GridStart+1, GridStart+2, - RegionDim, RegionDim+1, RegionDim+2, - Refinement, Refinement+1, Refinement+2); + + float *input_density = BaryonField[DensNum]; + float *input_velx = BaryonField[Vel1Num]; + float *input_vely = BaryonField[Vel2Num]; + float *input_velz = BaryonField[Vel3Num]; + + float *av_dens = NULL; + // supply zero velocity field and the time averaged density to the deposit routine for second step + if (RK2SecondStepBaryonDeposit == 1) { + int current_size = 1; + dt = 0; + for (dim = 0; dim < GridRank; dim++) + current_size *= GridDimension[dim]; + + if (OldBaryonField[DensNum] != NULL) { + av_dens = new float[current_size]; + for (int i=0; i #endif /* USE_MPI */ @@ -75,182 +80,7 @@ int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, extern int CopyPotentialFieldAverage; #define GRIDS_PER_LOOP 20000 - -/* ======================================================================= */ -/* This routine sets all the boundary conditions for Grids by either - interpolating from their parents or copying from sibling grids. */ - - - -#ifdef FAST_SIB -int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, - SiblingGridList SiblingList[], - int level, TopGridData *MetaData, - ExternalBoundary *Exterior, LevelHierarchyEntry *Level) -#else -int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, - int level, TopGridData *MetaData, - ExternalBoundary *Exterior, LevelHierarchyEntry *Level) -#endif -{ - - int loopEnd = (ShearingBoundaryDirection != 1) ? 2 : 1; - - int grid1, grid2, StartGrid, EndGrid, loop; - - JBPERF_START("SetBoundaryConditions"); - - for (loop = 0; loop < loopEnd; loop++){ - -#ifdef FORCE_MSG_PROGRESS - CommunicationBarrier(); -#endif - - TIME_MSG("Interpolating boundaries from parent"); - - for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { - - if (traceMPI) fprintf(tracePtr, "SBC loop\n"); - - EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); - - /* -------------- FIRST PASS ----------------- */ - /* Here, we just generate the calls to generate the receive buffers, - without actually doing anything. */ - - CommunicationDirection = COMMUNICATION_POST_RECEIVE; - CommunicationReceiveIndex = 0; - - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* a) Interpolate boundaries from the parent grid or set external - boundary conditions. */ - - CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - - if (level == 0) { - if (loop == 0) - Grids[grid1]->GridData->SetExternalBoundaryValues(Exterior); - } else { - Grids[grid1]->GridData->InterpolateBoundaryFromParent - (Grids[grid1]->ParentGrid->GridData); - } - - } // ENDFOR grids - - /* -------------- SECOND PASS ----------------- */ - /* Now we generate all the sends, and do all the computation - for grids which are on the same processor as well. */ - - CommunicationDirection = COMMUNICATION_SEND; - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* a) Interpolate boundaries from the parent grid or set - external boundary conditions. */ - - if (level > 0) - Grids[grid1]->GridData->InterpolateBoundaryFromParent - (Grids[grid1]->ParentGrid->GridData); - - - } // ENDFOR grids - //Grids[StartGrid]->GridData->PrintToScreenBoundaries(0); - - /* -------------- THIRD PASS ----------------- */ - /* In this final step, we get the messages as they come in and - then match them to the methods which generate the receive - handle. */ - - if (CommunicationReceiveHandler() == FAIL) - ENZO_FAIL(""); - - } // ENDFOR grid batches - - TIME_MSG("Copying zones in SetBoundaryConditions"); - for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { - EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); - - /* -------------- FIRST PASS ----------------- */ - /* b) Copy any overlapping zones for sibling grids. */ - - CommunicationDirection = COMMUNICATION_POST_RECEIVE; - CommunicationReceiveIndex = 0; - -#ifdef FAST_SIB - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) - for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++) - Grids[grid1]->GridData-> - CheckForOverlap(SiblingList[grid1].GridList[grid2], - MetaData->LeftFaceBoundaryCondition, - MetaData->RightFaceBoundaryCondition, - &grid::CopyZonesFromGrid); -#else - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) - for (grid2 = 0; grid2 < NumberOfGrids; grid2++) - Grids[grid1]->GridData-> - CheckForOverlap(Grids[grid2]->GridData, - MetaData->LeftFaceBoundaryCondition, - MetaData->RightFaceBoundaryCondition, - &grid::CopyZonesFromGrid); -#endif - - /* -------------- SECOND PASS ----------------- */ - /* b) Copy any overlapping zones for sibling grids. */ - - CommunicationDirection = COMMUNICATION_SEND; -#ifdef FAST_SIB - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) - for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++) - Grids[grid1]->GridData-> - CheckForOverlap(SiblingList[grid1].GridList[grid2], - MetaData->LeftFaceBoundaryCondition, - MetaData->RightFaceBoundaryCondition, - &grid::CopyZonesFromGrid); -#else - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) - for (grid2 = 0; grid2 < NumberOfGrids; grid2++) - Grids[grid1]->GridData-> - CheckForOverlap(Grids[grid2]->GridData, - MetaData->LeftFaceBoundaryCondition, - MetaData->RightFaceBoundaryCondition, - &grid::CopyZonesFromGrid); -#endif - - /* -------------- THIRD PASS ----------------- */ - - if (CommunicationReceiveHandler() == FAIL) - ENZO_FAIL(""); - - } // end loop over batchs of grids - - } // ENDFOR loop (for ShearinBox) - - /* c) Apply external reflecting boundary conditions, if needed. */ - - for (grid1 = 0; grid1 < NumberOfGrids; grid1++) - Grids[grid1]->GridData->CheckForExternalReflections - (MetaData->LeftFaceBoundaryCondition, - MetaData->RightFaceBoundaryCondition); - -#ifdef FORCE_MSG_PROGRESS - CommunicationBarrier(); -#endif - - CommunicationDirection = COMMUNICATION_SEND_RECEIVE; - - JBPERF_STOP("SetBoundaryConditions"); - - return SUCCESS; - -} - - - - /* ======================================================================= */ - /* This routine prepares the density field for all the grids on this level, - both particle and baryonic densities. It also calculates the potential - field if this is level 0 (since this involves communication). */ #ifdef FAST_SIB int PrepareDensityField(LevelHierarchyEntry *LevelArray[], @@ -712,316 +542,3 @@ int PrepareDensityField(LevelHierarchyEntry *LevelArray[], } - - /* ======================================================================= */ - /* This routines does the flux correction and project for all grids on this - level from the list of subgrids. */ - -#ifdef FLUX_FIX -int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids, - int NumberOfSubgrids[], - fluxes **SubgridFluxesEstimate[], - LevelHierarchyEntry* SUBlingList[], - TopGridData *MetaData) -#else -int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids, - int NumberOfSubgrids[], - fluxes **SubgridFluxesEstimate[]) -#endif - -{ - - int grid1, subgrid, StartGrid, EndGrid; - HierarchyEntry *NextGrid; - -#ifdef FLUX_FIX - int SUBlingGrid; - LevelHierarchyEntry *NextEntry; -#endif - - /* Define a temporary flux holder for the refined fluxes. */ - - fluxes SubgridFluxesRefined; - - /* For each grid, - (a) project the subgrid's solution into this grid (step #18) - (a1) project the neighboring subgrid's solution into this grid - (b) correct for the difference between this grid's fluxes and the - subgrid's fluxes. (step #19) */ - -#ifdef FORCE_MSG_PROGRESS - CommunicationBarrier(); -#endif - - TIME_MSG("UpdateFromFinerGrids"); - for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { - EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); - - /* -------------- FIRST PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_POST_RECEIVE; - CommunicationReceiveIndex = 0; - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* Loop over subgrids for this grid. */ - - NextGrid = Grids[grid1]->NextGridNextLevel; - subgrid = 0; - CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - - while (NextGrid != NULL && FluxCorrection) { - - /* Project subgrid's refined fluxes to the level of this grid. */ - -#ifdef USE_MPI - CommunicationReceiveArgumentInt[0][CommunicationReceiveIndex] = grid1; - CommunicationReceiveArgumentInt[1][CommunicationReceiveIndex] = subgrid; -#endif /* USE_MPI */ - - NextGrid->GridData-> - GetProjectedBoundaryFluxes(Grids[grid1]->GridData, SubgridFluxesRefined); - - NextGrid = NextGrid->NextGridThisLevel; - subgrid++; - } // ENDWHILE subgrids - - } // ENDFOR grids - - /* -------------- SECOND PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_SEND; - - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* Loop over subgrids for this grid. */ - - NextGrid = Grids[grid1]->NextGridNextLevel; - subgrid = 0; - while (NextGrid != NULL && FluxCorrection) { - - /* Project subgrid's refined fluxes to the level of this grid. */ - - NextGrid->GridData-> - GetProjectedBoundaryFluxes(Grids[grid1]->GridData, SubgridFluxesRefined); - - /* Correct this grid for the refined fluxes (step #19) - (this also deletes the fields in SubgridFluxesRefined). - (only call it if the grid and sub-grid are on the same - processor, otherwise handled in CommunicationReceiveHandler.) */ - -#ifdef FLUX_FIX - if (NextGrid->GridData->ReturnProcessorNumber() == - Grids[grid1]->GridData->ReturnProcessorNumber()) - Grids[grid1]->GridData->CorrectForRefinedFluxes - (SubgridFluxesEstimate[grid1][subgrid], &SubgridFluxesRefined, - SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], - FALSE, MetaData); -#else - if (NextGrid->GridData->ReturnProcessorNumber() == - Grids[grid1]->GridData->ReturnProcessorNumber()) - Grids[grid1]->GridData->CorrectForRefinedFluxes - (SubgridFluxesEstimate[grid1][subgrid], &SubgridFluxesRefined, - SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1]); -#endif /* FLUX_FIX */ - - NextGrid = NextGrid->NextGridThisLevel; - subgrid++; - - } // ENDWHILE subgrids - - } // ENDFOR grids - - /* -------------- THIRD PASS ----------------- */ - -#ifdef FLUX_FIX - CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids, - FALSE, MetaData); -#else - CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids); -#endif - - } // ENDFOR grid batches - - /************************************************************************ - (a1) project the neighboring subgrid's solution into this grid - (SUBlings) - ************************************************************************/ - -#ifdef FLUX_FIX - TIME_MSG("Projecting neighboring subgrid solution (FLUX_FIX)"); - for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { - EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); - - /* -------------- FIRST PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_POST_RECEIVE; - CommunicationReceiveIndex = 0; - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* Loop over subgrids for this grid. */ - - CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - NextEntry = SUBlingList[grid1]; - - while (NextEntry != NULL && FluxCorrection) { - - /* make sure this isn't a "proper" subgrid */ - - if (NextEntry->GridHierarchyEntry->ParentGrid != Grids[grid1]) { - - /* Project subgrid's refined fluxes to the level of this grid. */ - -#ifdef USE_MPI - CommunicationReceiveArgumentInt[0][CommunicationReceiveIndex] = grid1; - CommunicationReceiveArgumentInt[1][CommunicationReceiveIndex] = - NumberOfSubgrids[grid1]-1; -#endif /* USE_MPI */ - - NextEntry->GridData->GetProjectedBoundaryFluxes - (Grids[grid1]->GridData, SubgridFluxesRefined); - } - - NextEntry = NextEntry->NextGridThisLevel; - - } // ENDWHILE subgrids - } // ENDFOR grids - - /* -------------- SECOND PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_SEND; - - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - NextEntry = SUBlingList[grid1]; - - while (NextEntry != NULL && FluxCorrection) { - - /* make sure this isn't a "proper" subgrid */ - - if (NextEntry->GridHierarchyEntry->ParentGrid != Grids[grid1]) { - - /* Project subgrid's refined fluxes to the level of this grid. */ - - NextEntry->GridData->GetProjectedBoundaryFluxes - (Grids[grid1]->GridData, SubgridFluxesRefined); - - /* Correct this grid for the refined fluxes (step #19) - (this also deletes the fields in SubgridFluxesRefined). */ - - if (NextEntry->GridData->ReturnProcessorNumber() == - Grids[grid1]->GridData->ReturnProcessorNumber()) - Grids[grid1]->GridData->CorrectForRefinedFluxes - (SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], - &SubgridFluxesRefined, - SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], - TRUE, MetaData); - } - - NextEntry = NextEntry->NextGridThisLevel; - - } // ENDWHILE subgrids - } // ENDFOR grids - - /* -------------- THIRD PASS ----------------- */ - - CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids, - TRUE, MetaData); - - } // ENDFOR grid batches -#endif /* FLUX_FIX */ - - /************************************************************************ - (b) correct for the difference between this grid's fluxes and the - subgrid's fluxes. (step #19) - ************************************************************************/ - - TIME_MSG("Projecting solution to parent"); - for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { - EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); - - /* -------------- FIRST PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_POST_RECEIVE; - CommunicationReceiveIndex = 0; - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* Loop over subgrids for this grid: replace solution. */ - - CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - NextGrid = Grids[grid1]->NextGridNextLevel; - while (NextGrid != NULL) { - - /* Project the subgrid solution into this grid. */ - - NextGrid->GridData->ProjectSolutionToParentGrid(*Grids[grid1]->GridData); - NextGrid = NextGrid->NextGridThisLevel; - } // ENDWHILE subgrids - } // ENDFOR grids - - /* -------------- SECOND PASS ----------------- */ - - CommunicationDirection = COMMUNICATION_SEND; - for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { - - /* Loop over subgrids for this grid: replace solution. */ - - CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; - NextGrid = Grids[grid1]->NextGridNextLevel; - while (NextGrid != NULL) { - - /* Project the subgrid solution into this grid. */ - - NextGrid->GridData->ProjectSolutionToParentGrid(*Grids[grid1]->GridData); - NextGrid = NextGrid->NextGridThisLevel; - } // ENDWHILE subgrids - } // ENDFOR grids - - /* -------------- THIRD PASS ----------------- */ - - CommunicationReceiveHandler(); - - } // ENDFOR grid batches - - -#ifdef FORCE_MSG_PROGRESS - CommunicationBarrier(); -#endif - - CommunicationDirection = COMMUNICATION_SEND_RECEIVE; - - return SUCCESS; -} - - - -/* ======================================================================= */ -/* This routine simply converts a linked list of grids into an array of - pointers. */ - -int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, - HierarchyEntry **Grids[]) -{ - - /* Count the number of grids on this level. */ - - int NumberOfGrids = 0, counter = 0; - LevelHierarchyEntry *Temp = LevelArray[level]; - while (Temp != NULL) { - NumberOfGrids++; - Temp = Temp->NextGridThisLevel; - } - - /* Create a list of pointers and number of subgrids (and fill it out). */ - - typedef HierarchyEntry* HierarchyEntryPointer; - *Grids = new HierarchyEntryPointer[NumberOfGrids]; - Temp = LevelArray[level]; - while (Temp != NULL) { - (*Grids)[counter++] = Temp->GridHierarchyEntry; - Temp = Temp->NextGridThisLevel; - } - - return NumberOfGrids; -} - - diff --git a/src/enzo/PrepareGravitatingMassField.C b/src/enzo/PrepareGravitatingMassField.C index b310fa1b1..538ebfab0 100644 --- a/src/enzo/PrepareGravitatingMassField.C +++ b/src/enzo/PrepareGravitatingMassField.C @@ -101,7 +101,7 @@ int PrepareGravitatingMassField2(HierarchyEntry *Grid, TopGridData *MetaData, // if (CurrentGrid->AddBaryonsToGravitatingMassField() == FAIL) { //fprintf(stderr, " PGMF - DepositBaryons\n"); if (DepositBaryons(Grid, When) == FAIL) { - fprintf(stderr, "Error in grid->AddBaryonsToGravitatingMassField\n"); + fprintf(stderr, "Error in DepositBaryons\n"); ENZO_FAIL(""); } diff --git a/src/enzo/ProjectToPlane.C b/src/enzo/ProjectToPlane.C index bd969ef3a..e930633e5 100644 --- a/src/enzo/ProjectToPlane.C +++ b/src/enzo/ProjectToPlane.C @@ -18,11 +18,6 @@ #include #include #include - - - - - #include "ErrorExceptions.h" #include "macros_and_parameters.h" #include "typedefs.h" @@ -56,7 +51,6 @@ int CopyOverlappingParticleMassFields(grid* CurrentGrid, - int ProjectToPlane(TopGridData &MetaData, LevelHierarchyEntry *LevelArray[], int ProjectStartTemp[], int ProjectEndTemp[], FLOAT ProjectStartCoordinate[], diff --git a/src/enzo/ReadParameterFile.C b/src/enzo/ReadParameterFile.C index 4e49c73d3..4b1ecca20 100644 --- a/src/enzo/ReadParameterFile.C +++ b/src/enzo/ReadParameterFile.C @@ -830,32 +830,35 @@ int ReadParameterFile(FILE *fptr, TopGridData &MetaData, float *Initialdt) GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, &MassUnits, MetaData.Time); PressureUnits = DensityUnits*pow(VelocityUnits,2); - } - /* Change input physical parameters into code units */ - - double mh = 1.6726e-24; - double uheat = pow(VelocityUnits,2)*2.0*mh/TimeUnits; - PhotoelectricHeating /= uheat; - StarMakerOverDensityThreshold /= DensityUnits; - // StarEnergyFeedbackRate = StarEnergyFeedbackRate/pow(LengthUnits,2)*pow(TimeUnits,3); - - if (SinkMergeDistance > 1.0) - SinkMergeDistance /= LengthUnits; - //printf(" \n SinkMergeDistance = %"FSYM"\n \n", SinkMergeDistance); - SmallRho /= DensityUnits; - SmallP /= PressureUnits; - SmallT /= TemperatureUnits; - MaximumAlvenSpeed /= VelocityUnits; - float h, cs, dpdrho, dpde; - EOS(SmallP, SmallRho, SmallEint, h, cs, dpdrho, dpde, EOSType, 1); - if (debug && (HydroMethod == HD_RK || HydroMethod == MHD_RK)) - printf("smallrho=%g, smallp=%g, smalleint=%g, PressureUnits=%g, MaximumAlvenSpeed=%g\n", - SmallRho, SmallP, SmallEint, PressureUnits, MaximumAlvenSpeed); - for (int i = 0; i < MAX_FLAGGING_METHODS; i++) { - if (MinimumMassForRefinement[i] != FLOAT_UNDEFINED) { - MinimumMassForRefinement[i] /= MassUnits; + + /* Change input physical parameters into code units */ + + double mh = 1.6726e-24; + double uheat = pow(VelocityUnits,2)*2.0*mh/TimeUnits; + PhotoelectricHeating /= uheat; + StarMakerOverDensityThreshold /= DensityUnits; + // StarEnergyFeedbackRate = StarEnergyFeedbackRate/pow(LengthUnits,2)*pow(TimeUnits,3); + + if (SinkMergeDistance > 1.0) + SinkMergeDistance /= LengthUnits; + //printf(" \n SinkMergeDistance = %"FSYM"\n \n", SinkMergeDistance); + SmallRho /= DensityUnits; + SmallP /= PressureUnits; + SmallT /= TemperatureUnits; + MaximumAlvenSpeed /= VelocityUnits; + EOSSoundSpeed /= VelocityUnits; + float h, cs, dpdrho, dpde; + EOS(SmallP, SmallRho, SmallEint, h, cs, dpdrho, dpde, EOSType, 1); + if (debug && (HydroMethod == HD_RK || HydroMethod == MHD_RK)) + printf("smallrho=%g, smallp=%g, smalleint=%g, PressureUnits=%g, MaximumAlvenSpeed=%g\n", + SmallRho, SmallP, SmallEint, PressureUnits, MaximumAlvenSpeed); + for (int i = 0; i < MAX_FLAGGING_METHODS; i++) { + if (MinimumMassForRefinement[i] != FLOAT_UNDEFINED) { + MinimumMassForRefinement[i] /= MassUnits; + } } + } if (!ComovingCoordinates && UsePhysicalUnit) { diff --git a/src/enzo/SetBoundaryConditions.C b/src/enzo/SetBoundaryConditions.C new file mode 100644 index 000000000..f34f04e25 --- /dev/null +++ b/src/enzo/SetBoundaryConditions.C @@ -0,0 +1,215 @@ +/*********************************************************************** +/ +/ SET BOUNDARY CONDITIONS (CALLED BY EVOLVE LEVEL) +/ +/ written by: Greg Bryan +/ date: June, 1999 +/ modifiedN: Robert Harkness +/ date: February, 2008 +/ +/ PURPOSE: ======================================================================= +/ This routine sets all the boundary conditions for Grids by either +/ interpolating from their parents or copying from sibling grids. +/ +/ This is part of a collection of routines called by EvolveLevel. +/ These have been optimized for enhanced message passing +/ performance by performing two passes -- one which generates +/ sends and the second which receives them. +/ +/ modified: Robert Harkness, December 2007 +/ +************************************************************************/ + +#ifdef USE_MPI +#include +#endif /* USE_MPI */ + +#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" +#include "Hierarchy.h" +#include "TopGridData.h" +#include "LevelHierarchy.h" +#include "communication.h" +#include "CommunicationUtilities.h" + +/* function prototypes */ + +int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL, + int NumberOfSubgrids[] = NULL, + int FluxFlag = FALSE, + TopGridData* MetaData = NULL); + +#define GRIDS_PER_LOOP 20000 + + + +#ifdef FAST_SIB +int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, + SiblingGridList SiblingList[], + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, LevelHierarchyEntry *Level) +#else +int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, LevelHierarchyEntry *Level) +#endif +{ + + int loopEnd = (ShearingBoundaryDirection != 1) ? 2 : 1; + + int grid1, grid2, StartGrid, EndGrid, loop; + + JBPERF_START("SetBoundaryConditions"); + + for (loop = 0; loop < loopEnd; loop++){ + +#ifdef FORCE_MSG_PROGRESS + CommunicationBarrier(); +#endif + + TIME_MSG("Interpolating boundaries from parent"); + + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { + + if (traceMPI) fprintf(tracePtr, "SBC loop\n"); + + EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); + + /* -------------- FIRST PASS ----------------- */ + /* Here, we just generate the calls to generate the receive buffers, + without actually doing anything. */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; + CommunicationReceiveIndex = 0; + + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* a) Interpolate boundaries from the parent grid or set external + boundary conditions. */ + + CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; + + if (level == 0) { + if (loop == 0) + Grids[grid1]->GridData->SetExternalBoundaryValues(Exterior); + } else { + Grids[grid1]->GridData->InterpolateBoundaryFromParent + (Grids[grid1]->ParentGrid->GridData); + } + + } // ENDFOR grids + + /* -------------- SECOND PASS ----------------- */ + /* Now we generate all the sends, and do all the computation + for grids which are on the same processor as well. */ + + CommunicationDirection = COMMUNICATION_SEND; + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* a) Interpolate boundaries from the parent grid or set + external boundary conditions. */ + + if (level > 0) + Grids[grid1]->GridData->InterpolateBoundaryFromParent + (Grids[grid1]->ParentGrid->GridData); + + + } // ENDFOR grids + + //Grids[StartGrid]->GridData->PrintToScreenBoundaries(0); + + /* -------------- THIRD PASS ----------------- */ + /* In this final step, we get the messages as they come in and + then match them to the methods which generate the receive + handle. */ + + if (CommunicationReceiveHandler() == FAIL) + ENZO_FAIL(""); + + } // ENDFOR grid batches + + TIME_MSG("Copying zones in SetBoundaryConditions"); + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { + EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); + + /* -------------- FIRST PASS ----------------- */ + /* b) Copy any overlapping zones for sibling grids. */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; + CommunicationReceiveIndex = 0; + +#ifdef FAST_SIB + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) + for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++) + Grids[grid1]->GridData-> + CheckForOverlap(SiblingList[grid1].GridList[grid2], + MetaData->LeftFaceBoundaryCondition, + MetaData->RightFaceBoundaryCondition, + &grid::CopyZonesFromGrid); +#else + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) + for (grid2 = 0; grid2 < NumberOfGrids; grid2++) + Grids[grid1]->GridData-> + CheckForOverlap(Grids[grid2]->GridData, + MetaData->LeftFaceBoundaryCondition, + MetaData->RightFaceBoundaryCondition, + &grid::CopyZonesFromGrid); +#endif + + /* -------------- SECOND PASS ----------------- */ + /* b) Copy any overlapping zones for sibling grids. */ + + CommunicationDirection = COMMUNICATION_SEND; +#ifdef FAST_SIB + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) + for (grid2 = 0; grid2 < SiblingList[grid1].NumberOfSiblings; grid2++) + Grids[grid1]->GridData-> + CheckForOverlap(SiblingList[grid1].GridList[grid2], + MetaData->LeftFaceBoundaryCondition, + MetaData->RightFaceBoundaryCondition, + &grid::CopyZonesFromGrid); +#else + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) + for (grid2 = 0; grid2 < NumberOfGrids; grid2++) + Grids[grid1]->GridData-> + CheckForOverlap(Grids[grid2]->GridData, + MetaData->LeftFaceBoundaryCondition, + MetaData->RightFaceBoundaryCondition, + &grid::CopyZonesFromGrid); +#endif + + /* -------------- THIRD PASS ----------------- */ + + if (CommunicationReceiveHandler() == FAIL) + ENZO_FAIL(""); + + } // end loop over batchs of grids + + } // ENDFOR loop (for ShearinBox) + + /* c) Apply external reflecting boundary conditions, if needed. */ + + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) + Grids[grid1]->GridData->CheckForExternalReflections + (MetaData->LeftFaceBoundaryCondition, + MetaData->RightFaceBoundaryCondition); + +#ifdef FORCE_MSG_PROGRESS + CommunicationBarrier(); +#endif + + CommunicationDirection = COMMUNICATION_SEND_RECEIVE; + + JBPERF_STOP("SetBoundaryConditions"); + + return SUCCESS; + +} diff --git a/src/enzo/UpdateFromFinerGrids.C b/src/enzo/UpdateFromFinerGrids.C new file mode 100644 index 000000000..714b0bffa --- /dev/null +++ b/src/enzo/UpdateFromFinerGrids.C @@ -0,0 +1,328 @@ +/*********************************************************************** +/ +/ UPDATE FROM FINER GRIDS (CALLED BY EVOLVE LEVEL) +/ +/ written by: Greg Bryan +/ date: June, 1999 +/ modifiedN: Robert Harkness +/ date: February, 2008 +/ +/ PURPOSE: +/ ======================================================================= +/ This routines does the flux correction and project for all grids on this +/ level from the list of subgrids. +/ +/ This is part of a collection of routines called by EvolveLevel. +/ These have been optimized for enhanced message passing +/ performance by performing two passes -- one which generates +/ sends and the second which receives them. +/ +/ modified: Robert Harkness, December 2007 +/ +************************************************************************/ + +#ifdef USE_MPI +#include +#endif /* USE_MPI */ + +#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" +#include "Hierarchy.h" +#include "TopGridData.h" +#include "LevelHierarchy.h" +#include "communication.h" +#include "CommunicationUtilities.h" + +/* function prototypes */ + +int CommunicationReceiveHandler(fluxes **SubgridFluxesEstimate[] = NULL, + int NumberOfSubgrids[] = NULL, + int FluxFlag = FALSE, + TopGridData* MetaData = NULL); + +#define GRIDS_PER_LOOP 20000 + + +#ifdef FLUX_FIX +int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids, + int NumberOfSubgrids[], + fluxes **SubgridFluxesEstimate[], + LevelHierarchyEntry* SUBlingList[], + TopGridData *MetaData) +#else +int UpdateFromFinerGrids(int level, HierarchyEntry *Grids[], int NumberOfGrids, + int NumberOfSubgrids[], + fluxes **SubgridFluxesEstimate[]) +#endif + +{ + + int grid1, subgrid, StartGrid, EndGrid; + HierarchyEntry *NextGrid; + +#ifdef FLUX_FIX + int SUBlingGrid; + LevelHierarchyEntry *NextEntry; +#endif + + /* Define a temporary flux holder for the refined fluxes. */ + + fluxes SubgridFluxesRefined; + + /* For each grid, + (a) project the subgrid's solution into this grid (step #18) + (a1) project the neighboring subgrid's solution into this grid + (b) correct for the difference between this grid's fluxes and the + subgrid's fluxes. (step #19) */ + +#ifdef FORCE_MSG_PROGRESS + CommunicationBarrier(); +#endif + + TIME_MSG("UpdateFromFinerGrids"); + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { + EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); + + /* -------------- FIRST PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; + CommunicationReceiveIndex = 0; + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* Loop over subgrids for this grid. */ + + NextGrid = Grids[grid1]->NextGridNextLevel; + subgrid = 0; + CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; + + while (NextGrid != NULL && FluxCorrection) { + + /* Project subgrid's refined fluxes to the level of this grid. */ + +#ifdef USE_MPI + CommunicationReceiveArgumentInt[0][CommunicationReceiveIndex] = grid1; + CommunicationReceiveArgumentInt[1][CommunicationReceiveIndex] = subgrid; +#endif /* USE_MPI */ + + NextGrid->GridData-> + GetProjectedBoundaryFluxes(Grids[grid1]->GridData, SubgridFluxesRefined); + + NextGrid = NextGrid->NextGridThisLevel; + subgrid++; + } // ENDWHILE subgrids + + } // ENDFOR grids + + /* -------------- SECOND PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_SEND; + + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* Loop over subgrids for this grid. */ + + NextGrid = Grids[grid1]->NextGridNextLevel; + subgrid = 0; + while (NextGrid != NULL && FluxCorrection) { + + /* Project subgrid's refined fluxes to the level of this grid. */ + + NextGrid->GridData-> + GetProjectedBoundaryFluxes(Grids[grid1]->GridData, SubgridFluxesRefined); + + /* Correct this grid for the refined fluxes (step #19) + (this also deletes the fields in SubgridFluxesRefined). + (only call it if the grid and sub-grid are on the same + processor, otherwise handled in CommunicationReceiveHandler.) */ + +#ifdef FLUX_FIX + if (NextGrid->GridData->ReturnProcessorNumber() == + Grids[grid1]->GridData->ReturnProcessorNumber()) + Grids[grid1]->GridData->CorrectForRefinedFluxes + (SubgridFluxesEstimate[grid1][subgrid], &SubgridFluxesRefined, + SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], + FALSE, MetaData); +#else + if (NextGrid->GridData->ReturnProcessorNumber() == + Grids[grid1]->GridData->ReturnProcessorNumber()) + Grids[grid1]->GridData->CorrectForRefinedFluxes + (SubgridFluxesEstimate[grid1][subgrid], &SubgridFluxesRefined, + SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1]); +#endif /* FLUX_FIX */ + + NextGrid = NextGrid->NextGridThisLevel; + subgrid++; + + } // ENDWHILE subgrids + + } // ENDFOR grids + + /* -------------- THIRD PASS ----------------- */ + +#ifdef FLUX_FIX + CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids, + FALSE, MetaData); +#else + CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids); +#endif + + } // ENDFOR grid batches + + /************************************************************************ + (a1) project the neighboring subgrid's solution into this grid + (SUBlings) + ************************************************************************/ + +#ifdef FLUX_FIX + TIME_MSG("Projecting neighboring subgrid solution (FLUX_FIX)"); + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { + EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); + + /* -------------- FIRST PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; + CommunicationReceiveIndex = 0; + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* Loop over subgrids for this grid. */ + + CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; + NextEntry = SUBlingList[grid1]; + + while (NextEntry != NULL && FluxCorrection) { + + /* make sure this isn't a "proper" subgrid */ + + if (NextEntry->GridHierarchyEntry->ParentGrid != Grids[grid1]) { + + /* Project subgrid's refined fluxes to the level of this grid. */ + +#ifdef USE_MPI + CommunicationReceiveArgumentInt[0][CommunicationReceiveIndex] = grid1; + CommunicationReceiveArgumentInt[1][CommunicationReceiveIndex] = + NumberOfSubgrids[grid1]-1; +#endif /* USE_MPI */ + + NextEntry->GridData->GetProjectedBoundaryFluxes + (Grids[grid1]->GridData, SubgridFluxesRefined); + } + + NextEntry = NextEntry->NextGridThisLevel; + + } // ENDWHILE subgrids + } // ENDFOR grids + + /* -------------- SECOND PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_SEND; + + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + NextEntry = SUBlingList[grid1]; + + while (NextEntry != NULL && FluxCorrection) { + + /* make sure this isn't a "proper" subgrid */ + + if (NextEntry->GridHierarchyEntry->ParentGrid != Grids[grid1]) { + + /* Project subgrid's refined fluxes to the level of this grid. */ + + NextEntry->GridData->GetProjectedBoundaryFluxes + (Grids[grid1]->GridData, SubgridFluxesRefined); + + /* Correct this grid for the refined fluxes (step #19) + (this also deletes the fields in SubgridFluxesRefined). */ + + if (NextEntry->GridData->ReturnProcessorNumber() == + Grids[grid1]->GridData->ReturnProcessorNumber()) + Grids[grid1]->GridData->CorrectForRefinedFluxes + (SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], + &SubgridFluxesRefined, + SubgridFluxesEstimate[grid1][NumberOfSubgrids[grid1] - 1], + TRUE, MetaData); + } + + NextEntry = NextEntry->NextGridThisLevel; + + } // ENDWHILE subgrids + } // ENDFOR grids + + /* -------------- THIRD PASS ----------------- */ + + CommunicationReceiveHandler(SubgridFluxesEstimate, NumberOfSubgrids, + TRUE, MetaData); + + } // ENDFOR grid batches +#endif /* FLUX_FIX */ + + /************************************************************************ + (b) correct for the difference between this grid's fluxes and the + subgrid's fluxes. (step #19) + ************************************************************************/ + + TIME_MSG("Projecting solution to parent"); + for (StartGrid = 0; StartGrid < NumberOfGrids; StartGrid += GRIDS_PER_LOOP) { + EndGrid = min(StartGrid + GRIDS_PER_LOOP, NumberOfGrids); + + /* -------------- FIRST PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_POST_RECEIVE; + CommunicationReceiveIndex = 0; + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* Loop over subgrids for this grid: replace solution. */ + + CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; + NextGrid = Grids[grid1]->NextGridNextLevel; + while (NextGrid != NULL) { + + /* Project the subgrid solution into this grid. */ + + NextGrid->GridData->ProjectSolutionToParentGrid(*Grids[grid1]->GridData); + NextGrid = NextGrid->NextGridThisLevel; + } // ENDWHILE subgrids + } // ENDFOR grids + + /* -------------- SECOND PASS ----------------- */ + + CommunicationDirection = COMMUNICATION_SEND; + for (grid1 = StartGrid; grid1 < EndGrid; grid1++) { + + /* Loop over subgrids for this grid: replace solution. */ + + CommunicationReceiveCurrentDependsOn = COMMUNICATION_NO_DEPENDENCE; + NextGrid = Grids[grid1]->NextGridNextLevel; + while (NextGrid != NULL) { + + /* Project the subgrid solution into this grid. */ + + NextGrid->GridData->ProjectSolutionToParentGrid(*Grids[grid1]->GridData); + NextGrid = NextGrid->NextGridThisLevel; + } // ENDWHILE subgrids + } // ENDFOR grids + + /* -------------- THIRD PASS ----------------- */ + + CommunicationReceiveHandler(); + + } // ENDFOR grid batches + + +#ifdef FORCE_MSG_PROGRESS + CommunicationBarrier(); +#endif + + CommunicationDirection = COMMUNICATION_SEND_RECEIVE; + + return SUCCESS; +} + diff --git a/src/enzo/hydro_rk/Collapse3DInitialize.C b/src/enzo/hydro_rk/Collapse3DInitialize.C index 84d6ad5d5..2a171eaff 100644 --- a/src/enzo/hydro_rk/Collapse3DInitialize.C +++ b/src/enzo/hydro_rk/Collapse3DInitialize.C @@ -132,10 +132,10 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr, ret += sscanf(line, "SphereType[%d] = %d", &sphere, &SphereType[sphere]); if (sscanf(line, "SphereRadius[%d]", &sphere) > 0) - ret += sscanf(line, "SphereRadius[%d] = %"FSYM, &sphere, + ret += sscanf(line, "SphereRadius[%d] = %"PSYM, &sphere, &SphereRadius[sphere]); if (sscanf(line, "SphereCoreRadius[%d]", &sphere) > 0) - ret += sscanf(line, "SphereCoreRadius[%d] = %"FSYM, &sphere, + ret += sscanf(line, "SphereCoreRadius[%d] = %"PSYM, &sphere, &SphereCoreRadius[sphere]); if (sscanf(line, "SphereDensity[%d]", &sphere) > 0) ret += sscanf(line, "SphereDensity[%d] = %f", &sphere, @@ -147,7 +147,7 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr, ret += sscanf(line, "SphereSoundVelocity[%d] = %f", &sphere, &SphereSoundVelocity[sphere]); if (sscanf(line, "SpherePosition[%d]", &sphere) > 0) - ret += sscanf(line, "SpherePosition[%d] = %"FSYM" %"FSYM" %"FSYM, + ret += sscanf(line, "SpherePosition[%d] = %"PSYM" %"PSYM" %"PSYM, &sphere, &SpherePosition[sphere][0], &SpherePosition[sphere][1], &SpherePosition[sphere][2]); @@ -354,25 +354,25 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr, /* Write parameters to parameter output file */ - /*if (MyProcessorNumber == ROOT_PROCESSOR) { + if (MyProcessorNumber == ROOT_PROCESSOR) { fprintf(Outfptr, "NumberOfSpheres = %d\n", n_sphere); fprintf(Outfptr, "RefineAtStart = %d\n", RefineAtStart); fprintf(Outfptr, "UseParticles = %d\n", UseParticles); - fprintf(Outfptr, "UseColour = %d\n", - UseColour); - fprintf(Outfptr, "InitialTemperature = %f\n", - InitialTemperature); + // fprintf(Outfptr, "UseColour = %d\n", + // UseColour); + // fprintf(Outfptr, "InitialTemperature = %f\n", + // InitialTemperature); fprintf(Outfptr, "UniformVelocity = %f %f %f\n", UniformVelocity[0], UniformVelocity[1], UniformVelocity[2]); fprintf(Outfptr, "LengthUnit = %f\n", - LengthUnit); + lenu); fprintf(Outfptr, "DensityUnit = %f\n", - DensityUnit); - for (sphere = 0; sphere < NumberOfSpheres; sphere++) { + rhou); + for (sphere = 0; sphere < n_sphere; sphere++) { fprintf(Outfptr, "SphereType[%d] = %d\n", sphere, SphereType[sphere]); fprintf(Outfptr, "SphereRadius[%d] = %"GOUTSYM"\n", sphere, @@ -381,16 +381,16 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr, SphereCoreRadius[sphere]); fprintf(Outfptr, "SphereDensity[%d] = %f\n", sphere, SphereDensity[sphere]); - fprintf(Outfptr, "SphereTemperature[%d] = %f\n", sphere, - SphereTemperature[sphere]); + fprintf(Outfptr, "SphereSoundVelocity[%d] = %f\n", sphere, + SphereSoundVelocity[sphere]); fprintf(Outfptr, "SpherePosition[%d] = ", sphere); WriteListOfFloats(Outfptr, MetaData.TopGridRank, SpherePosition[sphere]); fprintf(Outfptr, "SphereVelocity[%d] = ", sphere); WriteListOfFloats(Outfptr, MetaData.TopGridRank, SphereVelocity[sphere]); - fprintf(Outfptr, "FracKeplarianRot[%d] = %"GOUTSYM"\n", sphere, - FracKeplarianRot[sphere]); + // fprintf(Outfptr, "FracKeplarianRot[%d] = %"GOUTSYM"\n", sphere, + // FracKeplarianRot[sphere]); fprintf(Outfptr, "SphereTurbulence[%d] = %"GOUTSYM"\n", sphere, SphereTurbulence[sphere]); fprintf(Outfptr, "SphereCutOff[%d] = %"GOUTSYM"\n", sphere, @@ -402,7 +402,7 @@ int Collapse3DInitialize(FILE *fptr, FILE *Outfptr, fprintf(Outfptr, "SphereNumShells[%d] = %d\n", sphere, SphereNumShells[sphere]); } - }*/ + } return SUCCESS; diff --git a/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C b/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C new file mode 100644 index 000000000..c66dcecf0 --- /dev/null +++ b/src/enzo/hydro_rk/ComputeDednerWaveSpeeds.C @@ -0,0 +1,77 @@ +/*********************************************************************** +/ +/ COMPUTE DEDNER WAVE SPEEDS +/ +/ written by: Tom Abel +/ date: October 2009 +/ +/ ======================================================================= +/ This routine computes the wave speeds used for the Dedner formalism. +/ +/ Reference: e.g. Matsumoto, PASJ, 2007, 59, 905 +/ +/ Output: C_h and C_p global variables defined in global_data.h +/ +************************************************************************/ + +#include +#include +#include "macros_and_parameters.h" +#include "typedefs.h" +#include "global_data.h" +#include "Fluxes.h" +#include "GridList.h" +#include "ExternalBoundary.h" +#include "TopGridData.h" +#include "Grid.h" +#include "Hierarchy.h" +#include "LevelHierarchy.h" +#include "../hydro_rk/tools.h" + +int GetUnits(float *DensityUnits, float *LengthUnits, + float *TemperatureUnits, float *TimeUnits, + float *VelocityUnits, FLOAT Time); + +int ComputeDednerWaveSpeeds(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], + int level, FLOAT dt0) +{ + /* Count the number of grids on this level. */ + + if (HydroMethod != MHD_RK) + return SUCCESS; + + float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, + VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits; + double MassUnits; + GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, 1.0); + MassUnits = DensityUnits*pow(LengthUnits,3); + + int lmax; + LevelHierarchyEntry *Temp; + for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) { + Temp = LevelArray[lmax]; + if (Temp != NULL) + break; + } + + // lmax = 0; // <- Pengs version had lmax = 6 + FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0; + + dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0]; + dy0 = (MetaData->TopGridRank > 1) ? + (DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8; + dz0 = (MetaData->TopGridRank > 2) ? + (DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8; + h_min = my_MIN(dx0, dy0, dz0); + h_min /= pow(RefineBy, lmax); + C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0); + C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems) + if (EOSType == 3) // for isothermal runs just use the constant sound speed + C_h = EOSSoundSpeed; + + C_p = sqrt(0.18*DivBDampingLength*C_h); + + return SUCCESS; +} + + diff --git a/src/enzo/hydro_rk/EOS.h b/src/enzo/hydro_rk/EOS.h index ef4e8c4fc..287258236 100644 --- a/src/enzo/hydro_rk/EOS.h +++ b/src/enzo/hydro_rk/EOS.h @@ -42,7 +42,7 @@ inline void EOS(float &p, float &rho, float &e, float &h, float &cs, float &dpdr GetUnits(&denu, &lenu, &tempu, &tu, &velu, 1); double c_s = EOSSoundSpeed; double rho_cr = EOSCriticalDensity; - c_s /= velu; + // c_s /= velu; rho_cr /= denu; cs = c_s*sqrt(1.0 + EOSGamma*pow(rho/rho_cr, EOSGamma-1.0)); @@ -59,7 +59,7 @@ inline void EOS(float &p, float &rho, float &e, float &h, float &cs, float &dpdr GetUnits(&denu, &lenu, &tempu, &tu, &velu, 1); double c_s = EOSSoundSpeed; double rho_cr = EOSCriticalDensity; - c_s /= velu; + // c_s /= velu; rho_cr /= denu; if (rho <= rho_cr) { diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C index 92a4924eb..1e91c4f64 100644 --- a/src/enzo/hydro_rk/EvolveLevel_RK2.C +++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C @@ -2,7 +2,7 @@ / / EVOLVE LEVEL USING RUNGE-KUTTA2 METHOD / -/ written by: Peng Wang +/ written by: Peng Wang & Tom Abel / date: May, 2007 / modified1: / @@ -40,7 +40,10 @@ int StarParticleInitialize(LevelHierarchyEntry *LevelArray[], int ThisLevel, int StarParticleFinalize(HierarchyEntry *Grids[], TopGridData *MetaData, int NumberOfGrids, LevelHierarchyEntry *LevelArray[], int level, Star *&AllStars); - +int AdjustRefineRegion(LevelHierarchyEntry *LevelArray[], + TopGridData *MetaData, int EL_level); +int ComputeDednerWaveSpeeds(TopGridData *MetaData,LevelHierarchyEntry *LevelArray[], + int level, FLOAT dt0); #ifdef TRANSFER int EvolvePhotons(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], Star *AllStars, FLOAT GridTime, int level, int LoopTime = TRUE); @@ -69,13 +72,28 @@ float CommunicationMaxValue(float Value); int CommunicationBarrier(); int GenerateGridArray(LevelHierarchyEntry *LevelArray[], int level, HierarchyEntry **Grids[]); +int CallProblemSpecificRoutines(TopGridData * MetaData, HierarchyEntry *ThisGrid, + int GridNum, float *norm, float TopGridTimeStep, + int level, int LevelCycleCount[]); //moo + #ifdef FAST_SIB int PrepareDensityField(LevelHierarchyEntry *LevelArray[], SiblingGridList SiblingList[], int level, TopGridData *MetaData, FLOAT When); +int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids, + SiblingGridList SiblingList[], + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, + LevelHierarchyEntry * Level, + int CycleNumber); #else // !FAST_SIB int PrepareDensityField(LevelHierarchyEntry *LevelArray[], int level, TopGridData *MetaData, FLOAT When); +int SetAccelerationBoundary(HierarchyEntry *Grids[], int NumberOfGrids, + int level, TopGridData *MetaData, + ExternalBoundary *Exterior, + LevelHierarchyEntry * Level, + int CycleNumber); #endif // end FAST_SIB int SetBoundaryConditions(HierarchyEntry *Grids[], int NumberOfGrids, @@ -130,8 +148,8 @@ int FinalizeFluxes(HierarchyEntry *Grids[],fluxes **SubgridFluxesEstimate[], int NumberOfGrids,int NumberOfSubgrids[]); int RadiationFieldUpdate(LevelHierarchyEntry *LevelArray[], int level, TopGridData *MetaData); -//int WriteStreamData(LevelHierarchyEntry *LevelArray[], int level, -// TopGridData *MetaData, int *CycleCount, int open=FALSE); +int WriteStreamData(LevelHierarchyEntry *LevelArray[], int level, + TopGridData *MetaData, int *CycleCount, int open=FALSE); //int WriteMovieData(char *basename, int filenumber, // LevelHierarchyEntry *LevelArray[], TopGridData *MetaData, // FLOAT WriteTime); @@ -144,9 +162,6 @@ int FastSiblingLocatorFinalize(ChainingMeshStructure *Mesh); int FastSiblingLocatorInitializeStaticChainingMesh(ChainingMeshStructure *Mesh, int Rank, int TopGridDims[]); -int GetUnits(float *DensityUnits, float *LengthUnits, - float *TemperatureUnits, float *TimeUnits, - float *VelocityUnits, FLOAT Time); double ReturnWallTime(); int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData, int level); @@ -162,9 +177,9 @@ int CallPython(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData, static int LevelCycleCount[MAX_DEPTH_OF_HIERARCHY]; static int MovieCycleCount[MAX_DEPTH_OF_HIERARCHY]; -//double LevelWallTime[MAX_DEPTH_OF_HIERARCHY]; -//double LevelZoneCycleCount[MAX_DEPTH_OF_HIERARCHY]; -//double LevelZoneCycleCountPerProc[MAX_DEPTH_OF_HIERARCHY]; +static double LevelWallTime[MAX_DEPTH_OF_HIERARCHY]; +static double LevelZoneCycleCount[MAX_DEPTH_OF_HIERARCHY]; +static double LevelZoneCycleCountPerProc[MAX_DEPTH_OF_HIERARCHY]; static float norm = 0.0; //AK static float TopGridTimeStep = 0.0; //AK @@ -178,6 +193,8 @@ static int StaticLevelZero = 1; static int StaticLevelZero = 0; #endif +extern int RK2SecondStepBaryonDeposit; + /* EvolveGrid function */ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], @@ -203,12 +220,7 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], #endif FLOAT When; - float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1, TimeUnits, - VelocityUnits, CriticalDensity = 1, BoxLength = 1, MagneticUnits; - double MassUnits; - GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, - &TimeUnits, &VelocityUnits, 1.0); - MassUnits = DensityUnits*pow(LengthUnits,3); + /* Create an array (Grids) of all the grids. */ @@ -218,14 +230,23 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], int *NumberOfSubgrids = new int[NumberOfGrids]; fluxes ***SubgridFluxesEstimate = new fluxes **[NumberOfGrids]; + /* Initialize the chaining mesh used in the FastSiblingLocator. */ + + SiblingGridList *SiblingList = new SiblingGridList[NumberOfGrids]; CreateSiblingList(Grids, NumberOfGrids, SiblingList, StaticLevelZero, MetaData, level); - /*if (SetBoundaryConditions(Grids, NumberOfGrids,SiblingList,level, MetaData, - Exterior) == FAIL) { - ENZO_FAIL(""); - }*/ + /* On the top grid, adjust the refine region so that only the finest + particles are included. We don't want the more massive particles + to contaminate the high-resolution region. */ + + AdjustRefineRegion(LevelArray, MetaData, level); + + /* ================================================================== */ + /* For each grid: a) interpolate boundaries from its parent. + b) copy any overlapping zones. */ + #ifdef FAST_SIB if (SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, MetaData, Exterior, LevelArray[level]) == FAIL) @@ -235,17 +256,24 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], Exterior, LevelArray[level]) == FAIL) ENZO_FAIL(""); #endif - - - + /* Count the number of colours in the first grid (to define NColor) */ Grids[0]->GridData->SetNumberOfColours(); //fprintf(stdout, "EvolveLevel_RK2: NColor =%d, NSpecies = %d\n", NColor, NSpecies); + /* Clear the boundary fluxes for all Grids (this will be accumulated over + the subcycles below (i.e. during one current grid step) and used to by the + current grid to correct the zones surrounding this subgrid (step #18). */ + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) Grids[grid1]->GridData->ClearBoundaryFluxes(); + /* After we calculate the ghost zones, we can initialize streaming + data files (only on level 0) */ + + if (MetaData->FirstTimestepAfterRestart == TRUE && level == 0) + WriteStreamData(LevelArray, level, MetaData, MovieCycleCount); /* ================================================================== */ /* Loop over grid timesteps until the elapsed time equals the timestep @@ -256,6 +284,11 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], SetLevelTimeStep(Grids, NumberOfGrids, level, &dtThisLevelSoFar, &dtThisLevel, dtLevelAbove); + /* Streaming movie output (write after all parent grids are + updated) */ + + WriteStreamData(LevelArray, level, MetaData, MovieCycleCount); + /* Initialize the star particles */ Star *AllStars = NULL; @@ -266,7 +299,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], RadiativeTransferPrepare(LevelArray, level, MetaData, AllStars, dtLevelAbove); #endif /* TRANSFER */ - /* For each grid, compute the number of it's subgrids. */ for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { @@ -325,173 +357,131 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], } // end loop over grids (create Subgrid list) - /* compute wave speed - Reference: Matsumoto, PASJ, 2007, 59, 905 */ - - if (HydroMethod == MHD_RK) { - int lmax; - LevelHierarchyEntry *Temp; - for (lmax = MAX_DEPTH_OF_HIERARCHY-1; lmax >= 0; lmax--) { - Temp = LevelArray[lmax]; - if (Temp != NULL) { - break; - } - } - // lmax = 0; // <- Pengs version had lmax = 6 - // lmax = 1; - FLOAT dx0, dy0, dz0, h_min, DivBDampingLength = 1.0; - - dx0 = (DomainRightEdge[0] - DomainLeftEdge[0]) / MetaData->TopGridDims[0]; - dy0 = (MetaData->TopGridRank > 1) ? - (DomainRightEdge[1] - DomainLeftEdge[1]) / MetaData->TopGridDims[1] : 1e8; - dz0 = (MetaData->TopGridRank > 2) ? - (DomainRightEdge[2] - DomainLeftEdge[2]) / MetaData->TopGridDims[2] : 1e8; - h_min = my_MIN(dx0, dy0, dz0); - h_min /= pow(RefineBy, lmax); - C_h = 0.5*MetaData->CourantSafetyNumber*(h_min/dt0); - C_h = min( C_h, 1e6/VelocityUnits); // never faster than __ cm/s (for very small dt0 a problems) - if (EOSType == 3) // for isothermal runs just use the constant sound speed - C_h = EOSSoundSpeed; - - C_p = sqrt(0.18*DivBDampingLength*C_h); - // C_p = sqrt(0.18*DivBDampingLength)*C_h; - if (debug) - fprintf(stderr, "lengthscale %g timestep: %g C_h: %g C_p: %g\n ", - h_min, dt0, C_h, C_p); - } + /* compute wave speed Reference: Matsumoto, PASJ, 2007, 59, 905 */ -// if (SelfGravity && MetaData->TopGridRank == 3) { -// if (PrepareDensityField(LevelArray, SiblingList, level, MetaData) == FAIL) { -// // if (PrepareDensityField(LevelArray, level, MetaData) == FAIL) { -// fprintf(stderr, "Error in PrepareDensityField.\n"); -// ENZO_FAIL(""); -// } -// } + if (HydroMethod == MHD_RK) + ComputeDednerWaveSpeeds(MetaData, LevelArray, level, dt0); + + if (debug && HydroMethod == MHD_RK) + fprintf(stderr, "wave speeds: timestep: %g C_h: %g C_p: %g\n ", + dt0, C_h, C_p); - When = 0.5; + When = 0.5; + RK2SecondStepBaryonDeposit = 0; + if (SelfGravity) { #ifdef FAST_SIB - if (SelfGravity) - if (PrepareDensityField(LevelArray, SiblingList, - level, MetaData, When) == FAIL) { - fprintf(stderr, "Error in PrepareDensityField.\n"); - ENZO_FAIL(""); - } + PrepareDensityField(LevelArray, SiblingList, level, MetaData, When); #else // !FAST_SIB - if (SelfGravity) - if (PrepareDensityField(LevelArray, level, MetaData, When) == FAIL) { - fprintf(stderr, "Error in PrepareDensityField.\n"); - ENZO_FAIL(""); - } + PrepareDensityField(LevelArray, level, MetaData, When); #endif // end FAST_SIB + } /* Solve the radiative transfer */ #ifdef TRANSFER - //Grids[grid1]->GridData->SetTimePreviousTimestep(); FLOAT GridTime = Grids[0]->GridData->ReturnTime(); EvolvePhotons(MetaData, LevelArray, AllStars, GridTime, level); - //Grids[grid1]->GridData->SetTimeNextTimestep(); #endif /* TRANSFER */ /* Compute particle-particle acceleration */ - //if (NBodyDirectSummation == TRUE) - //for (grid = 0; grid < NumberOfGrids; grid++) - //Grids[grid1]->GridData->ComputeParticleParticleAcceleration(level); + /* if (NBodyDirectSummation == TRUE) + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) + Grids[grid1]->GridData->ComputeParticleParticleAcceleration(level); */ /* ------------------------------------------------------- */ - /* Evolve all grids by timestep dtThisLevel. */ + /* Evolve all grids by timestep dtThisLevel. Predictor Step*/ for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { + + CallProblemSpecificRoutines(MetaData, Grids[grid1], grid1, &norm, + TopGridTimeStep, level, LevelCycleCount); /* Gravity: compute acceleration field for grid and particles. */ - if (SelfGravity) { int Dummy; if (level <= MaximumGravityRefinementLevel) { if (level > 0) Grids[grid1]->GridData->SolveForPotential(level) ; - Grids[grid1]->GridData->ComputeAccelerations(level) ; } - // otherwise, interpolate potential from coarser grid, which is - // now done in PrepareDensity. + // otherwise, use interpolated potential from coarser grid, which is + // done in PrepareDensity. } // end: if (SelfGravity) Grids[grid1]->GridData->ComputeAccelerationFieldExternal() ; - #ifdef TRANSFER - /* Radiation Pressure: add to acceleration field */ - - if (RadiativeTransfer && RadiationPressure) - if (Grids[grid1]->GridData->AddRadiationPressureAcceleration() == FAIL) { - fprintf(stderr,"Error in grid->AddRadiationPressureAcceleration.\n"); - ENZO_FAIL(""); - } - + Grids[grid1]->GridData->AddRadiationPressureAcceleration(); #endif /* TRANSFER */ - Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField(); - if (UseHydro) { - + if (UseHydro) if (HydroMethod == HD_RK) Grids[grid1]->GridData->RungeKutta2_1stStep - (SubgridFluxesEstimate[grid1], - NumberOfSubgrids[grid1], level, Exterior); - - else if (HydroMethod == MHD_RK) { + (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior); + else if (HydroMethod == MHD_RK) Grids[grid1]->GridData->MHDRK2_1stStep - (SubgridFluxesEstimate[grid1], - NumberOfSubgrids[grid1], level, Exterior); - } - } // ENDIF UseHydro + (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior); - /* Do this here so that we can get the correct - time interpolated boundary condition */ + /* Do this here so that we can get the correct time interpolated boundary condition */ Grids[grid1]->GridData->SetTimeNextTimestep(); } // end loop over grids - /*if (SetBoundaryConditions(Grids, NumberOfGrids,SiblingList,level, MetaData, - Exterior) == FAIL) { - ENZO_FAIL(""); - }*/ + #ifdef FAST_SIB - SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, - level, MetaData, Exterior, LevelArray[level]); + SetBoundaryConditions(Grids, NumberOfGrids, SiblingList, level, MetaData, Exterior, LevelArray[level]); #else - SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData, - Exterior, LevelArray[level]); + SetBoundaryConditions(Grids, NumberOfGrids, level, MetaData, Exterior, LevelArray[level]); #endif + // Recompute potential and accelerations with time centered baryon Field + // this also does the particles again at the moment so could be made more efficient. + + RK2SecondStepBaryonDeposit = 0; // set this to (0/1) to (not use/use) this extra step + // printf("SECOND STEP\n"); + if (RK2SecondStepBaryonDeposit & SelfGravity && UseHydro) { + When = 0.5; +#ifdef FAST_SIB + PrepareDensityField(LevelArray, SiblingList, level, MetaData, When); +#else + PrepareDensityField(LevelArray, level, MetaData, When); +#endif // end FAST_SIB + } + for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { + /* Gravity: compute acceleration field for grid and particles. */ + if (RK2SecondStepBaryonDeposit) { + int Dummy; + if (level <= MaximumGravityRefinementLevel) { + if (level > 0) + Grids[grid1]->GridData->SolveForPotential(level) ; + Grids[grid1]->GridData->ComputeAccelerations(level) ; + } + Grids[grid1]->GridData->ComputeAccelerationFieldExternal() ; + } // end: if (SelfGravity) + if (UseHydro) { if (HydroMethod == HD_RK) Grids[grid1]->GridData->RungeKutta2_2ndStep - (SubgridFluxesEstimate[grid1], - NumberOfSubgrids[grid1], level, Exterior); + (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior); else if (HydroMethod == MHD_RK) { - + Grids[grid1]->GridData->MHDRK2_2ndStep - (SubgridFluxesEstimate[grid1], - NumberOfSubgrids[grid1], level, Exterior); + (SubgridFluxesEstimate[grid1], NumberOfSubgrids[grid1], level, Exterior); - if (UseAmbipolarDiffusion) { + if (UseAmbipolarDiffusion) Grids[grid1]->GridData->AddAmbipolarDiffusion(); - } - - if (UseResistivity) { + + if (UseResistivity) Grids[grid1]->GridData->AddResistivity(); - } - time1 = ReturnWallTime(); Grids[grid1]->GridData->PoissonSolver(level); @@ -501,9 +491,8 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], /* Add viscosity */ - if (UseViscosity) { + if (UseViscosity) Grids[grid1]->GridData->AddViscosity(); - } /* Solve the cooling and species rate equations. */ @@ -528,7 +517,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], if (UseFloor) Grids[grid1]->GridData->SetFloor(); - /* If using comoving co-ordinates, do the expansion terms now. */ @@ -579,7 +567,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], #ifdef USE_JBPERF // Update lcaperf "level" attribute - jbPerf.attribute ("level",&jb_level,JB_INT); #endif @@ -650,9 +637,9 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], for (grid1 = 0; grid1 < NumberOfGrids; grid1++) { Grids[grid1]->GridData->CollectGridInformation (GridMemory, GridVolume, NumberOfCells, AxialRatio, CellsTotal, Particles); - // LevelZoneCycleCount[level] += NumberOfCells; - //if (MyProcessorNumber == Grids[grid1]->GridData->ReturnProcessorNumber()) - // LevelZoneCycleCountPerProc[level] += NumberOfCells; + LevelZoneCycleCount[level] += NumberOfCells; + if (MyProcessorNumber == Grids[grid1]->GridData->ReturnProcessorNumber()) + LevelZoneCycleCountPerProc[level] += NumberOfCells; } #endif @@ -672,12 +659,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], /* Clean up. */ -#ifdef UNUSED - if (level > MaximumGravityRefinementLevel && - level == MaximumRefinementLevel) - ZEUSQuadraticArtificialViscosity /= 1; -#endif - delete [] NumberOfSubgrids; delete [] Grids; delete [] SubgridFluxesEstimate; diff --git a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C index d674e59fc..9adb96b29 100644 --- a/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_Collapse3DInitializeGrid.C @@ -214,7 +214,7 @@ int grid::Collapse3DInitializeGrid(int n_sphere, eint = pow(cs_sphere[sphere], 2)/(Gamma-1.0); vel[0] = -omega_sphere[sphere]*ypos; vel[1] = omega_sphere[sphere]*xpos; - //printf("cs=%g, vel[0]=%g ", cs, vel[0]); + printf("cs=%g, vel[0]=%g ", cs, vel[0]); double mach_turb = 0.0; vel[0] += mach_turb*Gaussian(cs); vel[1] += mach_turb*Gaussian(cs); diff --git a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C index a8911ad7d..739bc0cec 100644 --- a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C +++ b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C @@ -58,13 +58,6 @@ int grid::MHDRK2_2ndStep(fluxes *SubgridFluxes[], this->ReturnHydroRKPointers(Prim, false); this->ReturnOldHydroRKPointers(OldPrim, false); - // ok. this line should absolutely not be needed. - // however there is one piece of hardware i've compiled on where it will not work - // without. There is absoltuely no reasons I can see how this can happen ... argh. - // this same line is in ReturnHydroRKPointers but for some reason n the kolob cluster - // in Heidlberg that is not enough. On my laptop no problem. - Prim[0] = BaryonField[0]; - #ifdef ECUDADEBUG printf("in Grid_MHDRK_2ndStep.C.\n"); for (int j=30; j < 33; j++) diff --git a/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C index 1bb206533..a213e2e04 100644 --- a/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_MHDTurbulenceInitializeGrid.C @@ -143,7 +143,7 @@ int grid::MHDTurbulenceInitializeGrid(float rho_medium, float cs_medium, float m float eint, h, dpdrho, dpde, cs; eint = cs_medium*cs_medium/(Gamma-1.0); FLOAT xc = 0.5, yc = 0.5, zc = 0.5, x, y, z, r; - FLOAT rs = 0.3; + FLOAT rs = 1.; // 0.3; n=0; for (int k = 0; k < GridDimension[2]; k++) { for (int j = 0; j < GridDimension[1]; j++) { diff --git a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C index 53a0c504c..24b6bbb97 100644 --- a/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C +++ b/src/enzo/hydro_rk/Grid_TurbulenceInitializeGrid.C @@ -32,7 +32,7 @@ void Turbulence_Generator(float **vel, int dim0, int dim1, int dim2, int ind, int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FLOAT CloudRadius, float CloudMachNumber, float CloudAngularVelocity, float InitialBField, - int SetTurbulence, int CloudType, int TurbulenceSeed, int level, + int SetTurbulence, int CloudType, int TurbulenceSeed, int PutSink, int level, int SetBaryonFields) { @@ -100,16 +100,17 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL FieldType[kdissH2INum = NumberOfBaryonFields++] = kdissH2I; } } + + if (RadiationPressure) { + FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0; + FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1; + FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2; + } + + NumberOfPhotonPackages = 0; + PhotonPackages-> NextPackage= NULL; + } - - if (RadiationPressure && RadiativeTransfer) { - FieldType[RPresNum1 = NumberOfBaryonFields++] = RadPressure0; - FieldType[RPresNum2 = NumberOfBaryonFields++] = RadPressure1; - FieldType[RPresNum3 = NumberOfBaryonFields++] = RadPressure2; - } - - NumberOfPhotonPackages = 0; - PhotonPackages-> NextPackage= NULL; #endif int idrivex, idrivey, idrivez; @@ -192,7 +193,9 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL /* assume isothermal initially */ CloudPressure = CloudDensity * pow(CloudSoundSpeed,2); - CloudInternalEnergy = pow(CloudSoundSpeed,2) / (Gamma - 1.0); + // CloudInternalEnergy = pow(CloudSoundSpeed,2) / (Gamma - 1.0); + float h, cs, dpdrho, dpde; + EOS(CloudPressure, CloudDensity, CloudInternalEnergy, h, cs, dpdrho, dpde, EOSType, 1); float InitialFractionHII = 1.2e-5; float InitialFractionHeII = 1.0e-14; @@ -228,6 +231,11 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL Velx = Vely = Velz = 0.0; + // default values: + Density = CloudDensity; + eint = CloudInternalEnergy; + + if (r < CloudRadius) { /* Type 0: uniform cloud */ @@ -291,8 +299,8 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL } else { if (CloudType == 0) { - Density = CloudDensity/10.0; - eint = CloudInternalEnergy; + Density = CloudDensity/100.0; + eint = CloudInternalEnergy*100.; } if (CloudType == 1) { @@ -551,28 +559,33 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL /* Put a sink particle if we are studying massive star formation */ - int PutSink = 0; - if (PutSink == 1 && level == MaximumRefinementLevel) { + // if (PutSink == 1 && level == MaximumRefinementLevel) { + if (PutSink == 1 && level == 0) { // set it up on level zero and make it mustrefine - double mass_p = 20.0*1.989e33; + // double mass_p = 20.0*1.989e33; + double mass_p = 0.01*1.989e33; mass_p /= MassUnits; double dx = CellWidth[0][0]; double den_p = mass_p / pow(dx,3); double t_dyn = sqrt(3*M_PI/(6.672e-8*den_p*DensityUnits)); t_dyn /= TimeUnits; - double dxm = dx / pow(2.0, MaximumRefinementLevel); + double dxm = dx / pow(RefineBy, MaximumRefinementLevel); + + printf("Adding a Sink Particle. \n"); NumberOfParticles = 1; NumberOfStars = 1; // MaximumParticleNumber = 1; + if (StellarWindFeedback) NumberOfParticleAttributes = 6; this->AllocateNewParticles(NumberOfParticles); + ParticleMass[0] = den_p; ParticleNumber[0] = 0; ParticleType[0] = PARTICLE_TYPE_MUST_REFINE; - ParticlePosition[0][0] = 0.5+0.5*dx; - ParticlePosition[1][0] = 0.5+0.5*dx; - ParticlePosition[2][0] = 0.5+0.5*dx; + ParticlePosition[0][0] = 0.5;//+0.5*dx; + ParticlePosition[1][0] = 0.5;//+0.5*dx; + ParticlePosition[2][0] = 0.5;//+0.5*dx; ParticleVelocity[0][0] = 0.0; ParticleVelocity[1][0] = 0.0; @@ -580,6 +593,18 @@ int grid::TurbulenceInitializeGrid(float CloudDensity, float CloudSoundSpeed, FL ParticleAttribute[0][0] = 0.0; // creation time ParticleAttribute[1][0] = 0; ParticleAttribute[2][0] = mass_p; + + if (StellarWindFeedback) { + ParticleAttribute[3][0] = 1.0; + ParticleAttribute[4][0] = 0.0; + ParticleAttribute[5][0] = 0.0; + } + + for (i = 0; i< MAX_DIMENSION+1; i++){ + ParticleAcceleration[i] = NULL; + } + this->ClearParticleAccelerations(); + } /* printf("XXX PutSinkParticle = %d\n", PutSinkParticle); diff --git a/src/enzo/hydro_rk/MHDTurbulenceInitialize.C b/src/enzo/hydro_rk/MHDTurbulenceInitialize.C index be92b5506..17becd5d1 100644 --- a/src/enzo/hydro_rk/MHDTurbulenceInitialize.C +++ b/src/enzo/hydro_rk/MHDTurbulenceInitialize.C @@ -61,7 +61,7 @@ int MHDTurbulenceInitialize(FILE *fptr, FILE *Outfptr, int RefineAtStart = TRUE; int RandomSeed = 1; - FLOAT rho_medium=1.0, cs=1.0, mach=1.0, Bnaught=0.0; + float rho_medium=1.0, cs=1.0, mach=1.0, Bnaught=0.0; /* read input from file */ rewind(fptr); @@ -139,6 +139,7 @@ int MHDTurbulenceInitialize(FILE *fptr, FILE *Outfptr, v_rms = sqrt(v_rms/Volume); // actuall v_rms fac = cs*mach/v_rms; + CurrentGrid = &TopGrid; while (CurrentGrid != NULL) { if (CurrentGrid->GridData->NormalizeVelocities(fac) == FAIL) { diff --git a/src/enzo/hydro_rk/TurbulenceInitialize.C b/src/enzo/hydro_rk/TurbulenceInitialize.C index 9b0eccf47..b257d432f 100644 --- a/src/enzo/hydro_rk/TurbulenceInitialize.C +++ b/src/enzo/hydro_rk/TurbulenceInitialize.C @@ -89,6 +89,7 @@ int TurbulenceInitialize(FILE *fptr, FILE *Outfptr, /* set default parameters */ int RefineAtStart = TRUE; + int PutSink = FALSE; int SetTurbulence = TRUE; int RandomSeed = 52761; float CloudDensity=1.0, CloudSoundSpeed=1.0, CloudMachNumber=1.0, CloudAngularVelocity = 0.0, InitialBField = 0.0; @@ -102,6 +103,7 @@ int TurbulenceInitialize(FILE *fptr, FILE *Outfptr, ret = 0; ret += sscanf(line, "RefineAtStart = %"ISYM, &RefineAtStart); + ret += sscanf(line, "PutSink = %"ISYM, &PutSink); ret += sscanf(line, "Density = %"FSYM, &CloudDensity); ret += sscanf(line, "SoundVelocity = %"FSYM, &CloudSoundSpeed); ret += sscanf(line, "MachNumber = %"FSYM, &CloudMachNumber); @@ -143,7 +145,7 @@ printf("Plasma beta=%g\n", CloudDensity*CloudSoundSpeed*CloudSoundSpeed/(Initial while (CurrentGrid != NULL) { if (CurrentGrid->GridData->TurbulenceInitializeGrid( CloudDensity, CloudSoundSpeed, CloudRadius, CloudMachNumber, CloudAngularVelocity, InitialBField, - SetTurbulence, CloudType, RandomSeed, 0, SetBaryonFields) == FAIL) { + SetTurbulence, CloudType, RandomSeed, PutSink, 0, SetBaryonFields) == FAIL) { fprintf(stderr, "Error in TurbulenceInitializeGrid.\n"); return FAIL; } @@ -223,7 +225,7 @@ printf("Plasma beta=%g\n", CloudDensity*CloudSoundSpeed*CloudSoundSpeed/(Initial while (Temp != NULL) { if (Temp->GridData->TurbulenceInitializeGrid( CloudDensity, CloudSoundSpeed, CloudRadius, fac, CloudAngularVelocity, InitialBField, - SetTurbulence, CloudType, RandomSeed, level+1, SetBaryonFields) == FAIL) { + SetTurbulence, CloudType, RandomSeed, PutSink, level+1, SetBaryonFields) == FAIL) { fprintf(stderr, "Error in TurbulenceInitializeGrid.\n"); return FAIL; } diff --git a/src/enzo/macros_and_parameters.h b/src/enzo/macros_and_parameters.h index 5485ea515..821ec409d 100644 --- a/src/enzo/macros_and_parameters.h +++ b/src/enzo/macros_and_parameters.h @@ -71,7 +71,7 @@ #define MAX_STATIC_REGIONS 1000 -#ifdef WINDS +#ifdef WINDS #define MAX_NUMBER_OF_PARTICLE_ATTRIBUTES 6 #else #define MAX_NUMBER_OF_PARTICLE_ATTRIBUTES 3 diff --git a/src/enzo/star_maker8.C b/src/enzo/star_maker8.C index bb6b2089c..01616e6ba 100644 --- a/src/enzo/star_maker8.C +++ b/src/enzo/star_maker8.C @@ -323,8 +323,8 @@ int star_maker8(int *nx, int *ny, int *nz, int *size, float *d, float *te, float if (mpold[bb]*pow(*dx,3)*umass < StellarWindTurnOnMass && (*t - tcpold[bb])*(*t1) < 1e5*3.1557e7) continue; int first = 0; - if (dmold[bb] > 0.99*mpold[bb]*pow(*dx,3)) first = 1; - + // if (dmold[bb] > 0.99*mpold[bb]*pow(*dx,3)) first = 1; + if (nx_jet[bb]+ny_jet[bb]+nz_jet[bb] < 0.1) first = 1; /* Decide whether the current grid contains the whole supercell */