diff --git a/example/test_problem/ELBDM/DiskHeating/DensTableExample b/example/test_problem/ELBDM/DiskHeating/DensTableExample new file mode 100644 index 0000000000..6469688293 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/DensTableExample @@ -0,0 +1,388 @@ +# radius (kpc) density (g/cm^3) +5.967904692125181e-02 1.127398940587602e-23 +6.064884386046022e-02 1.137773414631666e-23 +8.508169350868788e-02 1.127887838430367e-23 +1.049166812363333e-01 1.119757126647969e-23 +1.212970873991709e-01 1.122934549415392e-23 +1.336147053154722e-01 1.115285190358454e-23 +1.471831670433211e-01 1.107595763294151e-23 +1.701625448514275e-01 1.102849232495445e-23 +1.814957928459816e-01 1.097292804190173e-23 +1.904883873381100e-01 1.097139003467226e-23 +1.999265390215758e-01 1.090308300370590e-23 +2.098323239736347e-01 1.083716591279447e-23 +2.167073721559990e-01 1.084775735779502e-23 +2.238076777563557e-01 1.078236772746672e-23 +2.425929741607858e-01 1.073987795746785e-23 +2.465351605050992e-01 1.067006225185465e-23 +2.546127582508699e-01 1.063076380888987e-23 +2.629550143326325e-01 1.060274539679119e-23 +2.672280880696475e-01 1.055897692796490e-23 +2.759836789451054e-01 1.049901548899937e-23 +2.804684711949819e-01 1.048616160891342e-23 +2.943648772203501e-01 1.038518124168192e-23 +2.991483677696454e-01 1.037995914042226e-23 +3.089498101935024e-01 1.032469393415659e-23 +3.139703089405758e-01 1.028079000918807e-23 +3.242573846442659e-01 1.021302952519688e-23 +3.295266346636028e-01 1.016160836121982e-23 +3.403234053792958e-01 1.010744482973136e-23 +3.458537315809412e-01 1.005415294443138e-23 +3.514739267349218e-01 1.003917219720165e-23 +3.571854512304263e-01 9.989523125124693e-24 +3.629897891882715e-01 9.951347493282803e-24 +3.688884488465480e-01 9.888549539644256e-24 +3.809748891609392e-01 9.811454063637114e-24 +3.871658104387383e-01 9.759096485848333e-24 +3.934573354764133e-01 9.714160473621187e-24 +3.998510991060104e-01 9.664139243731232e-24 +4.063487627259373e-01 9.619788351952369e-24 +4.196625709594866e-01 9.502304021624033e-24 +4.264821751222990e-01 9.447247826618909e-24 +4.334125992727770e-01 9.395727571792467e-24 +4.404556442587957e-01 9.338290925946234e-24 +4.476131401923853e-01 9.257224252197791e-24 +4.548869469252716e-01 9.216126304122215e-24 +4.622789545321641e-01 9.160730707568329e-24 +4.697910838018779e-01 9.096940264673882e-24 +4.774252867364464e-01 9.061145818100212e-24 +4.851835470583421e-01 8.948121453499290e-24 +4.930678807259413e-01 8.882501814828118e-24 +5.010803364573656e-01 8.828333095847913e-24 +5.092229962628262e-01 8.746922872153511e-24 +5.174979759856403e-01 8.681156313053058e-24 +5.259074258520172e-01 8.598387743111771e-24 +5.344535310297901e-01 8.541325267723124e-24 +5.431385121962243e-01 8.459358365703292e-24 +5.519646261150534e-01 8.386026508319772e-24 +5.609341662228933e-01 8.311080785430627e-24 +5.700494632251762e-01 8.215523683164733e-24 +5.793128857017893e-01 8.104268601359583e-24 +5.887268407225342e-01 8.016729065489138e-24 +5.982937744725979e-01 7.962919799311366e-24 +6.080161728881858e-01 7.889548349968674e-24 +6.178965623024846e-01 7.776490974947911e-24 +6.279375101021236e-01 7.671181509276225e-24 +6.381416253942922e-01 7.562955912685749e-24 +6.485115596847222e-01 7.480215044038845e-24 +6.590500075666630e-01 7.380787506758842e-24 +6.697597074210659e-01 7.307986373108950e-24 +6.806434421281424e-01 7.200643136636133e-24 +6.917040397904870e-01 7.083167796225677e-24 +7.029443744679533e-01 6.975073584772102e-24 +7.143673669244570e-01 6.899719532547228e-24 +7.259759853869447e-01 6.777657693552330e-24 +7.377732463166657e-01 6.685845496820097e-24 +7.497622151929929e-01 6.584446175493451e-24 +7.619460073099772e-01 6.448278890588671e-24 +7.743277885858466e-01 6.345848922469148e-24 +7.869107763856628e-01 6.230117219011781e-24 +7.996982403573298e-01 6.116194146039691e-24 +8.126935032812208e-01 5.997420739510584e-24 +8.258999419335784e-01 5.896763714751286e-24 +8.393209879639629e-01 5.765814812935475e-24 +8.529601287869537e-01 5.653898202832119e-24 +8.668209084883440e-01 5.536559110237282e-24 +8.809069287460622e-01 5.422530729610690e-24 +8.952218497660425e-01 5.297576977199071e-24 +9.097693912333382e-01 5.179734874450507e-24 +9.245533332786562e-01 5.073740682113218e-24 +9.395775174606138e-01 4.944744337502179e-24 +9.548458477639561e-01 4.833540189994916e-24 +9.703622916139936e-01 4.703419852970734e-24 +9.861308809075273e-01 4.582097220582799e-24 +1.002155713060510e+00 4.459261558102618e-24 +1.018440952072767e+00 4.335530219059499e-24 +1.034990829609986e+00 4.221411751270444e-24 +1.051809646103302e+00 4.102861181931818e-24 +1.068901771866752e+00 3.993609954595193e-24 +1.086271648232886e+00 3.875093074669796e-24 +1.103923788706836e+00 3.738644692002888e-24 +1.121862780139121e+00 3.633716199290756e-24 +1.140093283917547e+00 3.531863816044027e-24 +1.158620037178440e+00 3.421089911806486e-24 +1.177447854037577e+00 3.304986685875536e-24 +1.196581626841120e+00 3.191067739246340e-24 +1.216026327436871e+00 3.084750746943772e-24 +1.235787008466203e+00 2.982187368755914e-24 +1.255868804676948e+00 2.882927207494185e-24 +1.276276934257672e+00 2.779230812411474e-24 +1.297016700193590e+00 2.685290917162804e-24 +1.318093491644527e+00 2.582069845779673e-24 +1.339512785345282e+00 2.491457872690128e-24 +1.361280147028730e+00 2.400033703825308e-24 +1.383401232872075e+00 2.305154887785479e-24 +1.405881790966556e+00 2.214274516356826e-24 +1.428727662811113e+00 2.129919751333224e-24 +1.451944784830252e+00 2.041837175864289e-24 +1.475539189916613e+00 1.966187812523663e-24 +1.499517008998606e+00 1.885709130000226e-24 +1.523884472633485e+00 1.806793802693778e-24 +1.548647912626373e+00 1.742408240539953e-24 +1.573813763675541e+00 1.665345908634699e-24 +1.599388565044447e+00 1.598424807127638e-24 +1.625378962260940e+00 1.528230266314316e-24 +1.651791708844081e+00 1.464449780161086e-24 +1.678633668059029e+00 1.407537256652847e-24 +1.705911814700412e+00 1.345839259242587e-24 +1.733633236904740e+00 1.293695189688161e-24 +1.761805137992213e+00 1.235558136785217e-24 +1.790434838338484e+00 1.185823273343146e-24 +1.819529777276835e+00 1.129767373659536e-24 +1.849097515031261e+00 1.089066731976520e-24 +1.879145734680979e+00 1.037051219897031e-24 +1.909682244156824e+00 9.967877771107034e-25 +1.940714978270150e+00 9.509845740895015e-25 +1.972252000774645e+00 9.132425497621187e-25 +2.004301506461671e+00 8.720033319623635e-25 +2.036871823289655e+00 8.370898666396444e-25 +2.069971414548086e+00 8.006659598407444e-25 +2.103608881056676e+00 7.669046354150865e-25 +2.137792963400222e+00 7.350448753502249e-25 +2.172532544199870e+00 7.033052029533597e-25 +2.207836650421202e+00 6.756317470616627e-25 +2.243714455719868e+00 6.460801341875871e-25 +2.280175282825335e+00 6.199175007962817e-25 +2.317228605963363e+00 5.956722695893998e-25 +2.354884053317864e+00 5.709362107016651e-25 +2.393151409532715e+00 5.484029529639351e-25 +2.432040618254322e+00 5.270122540414737e-25 +2.471561784715405e+00 5.065223910527253e-25 +2.511725178360821e+00 4.867339169001474e-25 +2.552541235516054e+00 4.687960726120802e-25 +2.594020562099020e+00 4.509965885964898e-25 +2.636173936376039e+00 4.341902532717477e-25 +2.679012311762495e+00 4.188286053185525e-25 +2.722546819669050e+00 4.052541852965764e-25 +2.766788772394110e+00 3.897588606248798e-25 +2.811749666063283e+00 3.755759078501531e-25 +2.857441183616626e+00 3.638725331586690e-25 +2.903875197844368e+00 3.526014102106276e-25 +2.951063774472086e+00 3.389943230449233e-25 +2.999019175295901e+00 3.278102370061026e-25 +3.047753861368671e+00 3.160797881004001e-25 +3.097280496237959e+00 3.055260476974870e-25 +3.147611949236613e+00 2.962831980325706e-25 +3.198761298826834e+00 2.856992425073186e-25 +3.250741835998515e+00 2.760606072965308e-25 +3.303567067722942e+00 2.663768325398848e-25 +3.357250720462485e+00 2.567382689904185e-25 +3.411806743737387e+00 2.485939157766528e-25 +3.467249313750496e+00 2.401639624166846e-25 +3.523592837070912e+00 2.322215733924149e-25 +3.580851954377493e+00 2.251484312126170e-25 +3.639041544263137e+00 2.181536483381551e-25 +3.698176727101018e+00 2.117277883791743e-25 +3.758272868973505e+00 2.058164620111951e-25 +3.819345585665008e+00 2.003846712386836e-25 +3.881410746719684e+00 1.952369655713000e-25 +3.944484479565089e+00 1.910134605778926e-25 +4.008583173702847e+00 1.864957199867478e-25 +4.073723484967353e+00 1.833473314789739e-25 +4.139922339853823e+00 1.795459555776047e-25 +4.207196939916530e+00 1.765520594911228e-25 +4.275564766238588e+00 1.736803112556888e-25 +4.345043583974348e+00 1.703727922047143e-25 +4.415651446965621e+00 1.684475265648046e-25 +4.487406702432923e+00 1.652323104280428e-25 +4.560327995742871e+00 1.625360931375033e-25 +4.634434275253227e+00 1.598834884275338e-25 +4.709744797236490e+00 1.567745285333286e-25 +4.786279130883601e+00 1.538979412730620e-25 +4.864057163388928e+00 1.509427618650410e-25 +4.943099105117886e+00 1.482265804376793e-25 +5.023425494858559e+00 1.450891895205857e-25 +5.105057205158552e+00 1.426911522859353e-25 +5.188015447748777e+00 1.406491258372321e-25 +5.272321779055169e+00 1.383627551217531e-25 +5.357998105800071e+00 1.372051320925959e-25 +5.445066690694631e+00 1.353935215453449e-25 +5.533550158223683e+00 1.349155972100849e-25 +5.623471500524676e+00 1.339710379471805e-25 +5.714854083362012e+00 1.335239466521643e-25 +5.807721652198699e+00 1.326293280740091e-25 +5.902098338366438e+00 1.319572517351957e-25 +5.998008665336096e+00 1.303465636092941e-25 +6.095477555090048e+00 1.283960558344008e-25 +6.194530334598075e+00 1.259078525215507e-25 +6.295192742398511e+00 1.227843098584635e-25 +6.397490935286224e+00 1.194880682010985e-25 +6.501451495109521e+00 1.152078608687583e-25 +6.607101435677253e+00 1.107122521717306e-25 +6.714468209778292e+00 1.064729799930043e-25 +6.823579716315042e+00 1.023803234876508e-25 +6.934464307552883e+00 9.865855641339699e-26 +7.047150796487440e+00 9.508727988824301e-26 +7.161668464331427e+00 9.250777849557254e-26 +7.278047068123453e+00 9.031832870519979e-26 +7.396316848460167e+00 8.839527300608182e-26 +7.516508537354213e+00 8.697094866095315e-26 +7.638653366219840e+00 8.555318200123143e-26 +7.762783073988281e+00 8.411829143781733e-26 +7.888929915355043e+00 8.262438830153437e-26 +8.017126669161051e+00 8.100484452740191e-26 +8.147406646910325e+00 7.883655451852640e-26 +8.279803701425720e+00 7.675162472355235e-26 +8.414352235645524e+00 7.448442229576448e-26 +8.551087211562919e+00 7.218067389704076e-26 +8.690044159310769e+00 6.994034983118804e-26 +8.831259186394028e+00 6.800102491187607e-26 +8.974768987072016e+00 6.615821287731542e-26 +9.120610851893513e+00 6.457119597011420e-26 +9.268822677386451e+00 6.310263646106842e-26 +9.419442975905227e+00 6.161875137121399e-26 +9.572510885638055e+00 6.014009117533226e-26 +9.728066180776743e+00 5.852198466617293e-26 +9.886149281852116e+00 5.681977196685150e-26 +1.004680126623700e+01 5.520544146360909e-26 +1.021006387882008e+01 5.362244869810976e-26 +1.037597954285323e+01 5.233498687827249e-26 +1.054459137097484e+01 5.125675559477709e-26 +1.071594317641283e+01 5.028055157326684e-26 +1.089007948436911e+01 4.949456198320113e-26 +1.106704554358942e+01 4.845418259324849e-26 +1.124688733812101e+01 4.734270281406557e-26 +1.142965159926151e+01 4.610657459690618e-26 +1.161538581770188e+01 4.473119573289557e-26 +1.180413825586659e+01 4.336192961242806e-26 +1.199595796045466e+01 4.198049474345829e-26 +1.219089477518417e+01 4.060763005084075e-26 +1.238899935374398e+01 3.932063770563067e-26 +1.259032317295593e+01 3.798963530067686e-26 +1.279491854615097e+01 3.672347096670550e-26 +1.300283863676258e+01 3.552642105172632e-26 +1.321413747214098e+01 3.444079556809445e-26 +1.342886995759222e+01 3.346930803304456e-26 +1.364709189064495e+01 3.255066933126238e-26 +1.386885997554932e+01 3.165388239165717e-26 +1.409423183801128e+01 3.076613453743422e-26 +1.432326604016658e+01 2.985100234729342e-26 +1.455602209579786e+01 2.889132176786951e-26 +1.479256048579899e+01 2.785412071401189e-26 +1.503294267389112e+01 2.684663511811691e-26 +1.527723112259360e+01 2.592926354137510e-26 +1.552548930945474e+01 2.516777387560006e-26 +1.577778174354622e+01 2.452648354586315e-26 +1.603417398222556e+01 2.390753120262996e-26 +1.629473264817108e+01 2.314008956877138e-26 +1.655952544669328e+01 2.222202259065510e-26 +1.682862118332825e+01 2.119959204659193e-26 +1.710208978171630e+01 2.023495539466881e-26 +1.738000230177147e+01 1.946691809676393e-26 +1.766243095814619e+01 1.889124227939758e-26 +1.794944913899604e+01 1.843617050866870e-26 +1.824113142504945e+01 1.797177329015623e-26 +1.853755360898704e+01 1.743712472510303e-26 +1.883879271513647e+01 1.688696410941136e-26 +1.914492701948668e+01 1.637985215401873e-26 +1.945603607002775e+01 1.598514612906590e-26 +1.977220070742116e+01 1.566547971448914e-26 +2.009350308600595e+01 1.530781070216547e-26 +2.042002669514638e+01 1.480413695144952e-26 +2.075185638092601e+01 1.411632022577918e-26 +2.108907836819518e+01 1.337743654751412e-26 +2.143178028297591e+01 1.274994049247744e-26 +2.178005117523133e+01 1.234751531756609e-26 +2.213398154200500e+01 1.209380067360784e-26 +2.249366335093631e+01 1.182983592506946e-26 +2.285919006415800e+01 1.145296031871243e-26 +2.323065666258166e+01 1.102486272402146e-26 +2.360815967057877e+01 1.068970036919585e-26 +2.399179718106183e+01 1.045517785201164e-26 +2.438166888097361e+01 1.025062487259429e-26 +2.477787607719042e+01 1.000221257042979e-26 +2.518052172284646e+01 9.674445038399919e-27 +2.558971044408587e+01 9.300623142695025e-27 +2.600554856724911e+01 8.906317172215793e-27 +2.642814414650212e+01 8.510128573501161e-27 +2.685760699191344e+01 8.121599116127128e-27 +2.729404869798809e+01 7.730545400769217e-27 +2.773758267266503e+01 7.367538049892562e-27 +2.818832416678584e+01 7.082507083181793e-27 +2.864639030404239e+01 6.843180604965076e-27 +2.911190011141059e+01 6.595674431495126e-27 +2.958497455007985e+01 6.340990363125453e-27 +3.006573654688394e+01 6.132170852280757e-27 +3.055431102624332e+01 5.975856075625596e-27 +3.105082494262627e+01 5.818582948021805e-27 +3.155540731353762e+01 5.608675919859258e-27 +3.206818925304366e+01 5.346513089885441e-27 +3.258930400584116e+01 5.086166534995140e-27 +3.311888698188132e+01 4.875360214555895e-27 +3.365707579155516e+01 4.718752095466635e-27 +3.420401028145119e+01 4.569138897001720e-27 +3.475983257069409e+01 4.375847206135559e-27 +3.532468708787391e+01 4.156933428844510e-27 +3.589872060857554e+01 3.971498042612144e-27 +3.648208229351721e+01 3.803602419259295e-27 +3.707492372731035e+01 3.611928422686495e-27 +3.767739895784770e+01 3.460289024596762e-27 +3.828966453633240e+01 3.386226659651756e-27 +3.891187955795716e+01 3.289777060318703e-27 +3.954420570324477e+01 3.142321188421409e-27 +4.018680728006029e+01 3.019459301092026e-27 +4.083985126630536e+01 2.904682923610275e-27 +4.150350735330775e+01 2.764423344857696e-27 +4.217794798991447e+01 2.647998649665928e-27 +4.286334842730224e+01 2.582369525606335e-27 +4.355988676451584e+01 2.490145007696827e-27 +4.426774399474667e+01 2.341964323108344e-27 +4.498710405236331e+01 2.210445046801357e-27 +4.571815386070561e+01 2.119817407825343e-27 +4.646108338065716e+01 2.032356682975057e-27 +4.721608566000528e+01 1.937820604496433e-27 +4.798335688360402e+01 1.836455570270898e-27 +4.876309642435236e+01 1.736926647034524e-27 +4.955550689499984e+01 1.665225676911459e-27 +5.036079420079591e+01 1.606346048416841e-27 +5.117916759299304e+01 1.519969600481343e-27 +5.201083972322012e+01 1.430496726598301e-27 +5.285602669873948e+01 1.374777185800932e-27 +5.371494813860085e+01 1.321910023245289e-27 +5.458782723070987e+01 1.252800768409013e-27 +5.547489078982191e+01 1.184626270664887e-27 +5.637636931647936e+01 1.119715428563933e-27 +5.729249705690637e+01 1.052066862696217e-27 +5.822351206387709e+01 9.966622765617625e-28 +5.916965625857299e+01 9.579973613082531e-28 +6.013117549344442e+01 9.238749809346667e-28 +6.110831961609595e+01 8.809036787615531e-28 +6.210134253420753e+01 8.293155893608045e-28 +6.311050228151176e+01 7.783627781564644e-28 +6.413606108484312e+01 7.383185505788139e-28 +6.517828543227679e+01 7.094419515437592e-28 +6.623744614237489e+01 6.725554750598382e-28 +6.731381843455677e+01 6.258286797972389e-28 +6.840768200061545e+01 5.799209420857523e-28 +6.951932107739358e+01 5.451912732730001e-28 +7.064902452064167e+01 5.232840544195202e-28 +7.179708588007637e+01 4.914575236267251e-28 +7.296380347565827e+01 4.596456011954125e-28 +7.414948047510977e+01 4.311178289556114e-28 +7.535442497269101e+01 4.001678545619092e-28 +7.657895006925889e+01 3.769009036940164e-28 +7.782337395362417e+01 3.561114094170975e-28 +7.908801998523197e+01 3.359943167688709e-28 +8.037321677818576e+01 3.130280832472630e-28 +8.167929828663662e+01 2.867454640850293e-28 +8.300660389156050e+01 2.660362396435601e-28 +8.435547848894376e+01 2.479844027638787e-28 +8.572627257940547e+01 2.279941602789714e-28 +8.711934235927235e+01 2.116186954887996e-28 +8.853504981313563e+01 1.990129068024858e-28 +8.997376280791148e+01 1.840421175463048e-28 +9.143585518843011e+01 1.649478050149463e-28 +9.292170687457843e+01 1.477387056955990e-28 +9.443170396001943e+01 1.354890073452201e-28 +9.596623881251975e+01 1.231771039051641e-28 +9.752571017590347e+01 1.106539849253688e-28 +9.911052327366488e+01 9.999819153420658e-29 +1.007210899142644e+02 9.052755184696057e-29 +1.023578285981359e+02 8.084078157406028e-29 +1.040211646264324e+02 7.117077350413077e-29 +1.057115302115383e+02 6.216063675393502e-29 +1.074293645893806e+02 5.353696794261152e-29 +1.091751141335608e+02 4.525721060844072e-29 +1.109492324713447e+02 3.770983442005429e-29 +1.127521806015354e+02 3.102875412146339e-29 +1.145844270142628e+02 2.489310876250868e-29 diff --git a/example/test_problem/ELBDM/DiskHeating/DispTableExample b/example/test_problem/ELBDM/DiskHeating/DispTableExample new file mode 100644 index 0000000000..567a08795b --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/DispTableExample @@ -0,0 +1,1001 @@ +# radius (kpc) 1d velocity dispersion (km/s) +9.051444959215345e-02 9.991245314344237e+00 +1.665676403020703e-01 9.974770950373200e+00 +2.168297368851021e-01 9.967924939657628e+00 +2.579357120415600e-01 9.923881567789760e+00 +2.936667891184207e-01 9.909181752434900e+00 +3.259621754414989e-01 9.880184724377930e+00 +3.555403686481572e-01 9.858400815375209e+00 +3.831672366202105e-01 9.824941485562107e+00 +4.091418862056036e-01 9.787784424040263e+00 +4.337513397389973e-01 9.776855879291269e+00 +4.572413416428418e-01 9.736616033616709e+00 +4.797112340834107e-01 9.733217868916796e+00 +5.013028308130825e-01 9.676159978358191e+00 +5.221848955793882e-01 9.654856706290374e+00 +5.423839583953626e-01 9.655754055887940e+00 +5.619105846589300e-01 9.624842557014430e+00 +5.809637408456354e-01 9.575227190106119e+00 +5.995399913132514e-01 9.589805805555189e+00 +6.176695341218641e-01 9.564170067875601e+00 +6.353619067366739e-01 9.530739777996587e+00 +6.527297342476823e-01 9.494612409016778e+00 +6.697130268715661e-01 9.513824287857258e+00 +6.863941661246031e-01 9.482170742613858e+00 +7.027689912333419e-01 9.454796418522410e+00 +7.188446924516487e-01 9.420729979463745e+00 +7.345692932481849e-01 9.412663731652367e+00 +7.500827845424061e-01 9.370348006504232e+00 +7.654373591142728e-01 9.385312719377993e+00 +7.804999278784134e-01 9.339504446262865e+00 +7.953811293006449e-01 9.345420316581190e+00 +8.100397201210403e-01 9.338389191965502e+00 +8.244635053628488e-01 9.315276259459335e+00 +8.387264751506742e-01 9.255228783714800e+00 +8.528137419016346e-01 9.266941005610802e+00 +8.667800636562711e-01 9.258116533666671e+00 +8.805700161017986e-01 9.241237610296029e+00 +8.942246502571259e-01 9.222024007435367e+00 +9.076913718999982e-01 9.200273763680448e+00 +9.210333451320235e-01 9.190156422365469e+00 +9.342267878961147e-01 9.130663313299825e+00 +9.473282036083772e-01 9.159845716734345e+00 +9.602564801501545e-01 9.139629858066456e+00 +9.730754024950397e-01 9.116926888660963e+00 +9.858303523313331e-01 9.110481821338666e+00 +9.984906671579762e-01 9.072113241060105e+00 +1.010995065164045e+00 9.086279985429703e+00 +1.023407864780478e+00 9.063227946617056e+00 +1.035707379885587e+00 9.045744479990947e+00 +1.047958031253424e+00 9.024470893824738e+00 +1.060114738259678e+00 9.021836007243493e+00 +1.072175293676792e+00 9.019826247067181e+00 +1.084076769065261e+00 8.976305088229291e+00 +1.095932037662822e+00 8.963177049471254e+00 +1.107703662875320e+00 8.951629743138472e+00 +1.119431621560145e+00 8.962875101146103e+00 +1.131054377662430e+00 8.961016741651424e+00 +1.142647108713709e+00 8.906825129744497e+00 +1.154162512767211e+00 8.900993076172583e+00 +1.165577212095682e+00 8.880261572359849e+00 +1.176947776759482e+00 8.891696317698349e+00 +1.188234271150456e+00 8.850774366001524e+00 +1.199458836258812e+00 8.859126537032051e+00 +1.210579995757309e+00 8.863669690080147e+00 +1.221671341673499e+00 8.852880925331977e+00 +1.232704023447025e+00 8.790609260295579e+00 +1.243667625393673e+00 8.783352784175660e+00 +1.254572338025936e+00 8.730794080770842e+00 +1.265408990306328e+00 8.769090154059002e+00 +1.276200350491230e+00 8.771066091949187e+00 +1.286970056494707e+00 8.749738671424270e+00 +1.297697581518465e+00 8.700815296796394e+00 +1.308389311359829e+00 8.710467026280680e+00 +1.319007322915593e+00 8.673441896453621e+00 +1.329570581306716e+00 8.655570007202479e+00 +1.340096938279985e+00 8.671804102706988e+00 +1.350546544861093e+00 8.643139166829620e+00 +1.361012985661630e+00 8.630794011317560e+00 +1.371442898126210e+00 8.622599864214271e+00 +1.381832478553920e+00 8.578731936879631e+00 +1.392171511069554e+00 8.599955504430127e+00 +1.402451600285972e+00 8.590516193851725e+00 +1.412685407845439e+00 8.530474704548379e+00 +1.422868874906741e+00 8.542868352367142e+00 +1.432994573480347e+00 8.534483786181079e+00 +1.443067858395854e+00 8.525789976068292e+00 +1.453101423971246e+00 8.489407222908676e+00 +1.463083946348457e+00 8.493135104228774e+00 +1.473062438780361e+00 8.519960808618055e+00 +1.483051250951952e+00 8.460396433425821e+00 +1.493010266912880e+00 8.479782703312120e+00 +1.502906984394867e+00 8.437813527461969e+00 +1.512733120187583e+00 8.446946665397959e+00 +1.522584185099444e+00 8.441494554820448e+00 +1.532402583601626e+00 8.416287932704822e+00 +1.542179743343738e+00 8.380436539953063e+00 +1.551947268130736e+00 8.350197626559977e+00 +1.561653970996335e+00 8.390253329492015e+00 +1.571353308265052e+00 8.361925477995605e+00 +1.581024856418998e+00 8.363211813403693e+00 +1.590631647727006e+00 8.333128505146457e+00 +1.600203587130945e+00 8.317683794089872e+00 +1.609805756475869e+00 8.306978534973876e+00 +1.619341866615469e+00 8.294721039230749e+00 +1.628853424065078e+00 8.283839910175541e+00 +1.638377919819222e+00 8.269880664799969e+00 +1.647843256089553e+00 8.274823210280525e+00 +1.657292296377826e+00 8.243372023317864e+00 +1.666705493968850e+00 8.228433790096524e+00 +1.676069831797134e+00 8.232840559516459e+00 +1.685472931023503e+00 8.229628933061854e+00 +1.694841078809259e+00 8.183429330855526e+00 +1.704154857830086e+00 8.191621012259018e+00 +1.713463839254789e+00 8.200987494747755e+00 +1.722741017433113e+00 8.180154817006043e+00 +1.732001291648819e+00 8.165203724680502e+00 +1.741230326149062e+00 8.148628362099178e+00 +1.750473518538592e+00 8.122148471995976e+00 +1.759676068848852e+00 8.130077650779143e+00 +1.768843570470882e+00 8.112870867813641e+00 +1.777976413479944e+00 8.110884681312490e+00 +1.787102863608611e+00 8.088284332144537e+00 +1.796244299958847e+00 8.052865624895301e+00 +1.805391337772920e+00 8.086877546045820e+00 +1.814474213560776e+00 8.084848525479904e+00 +1.823527596119607e+00 8.019307790688694e+00 +1.832601286160442e+00 8.043749017840216e+00 +1.841613944179411e+00 8.038914334160525e+00 +1.850656182881113e+00 8.020469326148413e+00 +1.859664238102065e+00 8.020624009049893e+00 +1.868629941357082e+00 8.014170216102398e+00 +1.877561297999507e+00 8.020700604150223e+00 +1.886483683817743e+00 7.983540975306535e+00 +1.895450404495590e+00 7.965629858874210e+00 +1.904367276352159e+00 7.940362661959291e+00 +1.913271683746916e+00 7.946254506659045e+00 +1.922161464076581e+00 7.909615065140586e+00 +1.931045087255394e+00 7.914249489894001e+00 +1.939940182133234e+00 7.895984997353528e+00 +1.948779124276933e+00 7.879651549188996e+00 +1.957624838731837e+00 7.879262761808786e+00 +1.966471882664920e+00 7.874318316117771e+00 +1.975299804771322e+00 7.839840782192348e+00 +1.984109574894776e+00 7.811255332364292e+00 +1.992886195769491e+00 7.822795912437560e+00 +2.001664282128357e+00 7.825913527695614e+00 +2.010432508647122e+00 7.814854843410613e+00 +2.019176146364370e+00 7.802430942256996e+00 +2.027905102470577e+00 7.774787728749876e+00 +2.036600009676353e+00 7.756400086907671e+00 +2.045278821184444e+00 7.759872170528842e+00 +2.053983603035677e+00 7.728151743278678e+00 +2.062672472888424e+00 7.733224695684336e+00 +2.071333113024306e+00 7.697226130631383e+00 +2.080013077752956e+00 7.689843294363090e+00 +2.088643855571413e+00 7.677536828432466e+00 +2.097254603532997e+00 7.677542114847082e+00 +2.105865109914667e+00 7.662523525325232e+00 +2.114466231849850e+00 7.652452172072873e+00 +2.123071200694739e+00 7.649438883850697e+00 +2.131632522236911e+00 7.619574138172604e+00 +2.140198305444741e+00 7.631375114126253e+00 +2.148747007139111e+00 7.607122685381352e+00 +2.157317601430382e+00 7.627721580706606e+00 +2.165852253169238e+00 7.561091729032874e+00 +2.174393632598873e+00 7.573994258626739e+00 +2.182916933550291e+00 7.560845881073616e+00 +2.191471504556234e+00 7.540601003899781e+00 +2.199970953701423e+00 7.510339994827286e+00 +2.208489688190068e+00 7.528279179547338e+00 +2.216989794986489e+00 7.531784855724363e+00 +2.225471856967510e+00 7.507171096923062e+00 +2.233954023274007e+00 7.476527111821166e+00 +2.242425446057462e+00 7.483195090702154e+00 +2.250878878374888e+00 7.472079789586532e+00 +2.259351283365691e+00 7.479042076041706e+00 +2.267801783086895e+00 7.447910674864783e+00 +2.276245721448047e+00 7.446598653057381e+00 +2.284682700832589e+00 7.425308530250501e+00 +2.293140148488713e+00 7.410482614875118e+00 +2.301584866965947e+00 7.387314613386135e+00 +2.310012258708539e+00 7.406262224653140e+00 +2.318428018997302e+00 7.382405565936730e+00 +2.326806961961323e+00 7.357158082804374e+00 +2.335183570942871e+00 7.369526310796542e+00 +2.343541751756464e+00 7.366978240664894e+00 +2.351906543804045e+00 7.348166849515945e+00 +2.360277228177180e+00 7.342588841032855e+00 +2.368658665785069e+00 7.324841086735391e+00 +2.377034128780185e+00 7.301419446033844e+00 +2.385380326443058e+00 7.297098680817219e+00 +2.393731403107788e+00 7.294085350492770e+00 +2.402078575366105e+00 7.291463834021896e+00 +2.410403279123224e+00 7.265827183192958e+00 +2.418712973813680e+00 7.247520145213514e+00 +2.427050814477540e+00 7.229690523742098e+00 +2.435378135474957e+00 7.236392741531482e+00 +2.443699285990923e+00 7.236892051686334e+00 +2.452033533018096e+00 7.227921983499942e+00 +2.460338270139747e+00 7.218505323924853e+00 +2.468651119860710e+00 7.198535601990642e+00 +2.476921153293524e+00 7.204656826773620e+00 +2.485213327702420e+00 7.165087726812131e+00 +2.493496654654376e+00 7.175685110423185e+00 +2.501786862504055e+00 7.184072849065688e+00 +2.510095434434051e+00 7.131211586668751e+00 +2.518413457556911e+00 7.136859972984735e+00 +2.526697456265067e+00 7.127468262539246e+00 +2.534945999203352e+00 7.127993082080171e+00 +2.543233098790593e+00 7.132116753954445e+00 +2.551503872053626e+00 7.096895676859104e+00 +2.559722589115218e+00 7.111635476162524e+00 +2.567995913792863e+00 7.091447268850341e+00 +2.576232190511472e+00 7.078394640459382e+00 +2.584481221632604e+00 7.074994009985541e+00 +2.592722378676825e+00 7.048491613550237e+00 +2.600959225457898e+00 7.038103553376057e+00 +2.609183994065232e+00 7.022390981493949e+00 +2.617411311385273e+00 7.033429226364427e+00 +2.625644235565350e+00 7.013837938356181e+00 +2.633882948163897e+00 6.989914717125933e+00 +2.642091942666962e+00 7.009461809926071e+00 +2.650305512980348e+00 6.982962186294321e+00 +2.658533246365642e+00 6.995243988888186e+00 +2.666746872682312e+00 6.990989866869168e+00 +2.674937392302140e+00 6.944173506889364e+00 +2.683150410954906e+00 6.963202361744078e+00 +2.691374207483112e+00 6.939473960572973e+00 +2.699572428959399e+00 6.933174500289426e+00 +2.707780726096840e+00 6.920847031086003e+00 +2.715986511752155e+00 6.925440543511925e+00 +2.724182211675575e+00 6.904073952700576e+00 +2.732344916724764e+00 6.890884558416359e+00 +2.740514960907789e+00 6.898957714479621e+00 +2.748696473940421e+00 6.902599060738130e+00 +2.756859586069139e+00 6.897139252191656e+00 +2.765018565014290e+00 6.854141195835196e+00 +2.773208021942693e+00 6.862680183697208e+00 +2.781392630925537e+00 6.843945952547908e+00 +2.789574826524269e+00 6.834240131786919e+00 +2.797720627166225e+00 6.801597620706762e+00 +2.805859120914128e+00 6.803250015292916e+00 +2.814014400809385e+00 6.797816865685114e+00 +2.822155126358587e+00 6.794610479098450e+00 +2.830314006296903e+00 6.786892277059615e+00 +2.838477934777611e+00 6.777580846030349e+00 +2.846660996334601e+00 6.778627037607439e+00 +2.854842516658457e+00 6.756274868911380e+00 +2.863008476912600e+00 6.745973570706161e+00 +2.871144783573757e+00 6.757733574373458e+00 +2.879289653251275e+00 6.738909149250871e+00 +2.887434640733716e+00 6.708406020893334e+00 +2.895569873444713e+00 6.710000923689988e+00 +2.903711340261766e+00 6.693768291037554e+00 +2.911908295134767e+00 6.685359214371700e+00 +2.920096338354381e+00 6.670608041030361e+00 +2.928232696022754e+00 6.659105345590580e+00 +2.936390842255041e+00 6.657812417668506e+00 +2.944560619079065e+00 6.622170276789201e+00 +2.952737389364393e+00 6.625827356530229e+00 +2.960880279570787e+00 6.610335246203304e+00 +2.969012770921461e+00 6.648910668810753e+00 +2.977141887360566e+00 6.588772640075120e+00 +2.985318736643642e+00 6.619080441608707e+00 +2.993511349322409e+00 6.593131961845067e+00 +3.001681937239850e+00 6.587913801261859e+00 +3.009842053367899e+00 6.571872010545192e+00 +3.017999446038180e+00 6.558342605615748e+00 +3.026177065799167e+00 6.554540082607050e+00 +3.034333209599952e+00 6.542646907603297e+00 +3.042482281682498e+00 6.528952350403667e+00 +3.050619066799459e+00 6.528435495725641e+00 +3.058759146148514e+00 6.512735368242679e+00 +3.066921831120243e+00 6.528781203226288e+00 +3.075107875067048e+00 6.499774043198638e+00 +3.083257604824291e+00 6.498123461253160e+00 +3.091378510907192e+00 6.486276410201403e+00 +3.099523900100231e+00 6.457603613096756e+00 +3.107707101900611e+00 6.472581387932846e+00 +3.115899853474430e+00 6.461509149905489e+00 +3.124060357544717e+00 6.427010358551749e+00 +3.132209481698509e+00 6.427060195142793e+00 +3.140357233432207e+00 6.455309960847099e+00 +3.148516452580783e+00 6.444426858558794e+00 +3.156669904365537e+00 6.389151029718747e+00 +3.164822414822372e+00 6.390982711226430e+00 +3.172995552918557e+00 6.387622899298519e+00 +3.181174994506634e+00 6.393465824144870e+00 +3.189354294700637e+00 6.375697101673302e+00 +3.197559576582689e+00 6.380841216242212e+00 +3.205712208474667e+00 6.370941169988106e+00 +3.213880567635618e+00 6.356783714053370e+00 +3.222051236624577e+00 6.359750176127593e+00 +3.230220468960430e+00 6.351763817138030e+00 +3.238419811324323e+00 6.296777232335423e+00 +3.246625583846536e+00 6.313101118393933e+00 +3.254806943901181e+00 6.336160215897265e+00 +3.262995596430245e+00 6.306112539813000e+00 +3.271160142736533e+00 6.297539064532894e+00 +3.279320979699928e+00 6.302851471295559e+00 +3.287504155256635e+00 6.287098995688323e+00 +3.295708805811539e+00 6.277885621240162e+00 +3.303914571221958e+00 6.260010269707073e+00 +3.312118256346885e+00 6.262133007535448e+00 +3.320296050344696e+00 6.247878429715598e+00 +3.328498363341823e+00 6.217171244236360e+00 +3.336699122031593e+00 6.242690238892128e+00 +3.344905476369670e+00 6.207488487913248e+00 +3.353109306902884e+00 6.207974684342064e+00 +3.361345255511544e+00 6.209998965011036e+00 +3.369572602638848e+00 6.189154619168844e+00 +3.377790205587923e+00 6.167400733733901e+00 +3.386003590971256e+00 6.171585031025693e+00 +3.394234377791271e+00 6.169684996729550e+00 +3.402459111003897e+00 6.147795946782023e+00 +3.410669232506410e+00 6.128594575851698e+00 +3.418926333994528e+00 6.145180138474571e+00 +3.427151507857932e+00 6.132999335138090e+00 +3.435360871355089e+00 6.114339987046916e+00 +3.443583739037654e+00 6.100391686175635e+00 +3.451784387280049e+00 6.094120131511723e+00 +3.460012762348004e+00 6.092642333315316e+00 +3.468234548719138e+00 6.076114616376823e+00 +3.476492000186071e+00 6.075112136796298e+00 +3.484731291659109e+00 6.066619898788500e+00 +3.492981933805859e+00 6.074125890057736e+00 +3.501232557361044e+00 6.031583223664412e+00 +3.509489740519445e+00 6.031029557528474e+00 +3.517758349272442e+00 6.047449265396898e+00 +3.526035905366625e+00 6.028682189613619e+00 +3.534314018553630e+00 6.023039962449616e+00 +3.542595956391936e+00 6.008333359930960e+00 +3.550867313870244e+00 6.021585002543975e+00 +3.559143469064073e+00 6.007139654187062e+00 +3.567391087349857e+00 5.956507262819867e+00 +3.575669373035081e+00 5.988930404341501e+00 +3.583953954027788e+00 5.963479067823802e+00 +3.592253187278612e+00 5.945324192139310e+00 +3.600547745042975e+00 5.942947617758270e+00 +3.608851379331711e+00 5.955495988472750e+00 +3.617175831979941e+00 5.940192903720446e+00 +3.625472304214326e+00 5.923788147442503e+00 +3.633764190621232e+00 5.920871312836272e+00 +3.642075357990406e+00 5.924743163068129e+00 +3.650401110722733e+00 5.914470534137101e+00 +3.658711323092430e+00 5.890086267018806e+00 +3.667028096349673e+00 5.895898913187962e+00 +3.675356844717340e+00 5.877799849120480e+00 +3.683710607464128e+00 5.853363310397377e+00 +3.692038106272393e+00 5.851106764640295e+00 +3.700361496907619e+00 5.866414118282698e+00 +3.708716332126331e+00 5.836065813789561e+00 +3.717080070104264e+00 5.836475213874331e+00 +3.725450916046979e+00 5.813329484521305e+00 +3.733842636404055e+00 5.777121785009736e+00 +3.742256340638015e+00 5.820344345668640e+00 +3.750639806736132e+00 5.791101270082161e+00 +3.759001352793136e+00 5.788900958879745e+00 +3.767388508861550e+00 5.791697541766860e+00 +3.775814239296701e+00 5.763180883667969e+00 +3.784224309892063e+00 5.750513821549282e+00 +3.792647280128901e+00 5.749937047936728e+00 +3.801061121887466e+00 5.745514677783294e+00 +3.809457463682138e+00 5.735052295150908e+00 +3.817866443523861e+00 5.728578899894710e+00 +3.826277988129487e+00 5.731474114310211e+00 +3.834718335655524e+00 5.706289245172696e+00 +3.843177875969209e+00 5.699642707864293e+00 +3.851600040321013e+00 5.678811710958494e+00 +3.860068217276394e+00 5.684268423210794e+00 +3.868521172276363e+00 5.672558636582172e+00 +3.876996420865674e+00 5.672401892949099e+00 +3.885418371418381e+00 5.654626698147651e+00 +3.893880144785535e+00 5.647901015550290e+00 +3.902362581633269e+00 5.641189728573530e+00 +3.910838394318595e+00 5.630696047562357e+00 +3.919346170513434e+00 5.612575254431053e+00 +3.927819217722589e+00 5.614732545631213e+00 +3.936294445076828e+00 5.608937364279955e+00 +3.944775422841003e+00 5.608364906797574e+00 +3.953257875639325e+00 5.606081811410886e+00 +3.961732938599100e+00 5.567204020184716e+00 +3.970244496213813e+00 5.572686037058609e+00 +3.978752949715724e+00 5.563607968528980e+00 +3.987273543226327e+00 5.567584958690826e+00 +3.995822413309198e+00 5.539127048651292e+00 +4.004358284083194e+00 5.545317805230821e+00 +4.012900565739628e+00 5.535605338375364e+00 +4.021464713842511e+00 5.540213257741856e+00 +4.030039466810087e+00 5.530602806265496e+00 +4.038632176028856e+00 5.503882373185770e+00 +4.047236896494178e+00 5.492071343806506e+00 +4.055845886992667e+00 5.485695833382739e+00 +4.064427700379603e+00 5.488329570869939e+00 +4.073005666551637e+00 5.485833557219857e+00 +4.081629959946633e+00 5.474381712548883e+00 +4.090232222949410e+00 5.463701392727389e+00 +4.098820779152379e+00 5.461466928218084e+00 +4.107399820926972e+00 5.449451295731829e+00 +4.116009600562549e+00 5.450277950209292e+00 +4.124617164996230e+00 5.425281153768283e+00 +4.133263358140320e+00 5.424148701580206e+00 +4.141907316933409e+00 5.442616402845466e+00 +4.150568036060033e+00 5.421989005677480e+00 +4.159221468867925e+00 5.378517745418850e+00 +4.167897469550966e+00 5.386390808069870e+00 +4.176584654274534e+00 5.386331193037437e+00 +4.185242164115321e+00 5.371954561681155e+00 +4.193931551261469e+00 5.375467405498643e+00 +4.202622533251951e+00 5.368480096826873e+00 +4.211311149948073e+00 5.365167805346929e+00 +4.220002043483638e+00 5.336262284943715e+00 +4.228686334878920e+00 5.327957555607427e+00 +4.237411288377338e+00 5.329814095034768e+00 +4.246151832155660e+00 5.294817188182140e+00 +4.254901090794918e+00 5.298202313842639e+00 +4.263664224394870e+00 5.309325801766002e+00 +4.272434295977463e+00 5.310894187401803e+00 +4.281179133596256e+00 5.308177026009643e+00 +4.289933101117387e+00 5.274037131533065e+00 +4.298707920825039e+00 5.274001441944628e+00 +4.307495079708110e+00 5.275252338826882e+00 +4.316298002165444e+00 5.255275573282980e+00 +4.325090768429716e+00 5.255701015002868e+00 +4.333882666801149e+00 5.232701959074258e+00 +4.342688585350415e+00 5.229627949067891e+00 +4.351490178001815e+00 5.226271350114633e+00 +4.360301197159095e+00 5.216522580279598e+00 +4.369130701508078e+00 5.217157113714349e+00 +4.378000668845476e+00 5.208127274005363e+00 +4.386855924316897e+00 5.196374143321233e+00 +4.395715769297293e+00 5.184254548258965e+00 +4.404546478259401e+00 5.192365382912064e+00 +4.413395412956774e+00 5.150760121280335e+00 +4.422284592755359e+00 5.176718870524613e+00 +4.431173737431078e+00 5.166754399164956e+00 +4.440072004637432e+00 5.155916849897993e+00 +4.449010242127653e+00 5.158934281625097e+00 +4.457929922439059e+00 5.137117517801758e+00 +4.466831002323119e+00 5.135967613376123e+00 +4.475752437663731e+00 5.121977546364575e+00 +4.484706731190525e+00 5.115655149734904e+00 +4.493668179762668e+00 5.118254676788109e+00 +4.502640914966831e+00 5.091449052202423e+00 +4.511608688247889e+00 5.109578022074562e+00 +4.520571438145049e+00 5.090007416012543e+00 +4.529558720313506e+00 5.077083683279709e+00 +4.538552734241494e+00 5.076086250966188e+00 +4.547579762521288e+00 5.059528707445230e+00 +4.556598215724882e+00 5.049910738655414e+00 +4.565635392304602e+00 5.038614272510560e+00 +4.574675171292834e+00 5.040159169903469e+00 +4.583736083399470e+00 5.040270092617535e+00 +4.592804581618599e+00 5.012833782133060e+00 +4.601872232730690e+00 5.031726480623454e+00 +4.610959948010735e+00 5.008479608492928e+00 +4.620047801278620e+00 5.000083574479217e+00 +4.629125921466623e+00 5.004293068014920e+00 +4.638226400852952e+00 4.979163958439991e+00 +4.647335080506092e+00 4.976879921518334e+00 +4.656425421145590e+00 4.965170701962748e+00 +4.665559981544257e+00 4.971501315409909e+00 +4.674695827545593e+00 4.981398293091504e+00 +4.683841261396431e+00 4.944994852649408e+00 +4.692999054161809e+00 4.926018895827054e+00 +4.702182020274011e+00 4.934480214505056e+00 +4.711398158763521e+00 4.919917304298099e+00 +4.720597919972222e+00 4.917677535261827e+00 +4.729811853271841e+00 4.905730629211804e+00 +4.739042472138693e+00 4.913822311556495e+00 +4.748259040691588e+00 4.896776900406905e+00 +4.757492914311331e+00 4.878581277731388e+00 +4.766725352930628e+00 4.881030016454018e+00 +4.775972896585908e+00 4.852164547430513e+00 +4.785232560867947e+00 4.874912399621820e+00 +4.794506320662743e+00 4.851335779694763e+00 +4.803770100913604e+00 4.842258636417895e+00 +4.813079943463439e+00 4.829743536564853e+00 +4.822375856599675e+00 4.836507327414035e+00 +4.831702863948692e+00 4.805072850888902e+00 +4.841062136910322e+00 4.796855462690245e+00 +4.850424349465538e+00 4.796165860042131e+00 +4.859761405165313e+00 4.795787491703195e+00 +4.869121276832495e+00 4.774669667803403e+00 +4.878510911924899e+00 4.781814022099654e+00 +4.887916273838806e+00 4.764720554962618e+00 +4.897349996616254e+00 4.773298881391852e+00 +4.906784472846828e+00 4.755226407709046e+00 +4.916217732070275e+00 4.750157592679571e+00 +4.925682410215659e+00 4.742510720661729e+00 +4.935124794979468e+00 4.739741448143403e+00 +4.944616624121466e+00 4.740618582219865e+00 +4.954104676801201e+00 4.711137374999883e+00 +4.963566843768075e+00 4.705749503263436e+00 +4.973071339956582e+00 4.686342159934660e+00 +4.982572864662687e+00 4.703775946351056e+00 +4.992086260905965e+00 4.687943329808800e+00 +5.001604908881880e+00 4.682017970855279e+00 +5.011114119275693e+00 4.682988517847212e+00 +5.020677004353428e+00 4.675026837513317e+00 +5.030213056057235e+00 4.668576588278317e+00 +5.039795585267298e+00 4.660036298441026e+00 +5.049386756716143e+00 4.653250345678063e+00 +5.059004513087085e+00 4.648107212457630e+00 +5.068649715405418e+00 4.640870185160491e+00 +5.078305337455581e+00 4.630313385276423e+00 +5.087951265598921e+00 4.618887242233978e+00 +5.097579717210087e+00 4.617276139345906e+00 +5.107232409848767e+00 4.600385082228846e+00 +5.116898276480450e+00 4.619817069104812e+00 +5.126598840406791e+00 4.596163338636813e+00 +5.136289281213891e+00 4.580727261943482e+00 +5.145997827335001e+00 4.587964108400219e+00 +5.155714137546062e+00 4.594043145500852e+00 +5.165451438050630e+00 4.580772572196066e+00 +5.175173215542804e+00 4.561510705457366e+00 +5.184895999236295e+00 4.556490367751648e+00 +5.194650629585702e+00 4.551853217111221e+00 +5.204411237391924e+00 4.558323917451319e+00 +5.214177860834836e+00 4.530017738672282e+00 +5.223998907254732e+00 4.536692236889771e+00 +5.233824512836004e+00 4.518127143401045e+00 +5.243674688703988e+00 4.522233641566471e+00 +5.253523520673234e+00 4.524616024354042e+00 +5.263374300562907e+00 4.493678738773324e+00 +5.273272747419278e+00 4.484919161838895e+00 +5.283210217370249e+00 4.483897144598834e+00 +5.293175432905500e+00 4.487273213793420e+00 +5.303137254568096e+00 4.462619884504717e+00 +5.313088024258258e+00 4.466933937369546e+00 +5.323049394255925e+00 4.457223812254113e+00 +5.333032191086468e+00 4.457018686207477e+00 +5.343026445074742e+00 4.446116646122910e+00 +5.353054264914664e+00 4.446513114394743e+00 +5.363094257083325e+00 4.428586745988746e+00 +5.373112894536689e+00 4.431039527821000e+00 +5.383132838183279e+00 4.428500608686486e+00 +5.393184753807840e+00 4.416071571481235e+00 +5.403277217087020e+00 4.386892264325899e+00 +5.413408681065474e+00 4.390408399078179e+00 +5.423521932024232e+00 4.393324491515308e+00 +5.433601762420428e+00 4.377772785394577e+00 +5.443721244311186e+00 4.362181425635068e+00 +5.453873215964512e+00 4.360691133339622e+00 +5.464046783147966e+00 4.364308010844866e+00 +5.474234603343787e+00 4.351496394897069e+00 +5.484399077201892e+00 4.350923744384971e+00 +5.494592280467954e+00 4.339980787976829e+00 +5.504831092698558e+00 4.321800879707885e+00 +5.515095595742240e+00 4.315707676033429e+00 +5.525362159830939e+00 4.315028268156754e+00 +5.535623236182191e+00 4.311569153735311e+00 +5.545929300969305e+00 4.307071394712189e+00 +5.556217629370403e+00 4.302366664119734e+00 +5.566529093853396e+00 4.287166739060715e+00 +5.576892170070350e+00 4.270887937084442e+00 +5.587263918344445e+00 4.282682830077253e+00 +5.597643046175391e+00 4.238171031655002e+00 +5.608027632433751e+00 4.264841539358067e+00 +5.618455629792979e+00 4.245642393659355e+00 +5.628844273246369e+00 4.241949425233143e+00 +5.639271311521084e+00 4.232622628955517e+00 +5.649744738418186e+00 4.228054675311080e+00 +5.660187736499131e+00 4.212700078049735e+00 +5.670670141784420e+00 4.208740522485945e+00 +5.681198502276006e+00 4.201062745367489e+00 +5.691775408827341e+00 4.183294694361227e+00 +5.702344654388274e+00 4.203908358668071e+00 +5.712923227070999e+00 4.181191314016644e+00 +5.723498053079581e+00 4.174671475274480e+00 +5.734105347429111e+00 4.166529420298161e+00 +5.744713826034659e+00 4.145359755128664e+00 +5.755347752168925e+00 4.140083160979465e+00 +5.766020358978351e+00 4.135408826683570e+00 +5.776679189107777e+00 4.118826813554116e+00 +5.787343653947743e+00 4.111715760253884e+00 +5.798086506056043e+00 4.114878616185942e+00 +5.808826618129560e+00 4.104633352353054e+00 +5.819558092283962e+00 4.097525332143333e+00 +5.830305252217725e+00 4.094188224219189e+00 +5.841058097431667e+00 4.087743995939009e+00 +5.851875923687886e+00 4.077281304476455e+00 +5.862714314089049e+00 4.069570006436522e+00 +5.873550734935986e+00 4.064410321185776e+00 +5.884409261865312e+00 4.070857378791387e+00 +5.895321586829691e+00 4.048405293873888e+00 +5.906217192346527e+00 4.034720400370873e+00 +5.917132628894704e+00 4.038426614042643e+00 +5.928074406672715e+00 4.028283562531364e+00 +5.939038293813677e+00 4.017970787539179e+00 +5.950057961466672e+00 4.005217508948509e+00 +5.961034005994071e+00 4.011435747635095e+00 +5.972061882761555e+00 3.993330964943866e+00 +5.983116706955193e+00 3.982999950854193e+00 +5.994181016035154e+00 3.988102212664813e+00 +6.005249430749426e+00 3.971027047731904e+00 +6.016367356999431e+00 3.969484295395704e+00 +6.027518880963474e+00 3.976506681134188e+00 +6.038661341953857e+00 3.959692703324433e+00 +6.049839686373002e+00 3.946899767555518e+00 +6.061064807161273e+00 3.934828326863590e+00 +6.072266551293311e+00 3.926213968901975e+00 +6.083523457123882e+00 3.923250943499193e+00 +6.094780048272005e+00 3.934329662296955e+00 +6.106044436908094e+00 3.913287550885834e+00 +6.117294302736195e+00 3.894921658596802e+00 +6.128627491427410e+00 3.897829861847746e+00 +6.140015424554287e+00 3.892053482710369e+00 +6.151406629531670e+00 3.887467831796218e+00 +6.162773073531773e+00 3.881531122604892e+00 +6.174158206700111e+00 3.849038971658570e+00 +6.185537937527248e+00 3.857634632155060e+00 +6.196963490298209e+00 3.848824044640410e+00 +6.208436260258430e+00 3.851828455597820e+00 +6.219965890710361e+00 3.836317275893189e+00 +6.231475437997032e+00 3.828929443504019e+00 +6.243017897834330e+00 3.821974450331766e+00 +6.254549192562205e+00 3.811442548545233e+00 +6.266110505690510e+00 3.802176043632312e+00 +6.277724953484727e+00 3.796961261388644e+00 +6.289360780343508e+00 3.785181301025694e+00 +6.300985941326545e+00 3.779938655739200e+00 +6.312596206790520e+00 3.789595549161791e+00 +6.324316135092385e+00 3.764974479900200e+00 +6.336086912594243e+00 3.760405167234205e+00 +6.347848073390830e+00 3.764123707896650e+00 +6.359592680163706e+00 3.759473869822676e+00 +6.371330616837175e+00 3.748180643270865e+00 +6.383153394093961e+00 3.725686483824433e+00 +6.394994658128362e+00 3.725454766999768e+00 +6.406828334162558e+00 3.719217942042828e+00 +6.418714048111982e+00 3.731801971669117e+00 +6.430598501025155e+00 3.706170743551907e+00 +6.442606887718113e+00 3.712373043236293e+00 +6.454631431745382e+00 3.684589614619088e+00 +6.466676304657094e+00 3.692669212766601e+00 +6.478687234267866e+00 3.683940282651618e+00 +6.490746893859312e+00 3.677124527539288e+00 +6.502831437156746e+00 3.675902130854609e+00 +6.514940906336868e+00 3.662807395603968e+00 +6.527040044835199e+00 3.651064962347081e+00 +6.539210499588886e+00 3.635077915945885e+00 +6.551412362792649e+00 3.624785641303959e+00 +6.563663907662078e+00 3.632143157997755e+00 +6.575913317597744e+00 3.619810662283532e+00 +6.588198286490552e+00 3.626786258060599e+00 +6.600550964003801e+00 3.606271748617919e+00 +6.612879814222400e+00 3.612541804473595e+00 +6.625249535557891e+00 3.597832857883881e+00 +6.637627551542804e+00 3.589990864400884e+00 +6.650041571608855e+00 3.586210631997906e+00 +6.662453647102328e+00 3.580031935638931e+00 +6.674944308261752e+00 3.564159599850784e+00 +6.687458821008684e+00 3.557481138867959e+00 +6.699973497007099e+00 3.556809946468211e+00 +6.712535782128611e+00 3.545766443237728e+00 +6.725098278510889e+00 3.541567805042965e+00 +6.737688462452954e+00 3.529624469542065e+00 +6.750333582499698e+00 3.530259642717706e+00 +6.763027242205489e+00 3.515493466038007e+00 +6.775700601414991e+00 3.515014580982985e+00 +6.788457909891222e+00 3.490243347163166e+00 +6.801227729565912e+00 3.483401110470984e+00 +6.814020648676300e+00 3.496892895717156e+00 +6.826829916581953e+00 3.479763410352272e+00 +6.839697749896733e+00 3.463988147804116e+00 +6.852590419132625e+00 3.464981306958875e+00 +6.865516880795735e+00 3.457694148081713e+00 +6.878479405327546e+00 3.453219137074385e+00 +6.891434287204564e+00 3.446343029382040e+00 +6.904445503159838e+00 3.445000005330413e+00 +6.917491693237846e+00 3.437630234298092e+00 +6.930602556801307e+00 3.416096272189539e+00 +6.943730462113793e+00 3.417968696465262e+00 +6.956855071770337e+00 3.407970942619714e+00 +6.970015777978787e+00 3.401138733451607e+00 +6.983223701637078e+00 3.390325529293037e+00 +6.996467242132024e+00 3.388455441762548e+00 +7.009757435628131e+00 3.378534890910256e+00 +7.023026002015287e+00 3.389052971449339e+00 +7.036347433039875e+00 3.376967902360148e+00 +7.049693087589527e+00 3.374563039810714e+00 +7.063103789147761e+00 3.359225768396887e+00 +7.076510691114111e+00 3.357062514749456e+00 +7.089983446381979e+00 3.341264204117231e+00 +7.103516685581478e+00 3.344801651327169e+00 +7.117070813556541e+00 3.320378136147368e+00 +7.130675556166977e+00 3.325347992679506e+00 +7.144348708940647e+00 3.315067091648295e+00 +7.157997911077988e+00 3.319095004357175e+00 +7.171713846345034e+00 3.298692262271332e+00 +7.185520555927099e+00 3.294004795751782e+00 +7.199352789276878e+00 3.296492584580011e+00 +7.213150352464662e+00 3.288689178859935e+00 +7.226933277116557e+00 3.274709962256319e+00 +7.240795032642382e+00 3.274648148013556e+00 +7.254712642923061e+00 3.263714738372886e+00 +7.268649384195200e+00 3.262457121231561e+00 +7.282639306327520e+00 3.250119470400902e+00 +7.296662972046594e+00 3.243701913355220e+00 +7.310719578453972e+00 3.241622957692827e+00 +7.324825346708137e+00 3.223906453135259e+00 +7.338959120098628e+00 3.222984262777465e+00 +7.353114157845634e+00 3.212457944089323e+00 +7.367309116070627e+00 3.218138424160495e+00 +7.381585742575831e+00 3.193270225198602e+00 +7.395837986469177e+00 3.197618970958155e+00 +7.410163118954403e+00 3.180786987330056e+00 +7.424555271853719e+00 3.185470799304978e+00 +7.438928914582668e+00 3.171454996451875e+00 +7.453367745513813e+00 3.158677442062837e+00 +7.467885085606354e+00 3.159600499014134e+00 +7.482456713182288e+00 3.148405601844239e+00 +7.497003259886452e+00 3.142541647137014e+00 +7.511680694519247e+00 3.138323682810002e+00 +7.526411491834713e+00 3.122884371444718e+00 +7.541082882746271e+00 3.121657257768797e+00 +7.555767066141364e+00 3.107898202274828e+00 +7.570566068647715e+00 3.094396285418785e+00 +7.585483002893878e+00 3.096515255015056e+00 +7.600404957188168e+00 3.090070121791309e+00 +7.615326452701997e+00 3.089287073196698e+00 +7.630305924401187e+00 3.081043190651930e+00 +7.645370879666919e+00 3.066527009575232e+00 +7.660417643131973e+00 3.052964649754162e+00 +7.675570274651377e+00 3.041446693449690e+00 +7.690736909180981e+00 3.038813503129471e+00 +7.705952836691163e+00 3.027987904071338e+00 +7.721203389870149e+00 3.022885647193073e+00 +7.736509213621817e+00 3.027800554221720e+00 +7.751867804722576e+00 3.018192923810358e+00 +7.767314319866564e+00 3.006624326937754e+00 +7.782783613108003e+00 2.999578062018823e+00 +7.798287410962590e+00 2.988309102608365e+00 +7.813836723253146e+00 2.981355359294452e+00 +7.829451971594067e+00 2.969953907124979e+00 +7.845154653468577e+00 2.976027125229647e+00 +7.860846523379060e+00 2.955437493873336e+00 +7.876609381993871e+00 2.959147276237673e+00 +7.892407076995961e+00 2.952527453760685e+00 +7.908282093620814e+00 2.934904315676216e+00 +7.924199672240439e+00 2.935671588057940e+00 +7.940117547607845e+00 2.926752115973458e+00 +7.956061650389660e+00 2.913819908714446e+00 +7.972112247388439e+00 2.915559175516016e+00 +7.988230053090938e+00 2.910848702066093e+00 +8.004387327964469e+00 2.900498667936511e+00 +8.020581442339507e+00 2.895886947158500e+00 +8.036822189012153e+00 2.889661459104354e+00 +8.053167406874495e+00 2.877645447714864e+00 +8.069594638246542e+00 2.851291559553903e+00 +8.086065562105691e+00 2.854703612437219e+00 +8.102581251064896e+00 2.848455824487961e+00 +8.119141568332681e+00 2.836860609557971e+00 +8.135740858003214e+00 2.831850484129236e+00 +8.152424794636607e+00 2.821343540160209e+00 +8.169206669873219e+00 2.809806878881043e+00 +8.185989338112380e+00 2.808367922999205e+00 +8.202878753866104e+00 2.797027537916351e+00 +8.219845899358495e+00 2.794984809305661e+00 +8.236871918505098e+00 2.783748364906226e+00 +8.253884664213098e+00 2.788465926541341e+00 +8.270966343595335e+00 2.770073158525793e+00 +8.288078930730448e+00 2.770535737575271e+00 +8.305279777981903e+00 2.758736181909892e+00 +8.322567842336857e+00 2.752487372629047e+00 +8.339967103584096e+00 2.741627804880288e+00 +8.357398818217588e+00 2.744461348240287e+00 +8.374838098525723e+00 2.737274118538594e+00 +8.392361808598423e+00 2.725478131928304e+00 +8.409955759633485e+00 2.713297302051984e+00 +8.427668201515749e+00 2.706080654565976e+00 +8.445405487785710e+00 2.695233698666157e+00 +8.463226811505059e+00 2.704061381726464e+00 +8.481141308562794e+00 2.674671510643962e+00 +8.499120405968828e+00 2.679813129031822e+00 +8.517107749700537e+00 2.673685767943390e+00 +8.535156814913350e+00 2.657247221786385e+00 +8.553332080364903e+00 2.652909254520416e+00 +8.571559056523148e+00 2.650898606589686e+00 +8.589877291409554e+00 2.639144537057880e+00 +8.608282525512919e+00 2.635516013997829e+00 +8.626719350260537e+00 2.632390720700676e+00 +8.645289600751155e+00 2.620700954131104e+00 +8.663965626593621e+00 2.612252413728478e+00 +8.682658157759494e+00 2.614126032661007e+00 +8.701445979197439e+00 2.597593287384506e+00 +8.720359422297363e+00 2.589314261699760e+00 +8.739184840003235e+00 2.579149494857749e+00 +8.758074298993101e+00 2.574244495672498e+00 +8.777160522226643e+00 2.573366210105512e+00 +8.796293606807650e+00 2.556282805880690e+00 +8.815562106304318e+00 2.552268933531098e+00 +8.834835601244366e+00 2.550887044879634e+00 +8.854189564687649e+00 2.536590892216023e+00 +8.873673751054733e+00 2.525899863210910e+00 +8.893248416182921e+00 2.518632394881933e+00 +8.912888405843276e+00 2.519248742264150e+00 +8.932591764325029e+00 2.506690561103154e+00 +8.952435986282028e+00 2.502486189662333e+00 +8.972357961949458e+00 2.493994159567435e+00 +8.992333251432614e+00 2.491282792530783e+00 +9.012440287007630e+00 2.476662753490056e+00 +9.032697918627564e+00 2.472633024315560e+00 +9.052897096278015e+00 2.460776466498469e+00 +9.073244818964366e+00 2.456616760969580e+00 +9.093756280848165e+00 2.444518386032422e+00 +9.114378891040934e+00 2.441329979721503e+00 +9.135028531983528e+00 2.432691945217595e+00 +9.155798816280653e+00 2.420086595613153e+00 +9.176660212753060e+00 2.409398229944333e+00 +9.197612570742654e+00 2.400733412883347e+00 +9.218595974446435e+00 2.397732204255574e+00 +9.239743356473054e+00 2.387762592615275e+00 +9.261033236974651e+00 2.384921872086804e+00 +9.282314054497562e+00 2.369260093385027e+00 +9.303654919984035e+00 2.363409436377947e+00 +9.325219064231849e+00 2.350868014406362e+00 +9.346867437287070e+00 2.348813790451135e+00 +9.368592280481810e+00 2.334661907179259e+00 +9.390461084950465e+00 2.331054483437840e+00 +9.412480441702488e+00 2.324942472378589e+00 +9.434609308465520e+00 2.310021077766146e+00 +9.456787744880494e+00 2.304363196151857e+00 +9.479132859779758e+00 2.291934168302471e+00 +9.501532273866964e+00 2.292906210231002e+00 +9.524074884720834e+00 2.281543996116447e+00 +9.546760175522760e+00 2.270811346196987e+00 +9.569662939507941e+00 2.260849311515461e+00 +9.592664577658210e+00 2.255008733022742e+00 +9.615751972517165e+00 2.245149227232055e+00 +9.638867596873414e+00 2.240448613601230e+00 +9.662121740420160e+00 2.235439636001920e+00 +9.685603315333678e+00 2.214555471204426e+00 +9.709137568744383e+00 2.212331508486050e+00 +9.732771065787128e+00 2.207324221547729e+00 +9.756546505348966e+00 2.194533992775143e+00 +9.780387840606352e+00 2.189828973006383e+00 +9.804419768315384e+00 2.184797523996344e+00 +9.828582133660241e+00 2.173201716131107e+00 +9.852925991717951e+00 2.174129558366568e+00 +9.877389832029518e+00 2.167607235642016e+00 +9.902009145422760e+00 2.154101463863629e+00 +9.926812707367047e+00 2.151189881182799e+00 +9.951668043124464e+00 2.123855781509024e+00 +9.976789344237023e+00 2.128871418600208e+00 +1.000193792146594e+01 2.114329173829787e+00 +1.002728297157961e+01 2.109363668873138e+00 +1.005273624707605e+01 2.095016648319172e+00 +1.007845610131479e+01 2.086317598939071e+00 +1.010425508516770e+01 2.080434477993742e+00 +1.013013614699210e+01 2.070748148102609e+00 +1.015613430196925e+01 2.066100239678576e+00 +1.018245065285365e+01 2.059366254736164e+00 +1.020885369000165e+01 2.048781467833260e+00 +1.023542733232621e+01 2.042676299492451e+00 +1.026215423165389e+01 2.036520048379253e+00 +1.028897659278438e+01 2.035379085874967e+00 +1.031598557409722e+01 2.022940909540991e+00 +1.034325640110906e+01 2.011561093858037e+00 +1.037070561220523e+01 2.012309279605455e+00 +1.039828479406370e+01 1.998446747026668e+00 +1.042601367265439e+01 1.993153061584416e+00 +1.045405236205431e+01 1.987534970893656e+00 +1.048223825780211e+01 1.978372501429631e+00 +1.051062946279310e+01 1.970922835192932e+00 +1.053922431181391e+01 1.964546271531353e+00 +1.056791894659611e+01 1.951547004853847e+00 +1.059681990914897e+01 1.947752319986311e+00 +1.062598184157560e+01 1.938822333494130e+00 +1.065520980325237e+01 1.931693153320327e+00 +1.068475675289162e+01 1.920962954547917e+00 +1.071457381922773e+01 1.920762603351708e+00 +1.074452218913663e+01 1.908283801870663e+00 +1.077470044587385e+01 1.899828093647898e+00 +1.080529519883718e+01 1.898026089907203e+00 +1.083591686155712e+01 1.882638782054732e+00 +1.086687059750680e+01 1.870984119630311e+00 +1.089800006536714e+01 1.861767050428116e+00 +1.092924936189662e+01 1.857769025796347e+00 +1.096069152773907e+01 1.852275104329806e+00 +1.099262670006797e+01 1.839674948083331e+00 +1.102482457585336e+01 1.830886593173797e+00 +1.105731619871706e+01 1.827033737341105e+00 +1.108987707187318e+01 1.817521594806416e+00 +1.112283089469190e+01 1.809485305285870e+00 +1.115606970854127e+01 1.799488825785699e+00 +1.118933945977239e+01 1.788087517568392e+00 +1.122285940573241e+01 1.787965691024092e+00 +1.125678966553714e+01 1.774239684881761e+00 +1.129109968243578e+01 1.764341792974305e+00 +1.132555474510823e+01 1.755567880755174e+00 +1.136046582409104e+01 1.750230401953144e+00 +1.139574483886010e+01 1.743240182994079e+00 +1.143106978440181e+01 1.735355423316995e+00 +1.146678606872380e+01 1.725598780373527e+00 +1.150288296881248e+01 1.718516113266582e+00 +1.153920715822614e+01 1.709753488661443e+00 +1.157604626834376e+01 1.709350723883968e+00 +1.161312028935976e+01 1.693667994960284e+00 +1.165028012742078e+01 1.678968553552141e+00 +1.168793289059454e+01 1.670550296975824e+00 +1.172610471209537e+01 1.657131024283915e+00 +1.176459260492638e+01 1.651268961350075e+00 +1.180336772075291e+01 1.636610141465336e+00 +1.184249359904616e+01 1.634491817281363e+00 +1.188194271004327e+01 1.619572033615788e+00 +1.192174979859032e+01 1.613905069167474e+00 +1.196209700035650e+01 1.612073607623486e+00 +1.200286471269033e+01 1.600197414383230e+00 +1.204403216146022e+01 1.589787025590657e+00 +1.208569449561884e+01 1.575241383432848e+00 +1.212772406981022e+01 1.561841496959909e+00 +1.217031748369286e+01 1.559056543964437e+00 +1.221331312560338e+01 1.544762360955713e+00 +1.225676792866482e+01 1.533449820590916e+00 +1.230064979482002e+01 1.524497863416844e+00 +1.234476550651174e+01 1.515097883130812e+00 +1.238961073579248e+01 1.507004397218803e+00 +1.243487468792004e+01 1.500061720460831e+00 +1.248066924011938e+01 1.500513033442246e+00 +1.252708843316362e+01 1.490386995264767e+00 +1.257414997469930e+01 1.488993419643792e+00 +1.262182492035476e+01 1.479529860339208e+00 +1.267009493352072e+01 1.468891328497975e+00 +1.271888368754871e+01 1.461499892986887e+00 +1.276828246118984e+01 1.451402354733088e+00 +1.281830420042961e+01 1.441069267030485e+00 +1.286888912042902e+01 1.434688492527573e+00 +1.291996260858576e+01 1.417643078485916e+00 +1.297190846976042e+01 1.406795643415921e+00 +1.302488385375759e+01 1.392432252907555e+00 +1.307820478032981e+01 1.377530938665278e+00 +1.313245303531817e+01 1.370448188327311e+00 +1.318736319606146e+01 1.358155264497709e+00 +1.324333171191776e+01 1.342321935095090e+00 +1.329986364716511e+01 1.333227122926989e+00 +1.335731802367998e+01 1.320974574249329e+00 +1.341562495807058e+01 1.310530892775097e+00 +1.347485234240245e+01 1.296799107859835e+00 +1.353509525855492e+01 1.286541354450607e+00 +1.359626006079933e+01 1.270326517075242e+00 +1.365843569018352e+01 1.249569524456653e+00 +1.372165680261726e+01 1.237734932380492e+00 +1.378589719320418e+01 1.213727944928627e+00 +1.385100535214515e+01 1.205980157063200e+00 +1.391747309463467e+01 1.192863578523194e+00 +1.398506845795805e+01 1.184789098814713e+00 +1.405363080218904e+01 1.178242093925171e+00 +1.412350505351254e+01 1.164452884725742e+00 +1.419485363793802e+01 1.159511171393126e+00 +1.426747025743661e+01 1.146265267693796e+00 +1.434142460993240e+01 1.134390556093717e+00 +1.441708047808282e+01 1.125565268863087e+00 +1.449435250885427e+01 1.115198215173279e+00 +1.457308480943307e+01 1.108006697975223e+00 +1.465329819966831e+01 1.089002051815754e+00 +1.473531320186686e+01 1.078620070673199e+00 +1.481914971229275e+01 1.061106606245876e+00 +1.490494291923437e+01 1.046010494406160e+00 +1.499293084576228e+01 1.028232627334566e+00 +1.508298399066317e+01 1.010213891705542e+00 +1.517541065745850e+01 1.008709912203584e+00 +1.527008219051617e+01 9.961635236125702e-01 +1.536745062502385e+01 9.903182559964580e-01 +1.546714308265249e+01 9.889532951294517e-01 +1.556952555032566e+01 9.866554248245296e-01 +1.567451843087067e+01 9.751022739757471e-01 +1.578258318133163e+01 9.581752004806124e-01 +1.589400105194650e+01 9.428478100725186e-01 +1.600883576639662e+01 9.315844261645112e-01 +1.612739873627202e+01 9.145155214678200e-01 +1.624994892403645e+01 9.033965869332343e-01 +1.637672258877143e+01 8.903349831849883e-01 +1.650779077679548e+01 8.815423628951109e-01 +1.664418992799417e+01 8.783163735036886e-01 +1.678554308959072e+01 8.616450519516042e-01 +1.693258133573742e+01 8.430254720373530e-01 +1.708569091866682e+01 8.197364551275157e-01 +1.724545509659271e+01 7.985487295339642e-01 +1.741298582570438e+01 7.731383750842359e-01 +1.758826809793949e+01 7.532898987246466e-01 +1.777225827066799e+01 7.324988226802176e-01 +1.796765617555977e+01 7.068536024157870e-01 +1.817428149684734e+01 6.998162255020254e-01 +1.839271680267985e+01 6.859530057856070e-01 +1.862568646838732e+01 6.774751369416968e-01 +1.887587394900956e+01 6.754708661651644e-01 +1.914386953833550e+01 6.455320858921243e-01 +1.943304403340177e+01 6.213828327528838e-01 +1.974863831988926e+01 5.979866618972404e-01 +2.009461352598835e+01 5.886938310784255e-01 +2.047803223982858e+01 5.757747375246646e-01 +2.090855495752863e+01 5.539528352816189e-01 +2.139949714153402e+01 5.102306067839727e-01 +2.197195509896472e+01 5.050121141028469e-01 +2.265608929794178e+01 4.457184212285230e-01 +2.351362565544672e+01 4.134252642197285e-01 +2.466187871583055e+01 3.929338773153763e-01 +2.642192176286753e+01 3.509561515969558e-01 +3.101828541027504e+01 2.806767148869000e-01 diff --git a/example/test_problem/ELBDM/DiskHeating/Input__Flag_NParPatch b/example/test_problem/ELBDM/DiskHeating/Input__Flag_NParPatch new file mode 100644 index 0000000000..b0f491a5e8 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/Input__Flag_NParPatch @@ -0,0 +1,12 @@ +# Level Number_of_particles_per_patch + 0 200 + 1 200 + 2 200 + 3 200 + 4 200 + 5 200 + 6 200 + 7 200 + 8 200 + 9 200 + 10 200 diff --git a/example/test_problem/ELBDM/DiskHeating/Input__Flag_Rho b/example/test_problem/ELBDM/DiskHeating/Input__Flag_Rho new file mode 100644 index 0000000000..41cc130a80 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/Input__Flag_Rho @@ -0,0 +1,5 @@ +# Level Density + 0 1.0e+3 + 1 7.0e+3 + 2 1.0e+5 + 3 4.0e+5 diff --git a/example/test_problem/ELBDM/DiskHeating/Input__Parameter b/example/test_problem/ELBDM/DiskHeating/Input__Parameter new file mode 100644 index 0000000000..fafe3b825d --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/Input__Parameter @@ -0,0 +1,318 @@ + + +# ================================================================================================================= +# NOTE: +# 1. Comment symbol: # +# 2. [*]: defaults +# 3. Parameters set to "auto" (usually by setting to a negative value) do not have deterministic default values +# and will be set according to the adopted compilation options and/or other runtime parameters +# 4. To add new parameters, please edit "Init/Init_Load_Parameter.cpp" +# 5. All dimensional variables should be set consistently with the code units (set by UNIT_L/M/T/V/D) unless +# otherwise specified (e.g., SF_CREATE_STAR_MIN_GAS_DENS & SF_CREATE_STAR_MIN_STAR_MASS) +# 6. For boolean options: 0/1 -> off/on +# ================================================================================================================= + + +# simulation scale +BOX_SIZE 1.60671536e-01 # box size along the longest side (in Mpc/h if COMOVING is adopted) +NX0_TOT_X 480 # number of base-level cells along x +NX0_TOT_Y 480 # number of base-level cells along y +NX0_TOT_Z 480 # number of base-level cells along z +OMP_NTHREAD -1 # number of OpenMP threads (<=0=auto) [-1] ##OPENMP ONLY## +END_T 2.5e-1 # end physical time (<0=auto -> must be set by test problems or restart) [-1.0] +END_STEP -1 # end step (<0=auto -> must be set by test problems or restart) [-1] + + +# test problems +TESTPROB_ID 1013 # test problem ID [0] + # 1003: ELBDM soliton merger (+GRAVITY) + +# code units (in cgs) +OPT__UNIT 1 # specify code units -> must set exactly 3 basic units below [0] ##USELESS FOR COMOVING## +UNIT_L 4.436632034507548e+24 # length unit (<=0 -> set to UNIT_V*UNIT_T or (UNIT_M/UNIT_D)^(1/3)) [-1.0] +UNIT_M -1.0 # mass unit (<=0 -> set to UNIT_D*UNIT_L^3) [-1.0] +UNIT_T -1.0 # time unit (<=0 -> set to UNIT_L/UNIT_V) [-1.0] +UNIT_V 1.00000000000000e+07 # velocity unit (<=0 -> set to UNIT_L/UNIT_T) [-1.0] +UNIT_D 2.5758579703106658e-30 # mass density unit (<=0 -> set to UNIT_M/UNIT_L^3) [-1.0] + + +# boundary conditions +OPT__BC_FLU_XM 4 # fluid boundary condition at the -x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_XP 4 # fluid boundary condition at the +x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YM 4 # fluid boundary condition at the -y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YP 4 # fluid boundary condition at the +y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZM 4 # fluid boundary condition at the -z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZP 4 # fluid boundary condition at the +z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_POT 2 # gravity boundary condition: (1=periodic, 2=isolated) +GFUNC_COEFF0 -1.0 # Green's function coefficient at the origin for the isolated BC (<0=auto) [-1.0] + + +# particle (PARTICLE only) +PAR_NPAR 20000000 # total number of particles (must be set for PAR_INIT==1/3; must be an integer) +PAR_INIT 1 # initialization option for particles: (1=FUNCTION, 2=RESTART, 3=FILE->"PAR_IC") +PAR_IC_FORMAT 1 # data format of PAR_IC: (1=[attribute][id], 2=[id][attribute]; row-major) [1] +PAR_IC_FLOAT8 -1 # floating-point precision for PAR_IC (<0: default, 0: single, 1: double) [default: same as FLOAT8_PAR] +PAR_IC_MASS -1.0 # mass of all particles for PAR_INIT==3 (<0=off) [-1.0] +PAR_IC_TYPE -1 # type of all particles for PAR_INIT==3 (<0=off) [-1] +PAR_INTERP 3 # particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [2] +PAR_INTEG 2 # particle integration scheme: (1=Euler, 2=KDK) [2] +PAR_TR_INTERP 3 # tracer particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [3] +PAR_TR_INTEG 2 # tracer particle integration scheme: (1=Euler, 2=RK2) [2] +PAR_IMPROVE_ACC 1 # improve force accuracy at patch boundaries [1] ##STORE_POT_GHOST and PAR_INTERP=2/3 ONLY## +PAR_PREDICT_POS 1 # predict particle position during mass assignment [1] +PAR_REMOVE_CELL -1.0 # remove particles X-root-cells from the boundaries (non-periodic BC only; <0=auto) [-1.0] +OPT__FREEZE_PAR 0 # do not update particles (except for tracers) [0] +PAR_TR_VEL_CORR 0 # correct tracer particle velocities in regions of discontinuous flow [0] + +# time-step +DT__MAX -1.0 # dt criterion: maximum allowed dt (<0=off) [-1.0] +DT__FLUID -1.0 # dt criterion: fluid solver CFL factor (<0=auto) [-1.0] +DT__FLUID_INIT -1.0 # dt criterion: DT__FLUID at the first step (<0=auto) [-1.0] +DT__SPEED_OF_LIGHT 0 # dt criterion: speed of light [0] ##SRHD ONLY## +DT__GRAVITY -1.0 # dt criterion: gravity solver safety factor (<0=auto) [-1.0] +DT__PHASE 0.0 # dt criterion: phase rotation safety factor (0=off) [0.0] ##ELBDM ONLY## +DT__HYBRID_CFL -1.0 # dt criterion: hybrid solver CFL factor (<0=auto) (diffusion) [-1.0] ## ELBDM_HYBRID ONLY## +DT__HYBRID_CFL_INIT -1.0 # dt criterion: DT__HYBRID_CFL in the first step (<0=auto) [-1.0] ## ELBDM_HYBRID ONLY## +DT__HYBRID_VELOCITY -1.0 # dt criterion: hybrid solver CFL factor (<0=auto) (Hamilton-Jacobi) [-1.0] ## ELBDM_HYBRID ONLY## +DT__HYBRID_VELOCITY_INIT -1.0 # dt criterion: DT__HYBRID_VELOCITY in the first step (<0=auto) [-1.0] ## ELBDM_HYBRID ONLY## +DT__PARVEL 0.5 # dt criterion: particle velocity safety factor [0.5] +DT__PARVEL_MAX -1.0 # dt criterion: maximum allowed dt from particle velocity (<0=off) [-1.0] +DT__PARACC 0.5 # dt criterion: particle acceleration safety factor (0=off) [0.5] ##STORE_PAR_ACC ONLY## +DT__CR_DIFFUSION 0.3 # dt criterion: cosmic-ray diffusion CFL factor [0.3] ##CR_DIFFUSION only## +DT__MAX_DELTA_A 0.01 # dt criterion: maximum variation of the cosmic scale factor [0.01] +DT__SYNC_PARENT_LV 0.1 # dt criterion: allow dt to adjust by (1.0+DT__SYNC_PARENT) in order to synchronize + # with the parent level (for OPT__DT_LEVEL==3 only) [0.1] +DT__SYNC_CHILDREN_LV 0.1 # dt criterion: allow dt to adjust by (1.0-DT__SYNC_CHILDREN) in order to synchronize + # with the children level (for OPT__DT_LEVEL==3 only; 0=off) [0.1] +OPT__DT_USER 0 # dt criterion: user-defined -> edit "Mis_GetTimeStep_UserCriteria.cpp" [0] +OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] +OPT__RECORD_DT 1 # record info of the dt determination [1] +AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] +AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] +AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] +AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## +AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## +AUTO_REDUCE_INT_MONO_FACTOR 0.8 # reduce INT_MONO_COEFF(_B) by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] +AUTO_REDUCE_INT_MONO_MIN 1.0e-2 # minimum allowed INT_MONO_COEFF(_B) after consecutive failures [1.0e-2] + + +# grid refinement (examples of Input__Flag_XXX tables are put at "example/input/") +REGRID_COUNT 4 # refine every REGRID_COUNT sub-step [4] +REFINE_NLEVEL 1 # number of new AMR levels to be created at once during refinement [1] +FLAG_BUFFER_SIZE 8 # number of buffer cells for the flag operation (0~PATCH_SIZE; <0=auto -> PATCH_SIZE) [-1] +FLAG_BUFFER_SIZE_MAXM1_LV 4 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-1 (<0=auto -> REGRID_COUNT) [-1] +FLAG_BUFFER_SIZE_MAXM2_LV -1 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-2 (<0=auto) [-1] +MAX_LEVEL 2 # maximum refinement level (0~NLEVEL-1) [NLEVEL-1] +OPT__FLAG_RHO 1 # flag: density (Input__Flag_Rho) [0] +OPT__FLAG_RHO_GRADIENT 0 # flag: density gradient (Input__Flag_RhoGradient) [0] +OPT__FLAG_PRES_GRADIENT 0 # flag: pressure gradient (Input__Flag_PresGradient) [0] ##HYDRO ONLY## +OPT__FLAG_LRTZ_GRADIENT 0 # flag: Lorentz factor gradient (Input__Flag_LrtzGradient) [0] ##SRHD ONLY## +OPT__FLAG_VORTICITY 0 # flag: vorticity (Input__Flag_Vorticity) [0] ##HYDRO ONLY## +OPT__FLAG_JEANS 0 # flag: Jeans length (Input__Flag_Jeans) [0] ##HYDRO ONLY## +OPT__FLAG_CURRENT 0 # flag: current density in MHD (Input__Flag_Current) [0] ##MHD ONLY## +OPT__FLAG_CRAY 0 # flag: cosmic-ray energy (Input__Flag_CRay) [0] ##COSMIC_RAY ONLY## +OPT__FLAG_ENGY_DENSITY 0 # flag: energy density (Input_Flag_EngyDensity) [0] ##ELBDM ONLY## +OPT__FLAG_INTERFERENCE 0 # flag: interference level (Input__Flag_Interference) [0] ##ELBDM ONLY## +OPT__FLAG_SPECTRAL 0 # flag: spectral refinement (Input__Flag_Spectral) [0] ##ELBDM ONLY## +OPT__FLAG_LOHNER_DENS 0 # flag: Lohner for mass density (Input__Flag_Lohner) [0] ##BOTH HYDRO AND ELBDM## +OPT__FLAG_LOHNER_ENGY 0 # flag: Lohner for energy density (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_PRES 0 # flag: Lohner for pressure (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_TEMP 0 # flag: Lohner for temperature (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_ENTR 0 # flag: Lohner for entropy (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_CRAY 0 # flag: Lohner for cosmic-ray energy (Input__Flag_Lohner) [0] ##COSMIC_RAY ONLY## +OPT__FLAG_LOHNER_FORM 4 # form of Lohner: (1=FLASH-1, 2=FLASH-2, 3=form-invariant-1, 4=form-invariant-2) [2] +OPT__FLAG_USER 0 # flag: user-defined (Input__Flag_User) -> edit "Flag_User.cpp" [0] +OPT__FLAG_USER_NUM 1 # number of threshold values in user-defined table (Input__Flag_User) [1] +OPT__FLAG_REGION 0 # flag: specify the regions **allowed** to be refined -> edit "Flag_Region.cpp" [0] +OPT__FLAG_NPAR_PATCH 2 # flag: # of particles per patch (Input__Flag_NParPatch): (0=off, 1=itself, 2=itself+siblings) [0] +OPT__FLAG_NPAR_CELL 0 # flag: # of particles per cell (Input__Flag_NParCell) [0] +OPT__FLAG_PAR_MASS_CELL 0 # flag: total particle mass per cell (Input__Flag_ParMassCell) [0] +OPT__NO_FLAG_NEAR_BOUNDARY 0 # flag: disallow refinement near the boundaries [0] +OPT__PATCH_COUNT 1 # record the # of patches at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__PARTICLE_COUNT 1 # record the # of particles at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__REUSE_MEMORY 2 # reuse patch memory to reduce memory fragmentation: (0=off, 1=on, 2=aggressive) [2] +OPT__MEMORY_POOL 0 # preallocate patches for OPT__REUSE_MEMORY=1/2 (Input__MemoryPool) [0] + + +# load balance (LOAD_BALANCE only) +LB_INPUT__WLI_MAX 0.1 # weighted-load-imbalance (WLI) threshold for redistributing all patches [0.1] +LB_INPUT__PAR_WEIGHT 1.0 # load-balance weighting of one particle over one cell [0.0] +OPT__RECORD_LOAD_BALANCE 1 # record the load-balance info [1] +OPT__MINIMIZE_MPI_BARRIER 1 # minimize MPI barriers to improve load balance, especially with particles [0] + # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0) + +# fluid solver in ELBDM (MODEL==ELBDM only) +ELBDM_MASS 4.0e-23 # particle mass in ev/c^2 (input unit is fixed even when OPT__UNIT or COMOVING is on) +ELBDM_PLANCK_CONST 1.0 # reduced Planck constant (will be overwritten if OPT__UNIT or COMOVING is on) +ELBDM_LAMBDA 1.0 # quartic self-interaction coefficient [1.0] ##QUARTIC_SELF_INTERACTION ONLY## +ELBDM_TAYLOR3_COEFF 0.166666667 # 3rd Taylor expansion coefficient [1.0/6.0] ##USELESS if ELBDM_TAYLOR3_AUTO is on## +ELBDM_TAYLOR3_AUTO 0 # Optimize ELBDM_TAYLOR3_COEFF automatically to minimize the damping at kmax [0] +ELBDM_REMOVE_MOTION_CM 0 # remove the motion of center-of-mass (must enable OPT__CK_CONSERVATION): + # (0=off, 1=init, 2=every step) [0] +ELBDM_BASE_SPECTRAL 0 # adopt the spectral method to evolve base-level wave function (must enable SUPPORT_FFTW) [0] +ELBDM_MATCH_PHASE 1 # match child phases with father phases during data restriction [1] ##ELBDM_HYBRID ONLY## +ELBDM_FIRST_WAVE_LEVEL -1 # level at which to switch to the wave solver (must >=1) [-1] ##ELBDM_HYBRID ONLY## + +# fluid solvers in all models +FLU_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU fluid solver (<=0=auto) [-1] +GPU_NSTREAM -1 # number of CUDA streams for the asynchronous memory copy in GPU (<=0=auto) [-1] +OPT__FIXUP_FLUX 1 # correct coarse grids by the fine-grid boundary fluxes [1] ##HYDRO and ELBDM ONLY## +OPT__FIXUP_ELECTRIC 1 # correct coarse grids by the fine-grid boundary electric field [1] ##MHD ONLY## +OPT__FIXUP_RESTRICT 1 # correct coarse grids by averaging the fine-grid data [1] +OPT__CORR_AFTER_ALL_SYNC -1 # apply various corrections after all levels are synchronized (see "Flu_CorrAfterAllSync"): + # (-1=auto, 0=off, 1=every step, 2=before dump) [-1] +OPT__NORMALIZE_PASSIVE 0 # ensure "sum(passive_scalar_density) == gas_density" [1] +OPT__INT_FRAC_PASSIVE_LR 1 # convert specified passive scalars to mass fraction during data reconstruction [1] +OPT__OVERLAP_MPI 0 # overlap MPI communication with CPU/GPU computations [0] ##NOT SUPPORTED YET## +OPT__RESET_FLUID 0 # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0] +OPT__RESET_FLUID_INIT -1 # reset fluid variables during initialization (<0=auto -> OPT__RESET_FLUID, 0=off, 1=on) [-1] +OPT__FREEZE_FLUID 0 # do not evolve fluid at all [0] +OPT__CHECK_PRES_AFTER_FLU -1 # check unphysical pressure at the end of the fluid solver (<0=auto) [-1] +OPT__LAST_RESORT_FLOOR 1 # apply floor values as the last resort when the fluid solver fails [1] ##HYDRO and MHD ONLY## +MIN_DENS 0.0 # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY## +MIN_PRES 0.0 # minimum pressure (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_EINT 0.0 # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_TEMP 0.0 # minimum temperature in K (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_ENTR 0.0 # minimum entropy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +JEANS_MIN_PRES 0 # minimum pressure estimated from the Jeans length [0] ##HYDRO/MHD and GRAVITY ONLY## +JEANS_MIN_PRES_LEVEL -1 # for JEANS_MIN_PRES; ensure Jeans length is resolved by JEANS_MIN_PRES_NCELL*dh[JEANS_MIN_PRES_LEVEL] + # (<0=auto -> MAX_LEVEL) [-1] +JEANS_MIN_PRES_NCELL 4 # for JEANS_MIN_PRES; see JEANS_MIN_PRES_LEVEL [4] + + +# gravity solvers in all models +NEWTON_G 1.0 # gravitational constant (will be overwritten if OPT__UNIT or COMOVING is on) +SOR_OMEGA -1.0 # over-relaxation parameter in SOR: (<0=auto) [-1.0] +SOR_MAX_ITER -1 # maximum number of iterations in SOR: (<0=auto) [-1] +SOR_MIN_ITER -1 # minimum number of iterations in SOR: (<0=auto) [-1] +MG_MAX_ITER -1 # maximum number of iterations in multigrid: (<0=auto) [-1] +MG_NPRE_SMOOTH -1 # number of pre-smoothing steps in multigrid: (<0=auto) [-1] +MG_NPOST_SMOOTH -1 # number of post-smoothing steps in multigrid: (<0=auto) [-1] +MG_TOLERATED_ERROR -1.0 # maximum tolerated error in multigrid (<0=auto) [-1.0] +POT_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Poisson solver (<=0=auto) [-1] +OPT__GRA_P5_GRADIENT 0 # 5-points gradient in the Gravity solver (must have GRA/USG_GHOST_SIZE_G>=2) [0] +OPT__SELF_GRAVITY 1 # add self-gravity [1] +OPT__EXT_ACC 0 # add external acceleration (0=off, 1=function, 2=table) [0] ##HYDRO ONLY## + # --> 2 (table) is not supported yet +OPT__EXT_POT 0 # add external potential (0=off, 1=function, 2=table) [0] + # --> for 2 (table), edit the corresponding parameters below too + +# initialization +OPT__INIT 3 # initialization option: (1=FUNCTION, 2=RESTART, 3=FILE->"UM_IC") +OPT__INIT_BFIELD_BYVECPOT 0 # initialize the magnetic field from vector potential + # (0=off, 1=external disk file named "B_IC", see tool/inits/gen_vec_pot.py for example, 2=function) [0] ##MHD ONLY## +RESTART_LOAD_NRANK 4 # number of parallel I/O (i.e., number of MPI ranks) for restart [1] +OPT__RESTART_RESET 0 # reset some simulation status parameters (e.g., current step and time) during restart [0] +OPT__UM_IC_LEVEL 0 # AMR level corresponding to UM_IC (must >= 0) [0] +OPT__UM_IC_NLEVEL 3 # number of AMR levels UM_IC [1] --> edit "Input__UM_IC" if >1 +OPT__UM_IC_NVAR 2 # number of variables in UM_IC: (1~NCOMP_TOTAL; <=0=auto) [HYDRO=5+passive/ELBDM=2] +OPT__UM_IC_FORMAT 1 # data format of UM_IC: (1=vzyx, 2=zyxv; row-major and v=field) [1] +OPT__UM_IC_FLOAT8 -1 # floating-point precision for UM_IC (<0: default, 0: single, 1: double) [default: same as FLOAT8] +OPT__UM_IC_DOWNGRADE 1 # downgrade UM_IC from level OPT__UM_IC_LEVEL to 0 [1] +OPT__UM_IC_REFINE 1 # refine UM_IC from level OPT__UM_IC_LEVEL to MAX_LEVEL [1] +OPT__UM_IC_LOAD_NRANK 4 # number of parallel I/O (i.e., number of MPI ranks) for loading UM_IC [1] +OPT__INIT_RESTRICT 1 # restrict all data during the initialization [1] +OPT__INIT_GRID_WITH_OMP 1 # enable OpenMP when assigning the initial condition of each grid patch [1] +OPT__GPUID_SELECT -1 # GPU ID selection mode: (-3=Laohu, -2=CUDA, -1=MPI rank, >=0=input) [-1] +INIT_SUBSAMPLING_NCELL 0 # perform sub-sampling during initialization: (0=off, >0=# of sub-sampling cells) [0] +OPT__FFTW_STARTUP -1 # initialise fftw plans: (-1=auto, 0=ESTIMATE, 1=MEASURE, 2=PATIENT (only FFTW3)) [-1] + +# interpolation schemes: (-1=auto, 1=MinMod-3D, 2=MinMod-1D, 3=vanLeer, 4=CQuad, 5=Quad, 6=CQuar, 7=Quar, 8=Spectral (##ELBDM & SUPPORT_SPECTRAL_INT ONLY##)) +OPT__INT_TIME 1 # perform "temporal" interpolation for OPT__DT_LEVEL == 2/3 [1] +OPT__INT_PRIM 1 # switch to primitive variables when the interpolation on conserved variables fails [1] ##HYDRO ONLY## +OPT__INT_PHASE 1 # interpolation on phase (does not support MinMod-1D) [1] ##ELBDM ONLY## +OPT__RES_PHASE 0 # restriction on phase [0] ##ELBDM ONLY## +OPT__FLU_INT_SCHEME -1 # ghost-zone fluid variables for the fluid solver [-1] +OPT__REF_FLU_INT_SCHEME -1 # newly allocated fluid variables during grid refinement [-1] +OPT__MAG_INT_SCHEME 4 # ghost-zone magnetic field for the MHD solver (2,3,4,6 only) [4] +OPT__REF_MAG_INT_SCHEME 4 # newly allocated magnetic field during grid refinement (2,3,4,6 only) [4] +OPT__POT_INT_SCHEME 5 # ghost-zone potential for the Poisson solver (only supports 4 & 5) [4] +OPT__RHO_INT_SCHEME 4 # ghost-zone mass density for the Poisson solver [4] +OPT__GRA_INT_SCHEME 5 # ghost-zone potential for the gravity solver (for UNSPLIT_GRAVITY as well) [4] +OPT__REF_POT_INT_SCHEME 5 # newly allocated potential during grid refinement [4] +INT_MONO_COEFF 2.0 # coefficient for ensuring the interpolation monotonicity (1.0~4.0) [2.0] +INT_MONO_COEFF_B 2.0 # coefficient for ensuring the interpolation monotonicity of B field (1.0~4.0) [2.0] ##MHD ONLY## +MONO_MAX_ITER 10 # maximum number of iterations to reduce INT_MONO_COEFF when interpolation fails (0=off) [10] +INT_OPP_SIGN_0TH_ORDER 0 # switch to 0th-order interpolation if adjacent values change signs [HYDRO:1; ELBDM:0] +SPEC_INT_TABLE_PATH ./ # path to tables for spectral interpolation ##ELBDM & SUPPORT_SPECTRAL_INT ONLY## +SPEC_INT_XY_INSTEAD_DEPHA 1 # interpolate x = density^0.5*cos( phase/SPEC_INT_WAVELENGTH_MAGNIFIER ), + # y = density^0.5*sin( phase/SPEC_INT_WAVELENGTH_MAGNIFIER ) instead of density and phase + # for the spectral interpolation, which has the advantage of being well-defined across vortices [1] +SPEC_INT_WAVELENGTH_MAGNIFIER 1.0e2 # stretching factor of wavelength for SPEC_INT_XY_INSTEAD_DEPHA + # --> for SPEC_INT_WAVELENGTH_MAGNIFIER=1, x=real part and y=imaginary part [1.0e2] + + +# data dump +OPT__OUTPUT_TOTAL 1 # output the simulation snapshot: (0=off, 1=HDF5, 2=C-binary) [1] +OPT__OUTPUT_PART 0 # output a single line or slice: (0=off, 1=xy, 2=yz, 3=xz, 4=x, 5=y, 6=z, 7=diag) [0] +OPT__OUTPUT_TEXT_FORMAT_FLT %24.16e # string format of output text files [%24.16e] +OPT__OUTPUT_USER 0 # output the user-specified data -> edit "Output_User.cpp" [0] +OPT__OUTPUT_PAR_MODE 0 # output the particle data: (0=off, 1=text-file, 2=C-binary) [0] ##PARTICLE ONLY## +OPT__OUTPUT_BASEPS 0 # output the base-level power spectrum [0] +OPT__OUTPUT_BASE 0 # only output the base-level data [0] ##OPT__OUTPUT_PART ONLY## +OPT__OUTPUT_POT 1 # output gravitational potential [1] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_PAR_DENS 1 # output the particle or total mass density on grids: + # (0=off, 1=particle mass density, 2=total mass density) [1] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_CC_MAG 1 # output **cell-centered** magnetic field (necessary for yt analysis) [1] ##MHD ONLY## +OPT__OUTPUT_PRES 0 # output gas pressure [0] ##HYDRO ONLY## +OPT__OUTPUT_TEMP 0 # output gas temperature [0 (HD) or 1 (SRHD)] ##HYDRO ONLY## +OPT__OUTPUT_ENTR 0 # output gas entropy [0] ##HYDRO ONLY## +OPT__OUTPUT_CS 0 # output sound speed [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVVEL 0 # output divergence(velocity) [0] ##HYDRO ONLY## +OPT__OUTPUT_MACH 0 # output mach number [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVMAG 0 # output |divergence(B)*dh/|B|| [0] ##MHD ONLY## +OPT__OUTPUT_LORENTZ 0 # output Lorentz factor [0] ##SRHD ONLY## +OPT__OUTPUT_3VELOCITY 0 # output 3-velocities [0] ##SRHD ONLY## +OPT__OUTPUT_ENTHALPY 1 # output reduced enthalpy [1] ##SRHD ONLY## +OPT__OUTPUT_USER_FIELD 0 # output user-defined derived fields [0] -> edit "Flu_DerivedField_User.cpp" +OPT__OUTPUT_MODE 2 # (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 +OPT__OUTPUT_RESTART 0 # output data immediately after restart [0] +OUTPUT_STEP 40 # output data every OUTPUT_STEP step ##OPT__OUTPUT_MODE==1 ONLY## +OUTPUT_DT 1.0e-2 # output data every OUTPUT_DT time interval ##OPT__OUTPUT_MODE==2 ONLY## +OUTPUT_WALLTIME -1.0 # output data every OUTPUT_WALLTIME walltime (<=0.0=off) [-1.0] +OUTPUT_WALLTIME_UNIT 0 # unit of OUTPUT_WALLTIME (0=second, 1=minute, 2=hour, 3=day) [0] +OUTPUT_PART_X -1.0 # x coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Y -1.0 # y coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Z -1.0 # z coordinate for OPT__OUTPUT_PART [-1.0] +INIT_DUMPID -1 # set the first dump ID (<0=auto) [-1] + + +# miscellaneous +OPT__VERBOSE 0 # output the simulation progress in detail [0] +OPT__TIMING_BARRIER -1 # synchronize before timing -> more accurate, but may slow down the run (<0=auto) [-1] +OPT__TIMING_BALANCE 0 # record the max/min elapsed time in various code sections for checking load balance [0] +OPT__TIMING_MPI 0 # record the MPI bandwidth achieved in various code sections [0] ##LOAD_BALANCE ONLY## +OPT__RECORD_NOTE 1 # take notes for the general simulation info [1] +OPT__RECORD_UNPHY 1 # record the number of cells with unphysical results being corrected [1] +OPT__RECORD_MEMORY 1 # record the memory consumption [1] +OPT__RECORD_PERFORMANCE 1 # record the code performance [1] +OPT__MANUAL_CONTROL 1 # support manually dump data, stop run, or pause run during the runtime + # (by generating the file DUMP_GAMER_DUMP, STOP_GAMER_STOP, PAUSE_GAMER_PAUSE, respectively) [1] +OPT__RECORD_CENTER 1 # record the position of maximum density, minimum potential, and center of mass [0] +COM_CEN_X -1.0 # x coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position x) [-1.0] +COM_CEN_Y -1.0 # y coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position y) [-1.0] +COM_CEN_Z -1.0 # z coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position z) [-1.0] +COM_MAX_R -1.0 # maximum radius for determining center of mass (<0=auto -> __FLT_MAX__) [-1.0] +COM_MIN_RHO 0.0 # minimum density for determining center of mass (must >= 0.0) [0.0] +COM_TOLERR_R -1.0 # maximum tolerated error of deviation in radius during the iterations of determining the center of mass (<0=auto -> amr->dh[MAX_LEVEL]) [-1.0] +COM_MAX_ITER 10 # maximum number of iterations for determining the center of mass (must >= 1) [10] +OPT__RECORD_USER 0 # record the user-specified info -> edit "Aux_Record_User.cpp" [0] +OPT__OPTIMIZE_AGGRESSIVE 0 # apply aggressive optimizations (experimental) [0] +OPT__SORT_PATCH_BY_LBIDX 1 # sort patches to improve bitwise reproducibility [SERIAL:0, LOAD_BALACNE:1] + + +# checks +OPT__CK_REFINE 0 # check the grid refinement [0] +OPT__CK_PROPER_NESTING 0 # check the proper-nesting condition [0] +OPT__CK_CONSERVATION 1 # check the conservation law [0] +ANGMOM_ORIGIN_X -1.0 # x coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +ANGMOM_ORIGIN_Y -1.0 # y coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +ANGMOM_ORIGIN_Z -1.0 # z coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +OPT__CK_NORMALIZE_PASSIVE 0 # check the normalization of passive scalars [0] ##OPT__NORMALIZE_PASSIVE ONLY## +OPT__CK_RESTRICT 0 # check the data restriction [0] +OPT__CK_FINITE 0 # check if all variables are finite [0] +OPT__CK_PATCH_ALLOCATE 0 # check if all patches are properly allocated [0] +OPT__CK_FLUX_ALLOCATE 0 # check if all flux arrays are properly allocated [0] ##HYDRO and ELBDM ONLY## +OPT__CK_NEGATIVE 0 # check the negative values: (0=off, 1=density, 2=pressure and entropy, 3=both) [0] ##HYDRO ONLY## +OPT__CK_MEMFREE 1.0 # check the free memory in GB (0=off, >0=threshold) [1.0] +OPT__CK_PARTICLE 0 # check the particle allocation [0] diff --git a/example/test_problem/ELBDM/DiskHeating/Input__TestProb b/example/test_problem/ELBDM/DiskHeating/Input__TestProb new file mode 100644 index 0000000000..57c1593806 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/Input__TestProb @@ -0,0 +1,31 @@ +# center of the halo (only useful when AddFixedHalo=1 or (AddParWhenRestart=1 and AddParWhenRestartByFile=0)) +CenX 7.9770907e-02 # center of the halo +CenY 7.9854590e-02 # center of the halo +CenZ 8.0858787e-02 # center of the halo + +# parameters for background fixed halo (disabled in this example) +AddFixedHalo 0 # add a fixed halo [0] ## must set OPT__INIT=1 and OPT__FREEZE_FLUID=1 ## +HaloUseTable 0 # use density table for the fixed halo [0] ## AddFixedHalo=1 ONLY## + +# parameters for analytical function ---------- HaloUseTable=0 ONLY ----------------------------------------- +# soliton parameters +#m_22 # m_22 of the analytical halo [0.4] +#CoreRadius # core radius (in kpc) [1.0] +# alpha-beta-gamma density profile parameters +#Rho_0 # halo rho_0 (in 1.0e+10 Msun*kpc^-3) [1.0] +#Rs # halo Rs (in kpc) [1.0] +#Alpha # dimensionless [1.0] +#Beta # dimensionless [1.0] +#Gamma # dimensionless [1.0] + +# density profile table parameters ------------ HaloUseTable=1 ONLY --------------------------------------------- +DensTableFile DensTableExample # density table name, radius should be in kpc and density should be in g/cm^3 + +# parameters for new disk (disabled in this example) +AddParWhenRestart 0 # add a new disk to an existing snapshot [0] ## must set OPT__INIT=2 and OPT__RESTART_RESET=1 ## +AddParWhenRestartByFile 0 # 0=add a thin disk (must enable SUPPORT_GSL), 1=add new disk via PAR_IC [1] ## AddParWhenRestart = 1 ONLY ## +AddParWhenRestartNPar 10000000 # total particles of new disk [0] ## AddParWhenRestart = 1 ONLY ## +NewDisk_RSeed 1002 # random seed for setting thin disk particle position and velocity [1002] ## AddParWhenRestartByFile = 0 ONLY ## +Disk_Mass 0.00139812 # thin disk total mass (code unit) [1.0] ## AddParWhenRestartByFile = 0 ONLY ## +Disk_R 0.0020865 # thin disk scale radius (code unit) [1.0] ## AddParWhenRestartByFile = 0 ONLY ## +DispTableFile DispTableExample # velocity dispersion table filename, radius should be in kpc and dispersion should be in km/s ## AddParWhenRestartByFile = 0 ONLY ## diff --git a/example/test_problem/ELBDM/DiskHeating/Input__UM_IC_RefineRegion b/example/test_problem/ELBDM/DiskHeating/Input__UM_IC_RefineRegion new file mode 100644 index 0000000000..8764515e0f --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/Input__UM_IC_RefineRegion @@ -0,0 +1,9 @@ +# For loading an AMR initial condition from a file --> OPT__INIT=3 && OPT__UM_IC_NLEVEL>1 +# +# dLv : target AMR level = OPT__UM_IC_LEVEL + dLv +# NP_Skip_xL: number of patches on the parent level to be skipped in the x direction +# from the left edge of the parent refinement region +# ================================================================================== +# dLv NP_Skip_xL NP_Skip_xR NP_Skip_yL NP_Skip_yR NP_Skip_zL NP_Skip_zR + 1 7 7 7 7 7 7 + 2 25 25 25 25 25 25 diff --git a/example/test_problem/ELBDM/DiskHeating/README b/example/test_problem/ELBDM/DiskHeating/README new file mode 100644 index 0000000000..4d19908d9c --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/README @@ -0,0 +1,117 @@ +Compilation flags: +======================================== +Enable : MODEL=ELBDM, GRAVITY, PARTICLE, STORE_PAR_ACC, SUPPORT_HDF5, + SUPPORT_GSL (optional, only useful for thin disk) + +Disable : COMOVING + + +Quick start: +======================================= +0. Copy "generate_make.sh" to the directory "src", execute "sh generate_make.sh" to generate "Makefile", + and then execute "make clean" and "make -j 4" to generate the executable "gamer" + +1. Download initial conditions of m_22=0.4, M_h=7e10 Msun halo and stellar disk with "sh download_ic.sh" + +2. Ensure PAR_INIT = 1 and OPT__INIT = 3 in Input__Parameter + +3. Default END_T is 2.5e-1 (about 3.5 Gyr) as in Yang et al. 2023 and OUTPUT_DT is 1.0e-2 (about 0.14 Gyr) + +4. To switch to a high-resolution run, command "ln -sf ic_files/PAR_IC_0.4_M7 PAR_IC" + Set PAR_NPAR=80000000, MAX_LEVEL=3, and change all values in Input__Flag_NParPatch to 800 + + +General initial condition setup: +======================================== +1. Disk + + a. Generate the disk via modified GALIC (https://github.com/HsunYeong/GALIC.git) + The snapshots have the filenames snap_XXX.hdf5 + + b. Set the filename, units in get_par_ic.py to match the GALIC set-up + + c. Set center to be location of the soliton in get_par_ic.py + + d. Execute get_par_ic.py, it will generate PAR_IC + +2. Halo + + a. If the data is binary file UM_IC + + * Set OPT__INIT = 3 and PAR_INIT = 1 + * Input__UM_IC_RefineRegion is required + + b. If the data is GAMER snapshot + + * Command "ln -s Data_XXXXXX RESTART" to create a soft link + * Set OPT__INIT = 2 and PAR_INIT = 2 in Input__Parameter + * Turn on AddParWhenRestart and AddParWhenRestartByFile in Input__TestProb + * Set AddParWhenRestartNPar in Input__TestProb + * Turn on OPT__RESTART_RESET in Input__Parameter + # Recommand to turn off AddParWhenRestart, AddParWhenRestartByFile, OPT__RESTART_RESET right after the simulation starts + + c. The code for FDM halo reconstruction (https://github.com/calab-ntu/psidm-halo-reconstruction) + +3. Thin disk (optional) + + a. Command "ln -s Data_XXXXXX RESTART" to create a soft link for restart + + b. Set Set OPT__INIT = 2 and PAR_INIT = 2 and turn on OPT__RESTART_RESET in Input__Parameter + + c. Turn on AddParWhenRestart in Input__TestProb + + d. Set AddParWhenRestartNPar, Disk_Mass, Disk_R, and DispTableFile in Input__TestProb + + +Analysis scripts: +======================================== +1. particle_proj.py, plot_halo_slice.py: + + * Plot the projection of disk particles or halo density slice + * Output files: particle_proj_*.png, Data_*_Slice_x_Dens.png + +2. plot_halo_density.py, plot_halo_potential.py: + + * Compute and plot shell-averaged halo density or gravitational potential profiles + * Output files: Data_*_1d-Profile_radius_Dens.png, Halo_Dens_Data_*.npy, + Data_*_1d-Profile_radius_Pote.png, Halo_Pote_Data_*.npy + +3. data_disk.py: + + * Compute the disk information (rotation speed, velocity dispersion, surface density, scale height, etc.) + * Output files: Data_Disk_*.npy + +4. data_halo.py: + + * Compute the halo information (enclosed mass, velocity dispersion, etc.) + * Required files: Halo_Dens_*.npy (generated by "plot_halo_density.py") and Halo_Pote_*.npy (by "plot_halo_potential.py") + * Output files: Data_Halo_*.npy + +5. vel_distribution.py + + * Compute the velocity distribution in 2-kpc-wide radial bins centered on R = 4, 6, 8, 10 kpc + * Output files: Vel_data_*.npz, vel_*.png + +6. get_heating.py + + * Compute the theoretical heating rate of the stellar disk as a function of radius + * Required files: Data_Disk_*.npy (generated by "data_disk.py"), Data_Halo_*.npy (by "data_halo.py") + * Output files: Heating_*.npz + +7. get_heating_rate_theory.py + + * Get the time-averaged theoretical heating rate + * Required files: Heating_*.npz (generated by "get_heating.py") + * Output files: printed plain text + +8. get_heating_rate_simulation.py + + * Compute ensemble- and time-averaged disk heating rates measured from simulation data + * Required files: Vel_data_*.npz (generated by "vel_distribution.py") + * Output files: sigma_z_sqr.png and printed plain text + +9. plot_data_example.py + + * An example script to plot the angle-averaged disk rotation curve and shell-averaged density profile + * Required files: Data_Disk_*.npy (generated by "data_disk.py") and/or Data_Halo_*.npy (by "data_halo.py") + * Output files: Rotation_Curve.png, Halo_Density_Profile.png diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/data_disk.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/data_disk.py new file mode 100644 index 0000000000..987bca944a --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/data_disk.py @@ -0,0 +1,172 @@ +import h5py +import math +import numpy as np +import argparse +import sys + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 500 # total number of data points sampled per disk rotation curve +ratio = 0.76159415595 # tanh(1) + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Get disk properties' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def SearchIndex(x, A, N): + i = 0 + j = N - 1 + while(i <= j): + mid = int(i + (j - i)/2) + if(A[mid] == x): + i = mid + break + elif(A[mid] > x): + j = mid - 1 + else: i = mid + 1 + return i + +f = h5py.File('../../Data_%06d'%idx_start, 'r') +Unit_L = f['Info']['InputPara']['Unit_L'] +Unit_T = f['Info']['InputPara']['Unit_T']*(3.16887646e-14) +Unit_V = f['Info']['InputPara']['Unit_V'] +Unit_M = f['Info']['InputPara']['Unit_M'] +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +Center = Center * Unit_L +f.close() +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + + +#----------------------------------------------------------------------------------------------------------------- +# analyze simulation data and generate processed output files +for idx in range(idx_start, idx_end+1, didx): + print('loading file Data_%06d ...'%idx) + with h5py.File('../../Data_%06d'%idx, 'r') as f: + disk_mass = np.array(f['Particle/ParMass'], dtype = np.float64) * Unit_M + disk_posx = np.array(f['Particle/ParPosX'], dtype = np.float64) * Unit_L + disk_posy = np.array(f['Particle/ParPosY'], dtype = np.float64) * Unit_L + disk_posz = np.array(f['Particle/ParPosZ'], dtype = np.float64) * Unit_L + disk_velx = np.array(f['Particle/ParVelX'], dtype = np.float64) * Unit_V + disk_vely = np.array(f['Particle/ParVelY'], dtype = np.float64) * Unit_V + disk_velz = np.array(f['Particle/ParVelZ'], dtype = np.float64) * Unit_V + disk_type = np.array(f['Particle/ParType'], dtype = np.int32) + current_step = f['Info/KeyInfo']['Step'] + time = f['Info/KeyInfo']['Time'][0]*Unit_T + # particle filter: disk particles have ParType=2 + disk_index = (disk_type==2) + disk_mass = disk_mass[disk_index] + disk_posx = disk_posx[disk_index] + disk_posy = disk_posy[disk_index] + disk_posz = disk_posz[disk_index] + disk_velx = disk_velx[disk_index] + disk_vely = disk_vely[disk_index] + disk_velz = disk_velz[disk_index] + disk_size = np.size(disk_mass) + print("number of disk particles = %d"%disk_size) + VCM = [ np.sum(disk_mass*disk_velx)/ np.sum(disk_mass), + np.sum(disk_mass*disk_vely)/ np.sum(disk_mass), + np.sum(disk_mass*disk_velz)/ np.sum(disk_mass) ] + center = Center[current_step,3:6] + disk_posx = disk_posx - center[0] + disk_posy = disk_posy - center[1] + disk_posz = disk_posz - center[2] + disk_velx = disk_velx - VCM[0] + disk_vely = disk_vely - VCM[1] + disk_velz = disk_velz - VCM[2] + # compute angular momentum + disk_pos = np.zeros((disk_size, 3)) + disk_vel = np.zeros((disk_size, 3)) + disk_pos[:,0] = disk_posx + disk_pos[:,1] = disk_posy + disk_pos[:,2] = disk_posz + disk_vel[:,0] = disk_velx + disk_vel[:,1] = disk_vely + disk_vel[:,2] = disk_velz + disk_L = np.cross(disk_pos, disk_vel) + tot_L = np.array([np.sum(disk_L[:,0]),np.sum(disk_L[:,1]),np.sum(disk_L[:,2])]) + tot_L = tot_L/(tot_L[0]**2+tot_L[1]**2+tot_L[2]**2)**0.5 + + vec = np.cross(tot_L, np.array([0,0,1])) + c = tot_L[2] + s = 1-c**2 + matric = np.array([ [0, -vec[2], vec[1]], + [vec[2], 0, -vec[0]], + [-vec[1], vec[0], 0] ]) + R = np.identity(3) + matric + np.matmul(matric, matric)/(1+c) + disk_pos_new = np.matmul(R, np.transpose(disk_pos)) + disk_vel_new = np.matmul(R, np.transpose(disk_vel)) + disk_posx = disk_pos_new[0,:] + disk_posy = disk_pos_new[1,:] + disk_posz = disk_pos_new[2,:] + disk_velx = disk_vel_new[0,:] + disk_vely = disk_vel_new[1,:] + disk_velz = disk_vel_new[2,:] + + disk_r = (disk_posx**2 + disk_posy**2)**0.5 + disk_velr = (disk_posx*disk_velx + disk_posy*disk_vely)/disk_r + disk_velp = (disk_posx*disk_vely - disk_posy*disk_velx)/disk_r + disk_sortR = np.sort(disk_r) + disk_indexR = np.argsort(disk_r) + Data = np.zeros((8,Div_disk)) + + r = 0 + disk_num = 0 + dr = 15.0*3.08568e+21/Div_disk + for j in range(Div_disk): + disk_num_pre = disk_num + Data[0,j] = r + 0.5 * dr + disk_num = SearchIndex( r+dr, disk_sortR, disk_size ) + mean_vr = np.mean( disk_velr[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vp = np.mean( disk_velp[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vz = np.mean( disk_velz[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vp2 = np.mean( disk_velp[ disk_indexR[ disk_num_pre:disk_num ] ]**2 ) + Data[1,j] = mean_vp # rotation curve + Data[2,j] = (np.mean(( disk_velr[ disk_indexR[ disk_num_pre:disk_num ] ] - mean_vr )**2))**0.5 # sigma_r + Data[3,j] = (np.mean(( disk_velz[ disk_indexR[ disk_num_pre:disk_num ] ] - mean_vz )**2))**0.5 # sigma_z + mass = disk_mass[0]*( disk_num-disk_num_pre ) + mass_total = disk_mass[0]*disk_num + area = (math.pi *((r + dr)**2 - r**2)) + Data[4,j] = mass/area # surface density Sigma + Data[5,j] = mean_vp2 # sigma_phi^2 + Data[6,j] = mass_total # enclosed mass + # get sacle height + CMZ = np.average(disk_posz[ disk_indexR[ disk_num_pre:disk_num ] ]) + disk_z = np.abs( disk_posz[ disk_indexR[ disk_num_pre:disk_num ] ] - CMZ ) + sortZ = np.sort(disk_z) + target_index = ratio*len(disk_z) + index_plus = int(target_index) + index_minus = int(target_index - 1) + target_z = sortZ[index_minus] + (target_index - index_minus)*(sortZ[index_plus] - sortZ[index_minus]) + Data[7,j] = target_z # scale height + + r = r + dr + + np.save('Data_Disk_%06d'%idx, Data) + print(' ') diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating.py new file mode 100644 index 0000000000..c2bae1df12 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating.py @@ -0,0 +1,157 @@ +import numpy as np +import argparse +import sys +import scipy.special as sp +from scipy.integrate import quad +import yt + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 500 # total number of evenly divided radial bins +dr = 15.0/Div_disk + +r_min = 5.0e-2 +r_max = 1.2e+2 +n_bin = 500 +ratio = (r_max/r_min)**(1.0/(n_bin-1.0)) + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Compute theoretical heating rate of each snapshot' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +# ------------------------------------------------------------------------------------------------------------------------- +# specify script unit setting +''' +unit_base = {'UnitLength_in_cm' : 3.08568e+21, + 'UnitMass_in_g' : 1.989e+43, + 'UnitVelocity_in_cm_per_s' : 100000} +''' +G = 6.6743e-8 # gravitational constant (cm**3/(g*s**2)) +hb = 1.0546e-27 # reduced planck const cm**2*g/s +kappa1 = 0.526 +kappa2 = 1.0 +ILambda1 = 4.82 +ILambda2 = 3.33 + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def Gx(x): + return (sp.erf(x)-2*x*np.exp(-x*x)/(np.pi)**0.5)/2/x**2 + +def interpolation(r, array_radius, array_value): + log_radius = np.log(array_radius) + count = int(np.log(r/r_min)/np.log(ratio)) + while(r > array_radius[count]): + count = count+1 + output = array_value[count-1]+(array_value[count]-array_value[count-1])*(np.log(r)-log_radius[count-1])/(log_radius[count]-log_radius[count-1]) + return output + +def phi_bg(R, x, h): + z = x*h + radius = (R*R+z*z)**0.5/3.08568e+21 + if radius < r_p[0]: + out = pote_data[0] + elif radius > r_p[-2]: + out = pote_data[-2] + else: + out = interpolation(radius, r_p, pote_data) + return out + +def phi_int(r, h, x, x1): + phi_diff = phi_bg(r, x, h) - phi_bg(r, x1, h) + if((phi_diff) < 1.0e-50): + return 0.0 + else: + return 2.0/(2.0*phi_diff)**0.5 + +def IP2(r, x, h): + return quad(lambda x1: phi_int(r, h, x, x1), 0, x)[0] + +def IP2_log_IP2_int(r, x, h): + ip2 = IP2(r, x, h) + if (ip2 < 1.0e-50): + return 0.0 + else: + return ip2*np.log(ip2)*2.0*np.exp(-x*x)/(np.pi**0.5) + + +#----------------------------------------------------------------------------------------------------------------- +# analyze simulation data and generate processed output files +ds = yt.load( '../../Data_%06d'%idx_start ) +ma = ds.parameters["ELBDM_Mass"]*ds.parameters["Unit_M"] # FDM particle mass [gram] + +for idx in range(idx_start, idx_end+1, didx): + data = np.load('Data_Disk_%06d.npy'%idx) + data_h = np.load('../halo/Data_Halo_%06d.npy'%idx) + + # data in cgs + r = data[0] + r_in_kpc = r/3.08568e+21 + v_c = data[1] + sigma_r = data[2] + sigma_z = data[3] + Sigma = data[4] + v_c2 = data[5] + enmass = data[6] + data_h[2] + dMdR = np.gradient(enmass, r) + rho_eff = dMdR/(4.0*np.pi*r**2)-enmass/(4.0*np.pi*r**3) + sigma_h = data_h[3] + rho_h = data_h[1] + rho_bg_eff = rho_h - rho_eff + rho_bg_eff = rho_bg_eff.clip(min=0) + nu = (((v_c2 + sigma_z**2 + sigma_r**2)/3.0)**0.5)/sigma_z + r_p = data_h[4] + pote_data = data_h[5] + h = data[7] + mean_P1 = sigma_z*4.82/(np.pi*G*Sigma) + mean_P2 = np.zeros(Div_disk) + + for i in range(Div_disk): + mean_P2[i] = 2*h[i]*IP2(r[i], 1, h[i]) +# np.save('mean_P2_%06d'%idx, mean_P2) +# mean_P2 = np.load('mean_P2_%06d.npy'%idx) + + X = nu*sigma_z*3**0.5/2**0.5/sigma_h + E = (4.0*X)/(3.0*nu**2*np.pi**0.5)*np.exp(-X**2) + (1-nu**(-2))*(sp.erf(X)-Gx(X)) + tau = hb/(2**0.5*ma*sigma_h**2) + M = (np.pi**(2.5)*G**2*hb**3*rho_h**2)/(6**0.5*ma**3*sigma_h**3*nu*sigma_z)*E + + heat_SGD = kappa1*M*np.log(mean_P1/tau) + heat_BGD = kappa2*M*np.log(mean_P2/tau) + + heat = np.zeros(Div_disk) + for i in range(Div_disk): + k = sigma_z[i]*(8*G*rho_bg_eff[i])**0.5/(np.pi*G*Sigma[i]) + if (k<0.5): + heat[i] = heat_SGD[i] + elif (k>1.5): + heat[i] = heat_BGD[i] + else: + heat[i] = (k-0.5)*heat_BGD[i] + (1.5-k)*heat_SGD[i] + + np.savez('Heating_%06d'%idx, a = r, b = heat) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_simulation.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_simulation.py new file mode 100644 index 0000000000..8b3484d262 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_simulation.py @@ -0,0 +1,106 @@ +import matplotlib +matplotlib.use('Agg') +import numpy as np +import argparse +import sys +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit +from matplotlib.pyplot import cm +from matplotlib.patches import Rectangle + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Compute the time-averaged simulating heating rate' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def fit(t, a, b): + return a*t+b + + +#----------------------------------------------------------------------------------------------------------------- +# compute disk heating rates from simulation data +delta_t = 140.59 # delta_t between snapshots (Myr) +dt = delta_t*1e6*365*24*3600 +t = [] +for idx in range(idx_start, idx_end+1): + t.append(idx*delta_t) +t = np.array(t) +mean_disp = np.zeros((4, len(t))) +i = 0 +for idx in range(idx_start, idx_end+1): + ParData = np.load('Vel_data_%06d.npz'%idx) + mean_disp[0,i] = ParData['disp'][0]**2. + mean_disp[1,i] = ParData['disp'][1]**2. + mean_disp[2,i] = ParData['disp'][2]**2. + mean_disp[3,i] = ParData['disp'][3]**2. + i = i + 1 + +heating = np.zeros(4) +const = np.zeros(4) +print('R(kpc) heating rate (km^2/(s^2*Gyr))') +for i in range(4): + xData = t + yData = mean_disp[i,:] + para_F, _ = curve_fit(fit, xData, yData) + heating[i] = para_F[0] + const[i] = para_F[1] + print('%2d %2.6f'%(i*2+4, heating[i]*1000 )) + +# plot result +lines=[] +color = [cm.Blues(np.linspace(0.4, 0.7, 2)), + cm.Oranges(np.linspace(0.3, 0.6, 2)), + cm.Greens(np.linspace(0.4, 0.7, 2)), + cm.Reds(np.linspace(0.5, 0.8, 2))] + +fig = plt.figure() +axs = fig.add_subplot(111) +plt.title('$\sigma_z^2$, bin size = 2 kpc') +for i in range(4): + l1, = axs.plot(t/1000., mean_disp[i,:],'-o', ms=2, c=color[i][1]) + l2, = axs.plot(t/1000., heating[i]*t+const[i],'--', dashes=(5,1), ms=2, c=color[i][0]) + lines.append([l1, l2]) + +# create legend +extra = Rectangle((0, 0), 1, 1.0, fc="w", fill=False, edgecolor='none', linewidth=0) + +legend_handle = [extra, extra , extra , extra , extra, + extra, lines[0][0], lines[1][0], lines[2][0], lines[3][0], + extra, lines[0][1], lines[1][1], lines[2][1], lines[3][1]] + +length = len(legend_handle) +legend_handle = np.array(legend_handle).reshape(int(length/5), 5).T.reshape(length) + +label_row1 = ['' , '$R$ = 4 kpc', '$R$ = 6 kpc', '$R$ = 8 kpc', '$R$ = 10 kpc'] +label_row2 = ['simulation', '' , '' , '' , '' ] +label_row3 = ['linear fit', '' , '' , '' , '' ] +legend_labels = np.concatenate([label_row1, label_row2, label_row3]).reshape(int(length/5), 5).T.reshape(length) + +axs.legend(legend_handle, legend_labels, loc='upper left',fontsize=7, ncol = 5, handletextpad = -2.5, handlelength=2.5, handleheight=1.5) + +plt.ylim(0, 1.1*max(mean_disp[0])) +plt.xlabel(r"$t_{\rm rel}$ (Gyr)") +plt.ylabel("$\sigma_z^2$ (km$^2$/s$^2$)") +plt.savefig("sigma_z_sqr.png", dpi = 140, bbox_inches="tight") +plt.close() diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_theory.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_theory.py new file mode 100644 index 0000000000..11270b7a8f --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/get_heating_rate_theory.py @@ -0,0 +1,87 @@ +import numpy as np +import argparse +import sys + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 500 # total number of evenly divided radial bins +delta_t = 140.59 # delta_t between snapshots (Myr) + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Compute the time-averaged theoretical heating rate' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def SearchIndex(x, A, N): + i = 0 + j = N - 1 + while(i <= j): + mid = int(i + (j - i)/2) + if(A[mid] == x): + i = mid + break + elif(A[mid] > x): + j = mid - 1 + else: i = mid + 1 + return i + + +#----------------------------------------------------------------------------------------------------------------- +# compute theoretical disk heating rates +t = [] +for idx in range(idx_start, idx_end+1): + t.append(idx*delta_t) +t = np.array(t) +h_theory = np.zeros((Div_disk, len(t))) +w = np.zeros((Div_disk, len(t))) +mean_theory = np.zeros((4, len(t))) +i = 0 +for idx in range(idx_start, idx_end+1): + F = np.load('Data_Disk_%06d.npy'%idx) + H = np.load('Heating_%06d.npz'%idx) + r = F[0]/3.08568e+21 + w[:,i] = r*F[4] ## weight = r*Sigma + h_theory[:,i] = H['b']*(1e9*365*24*3600)/1e10 + i = i + 1 +dt = delta_t*1e6*365*24*3600 + +h_theory[h_theory==-np.inf]=0 + +num_3 = SearchIndex( 3, r, Div_disk ) +num_5 = SearchIndex( 5, r, Div_disk ) +num_7 = SearchIndex( 7, r, Div_disk ) +num_9 = SearchIndex( 9, r, Div_disk ) +num_11 = SearchIndex( 11, r, Div_disk ) + +for i in range (len(t)): + mean_theory[0,i] = np.average(h_theory[num_3:num_5,i], weights=w[num_3:num_5,i]) + mean_theory[1,i] = np.average(h_theory[num_5:num_7,i], weights=w[num_5:num_7,i]) + mean_theory[2,i] = np.average(h_theory[num_7:num_9,i], weights=w[num_7:num_9,i]) + mean_theory[3,i] = np.average(h_theory[num_9:num_11,i], weights=w[num_9:num_11,i]) + + +print('R(kpc) heating rate (km^2/(s^2*Gyr))') +for i in range(4): + print('%2d %2.6f'%(i*2+4, np.mean(mean_theory[i,:]) )) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/particle_proj.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/particle_proj.py new file mode 100644 index 0000000000..11642e0998 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/particle_proj.py @@ -0,0 +1,139 @@ +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import numpy as np +import yt +from mpl_toolkits.axes_grid1 import AxesGrid +import argparse +import sys + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +FONT_SIZE = 24.0 +LINE_WIDTH = 0.5 +colormap = 'arbre' + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Plot the face-on and edge-on projection of stellar disk' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +# ------------------------------------------------------------------------------------------------------------------------- +# specify script unit and output figure settings +cm = 1/2.54 # centimeters in inches +field = ('Disk', 'ParMass') +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + + +plt.rcParams['font.size'] = FONT_SIZE +plt.rcParams['figure.titlesize'] = 2*FONT_SIZE + +plt.rcParams['axes.titlesize'] = 2*FONT_SIZE +plt.rcParams['axes.labelsize'] = FONT_SIZE +plt.rcParams['axes.labelpad'] = 0.05 +plt.rcParams['axes.linewidth'] = LINE_WIDTH + +plt.rcParams['legend.fontsize'] = 6.0 +plt.rcParams['lines.linewidth'] = LINE_WIDTH + +plt.rcParams['xtick.major.size'] = 2 +plt.rcParams['xtick.major.width'] = 0.5 +plt.rcParams['xtick.minor.size'] = 1 +plt.rcParams['xtick.minor.width'] = 0.25 +plt.rcParams['xtick.labelsize'] = FONT_SIZE + +plt.rcParams['ytick.major.size'] = 2 +plt.rcParams['ytick.major.width'] = 0.5 +plt.rcParams['ytick.minor.size'] = 1 +plt.rcParams['ytick.minor.width'] = 0.25 +plt.rcParams['ytick.labelsize'] = FONT_SIZE + +plt.rcParams['font.family'] = 'STIXGeneral' +plt.rcParams['mathtext.fontset'] = 'custom' +plt.rcParams['mathtext.rm'] = 'STIXGeneral:regular' +plt.rcParams['mathtext.it'] = 'STIXGeneral:italic' +plt.rcParams['mathtext.bf'] = 'STIXGeneral:italic:bold' + +# add yt particle filter +def Disk(pfilter, data): + filter = data['all', 'ParType'] == 2 + return filter +yt.add_particle_filter('Disk', function = Disk, filtered_type = 'all', requires = ['ParType']) + +# ------------------------------------------------------------------------------------------------------------------------- +# output figures +fig = plt.figure() +fig.dpi = 150 + +grid = AxesGrid( fig, (0.1, 0.05, 3.2, 2.7), nrows_ncols=(1, 2), axes_pad=(1.2,0.5), label_mode="all", share_all=True, cbar_location="right", cbar_mode="single", cbar_size="2%", cbar_pad="2%") + +for idx in range(idx_start, idx_end+1, didx): + ds = yt.load( '../../Data_%06d'%idx ) + ds.add_particle_filter('Disk') + if sys.version_info[0] == 2: + ds.periodicity=(True,True,True) + current_step = ds.parameters["Step"] + print("Current Simulation Time = %.5e [code units]"%ds.parameters["Time"][0]) + print("Current Simulation Step = %i"%current_step) + + parx = yt.ParticleProjectionPlot( ds, 'x', fields = field, center = Center[current_step,3:6], width =( (60, 'kpc'),(60,'kpc'))) + parx.set_background_color( field ) + parx.set_cmap( field, colormap ) + parx.set_colorbar_label(field, "Projected stellar mass (M$_\odot$)") + parx.set_font( {'size':FONT_SIZE} ) + parx.set_axes_unit( 'kpc' ) + parx.set_unit( field, 'Msun' ) + parx.set_zlim( field, 5.5e+5, 5.5e+2, dynamic_range=None) + parx.annotate_text([ 0.05, 0.92], 'edge-on', coord_system="axis",text_args={"size":FONT_SIZE,"color":"black","weight":"normal","bbox":dict(boxstyle="round",ec='white',fc='white',alpha=0.7)}) + plot = parx.plots[field] + plot.figure = fig + plot.axes = grid[0].axes + plot.cax = grid.cbar_axes[0] + parx._setup_plots() + + parz = yt.ParticleProjectionPlot( ds, 'z', fields = field, center = Center[current_step,3:6], width =( (60, 'kpc'),(60,'kpc'))) + parz.set_background_color( field ) + parz.set_cmap( field, colormap ) + parz.set_colorbar_label(field, "Projected stellar mass (M$_\odot$)") + parz.set_font( {'size':FONT_SIZE} ) + parz.set_axes_unit( 'kpc' ) + parz.set_unit( field, 'Msun' ) + parz.set_zlim( field, 5.5e+5, 5.5e+2, dynamic_range=None) + parz.annotate_text([ 0.05, 0.92], 'face-on', coord_system="axis",text_args={"size":FONT_SIZE,"color":"black","weight":"normal","bbox":dict(boxstyle="round",ec='white',fc='white',alpha=0.7)}) + parz.annotate_text([ 0.75, 0.92], r'$t_{\rm rel}$ = %2.1f Gyr'%(140.59*idx/1000.), coord_system="axis", text_args={"size":FONT_SIZE,"color":"white"}) + plot = parz.plots[field] + plot.figure = fig + plot.axes = grid[1].axes + plot.cax = grid.cbar_axes[1] + parz._setup_plots() + + fig.set_size_inches(18*cm, 8*cm) + fig.savefig("particle_proj_%06d.png"%idx, bbox_inches='tight',pad_inches=0.02) + #fig.savefig("particle_proj_%06d.pdf"%idx, bbox_inches='tight',pad_inches=0.02) + + print('\nparticle_proj_%06d.png completed\n'%idx) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/vel_distribution.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/vel_distribution.py new file mode 100644 index 0000000000..7c828213dc --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/disk/vel_distribution.py @@ -0,0 +1,217 @@ +import matplotlib +matplotlib.use('Agg') +import h5py +import numpy as np +import argparse +import sys +import matplotlib.pyplot as plt + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 249 # total number of evenly divided radial bins + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Get the velocity distribution of each 2-kpc bin' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def SearchIndex(x, A, N): + i = 0 + j = N - 1 + while(i <= j): + mid = int(i + (j - i)/2) + if(A[mid] == x): + i = mid + break + elif(A[mid] > x): + j = mid - 1 + else: i = mid + 1 + return i + +def Gauss(x, mean, disp): + return np.exp(-(x-mean)**2./2./disp**2.)/disp/(2.*np.pi)**0.5 + + +#----------------------------------------------------------------------------------------------------------------- +# load simulation units and parameters +''' +unit_base = {'UnitLength_in_cm' : 3.08568e+21, + 'UnitMass_in_g' : 1.989e+43, + 'UnitVelocity_in_cm_per_s' : 100000} +''' +f = h5py.File('../../Data_%06d'%idx_start, 'r') +Unit_L = f['Info']['InputPara']['Unit_L'] +Unit_T = f['Info']['InputPara']['Unit_T']*(3.16887646e-14) +Unit_V = f['Info']['InputPara']['Unit_V'] +Unit_M = f['Info']['InputPara']['Unit_M'] +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +Center = Center * Unit_L +f.close() +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + + +#----------------------------------------------------------------------------------------------------------------- +# analyze and plot input simulation data +label = ['$R$ = 4 kpc','$R$ = 6 kpc','$R$ = 8 kpc','$R$ = 10 kpc'] +for idx in range(idx_start, idx_end+1, didx): + print('loading file Data_%06d ...'%idx) + with h5py.File('../../Data_%06d'%idx, 'r') as f: + disk_mass = np.array(f['Particle/ParMass'], dtype = np.float64) * Unit_M + disk_posx = np.array(f['Particle/ParPosX'], dtype = np.float64) * Unit_L + disk_posy = np.array(f['Particle/ParPosY'], dtype = np.float64) * Unit_L + disk_posz = np.array(f['Particle/ParPosZ'], dtype = np.float64) * Unit_L + disk_velx = np.array(f['Particle/ParVelX'], dtype = np.float64) * Unit_V + disk_vely = np.array(f['Particle/ParVelY'], dtype = np.float64) * Unit_V + disk_velz = np.array(f['Particle/ParVelZ'], dtype = np.float64) * Unit_V + disk_type = np.array(f['Particle/ParType'], dtype = np.int32) + current_step = f['Info/KeyInfo']['Step'] + time = f['Info/KeyInfo']['Time'][0]*Unit_T + # particle filter: disk particles have ParType=2 + disk_index = (disk_type==2) + disk_mass = disk_mass[disk_index] + disk_posx = disk_posx[disk_index] + disk_posy = disk_posy[disk_index] + disk_posz = disk_posz[disk_index] + disk_velx = disk_velx[disk_index] + disk_vely = disk_vely[disk_index] + disk_velz = disk_velz[disk_index] + disk_size = np.size(disk_mass) + VCM = [ np.sum(disk_mass*disk_velx)/ np.sum(disk_mass), + np.sum(disk_mass*disk_vely)/ np.sum(disk_mass), + np.sum(disk_mass*disk_velz)/ np.sum(disk_mass) ] + center = Center[current_step,3:6] + disk_posx = disk_posx - center[0] + disk_posy = disk_posy - center[1] + disk_posz = disk_posz - center[2] + disk_velx = disk_velx - VCM[0] + disk_vely = disk_vely - VCM[1] + disk_velz = disk_velz - VCM[2] + # compute angular momentum + disk_pos = np.zeros((disk_size, 3)) + disk_vel = np.zeros((disk_size, 3)) + disk_pos[:,0] = disk_posx + disk_pos[:,1] = disk_posy + disk_pos[:,2] = disk_posz + disk_vel[:,0] = disk_velx + disk_vel[:,1] = disk_vely + disk_vel[:,2] = disk_velz + disk_L = np.cross(disk_pos, disk_vel) + tot_L = np.array([np.sum(disk_L[:,0]),np.sum(disk_L[:,1]),np.sum(disk_L[:,2])]) + tot_L = tot_L/(tot_L[0]**2+tot_L[1]**2+tot_L[2]**2)**0.5 + vec = np.cross(tot_L, np.array([0,0,1])) + c = tot_L[2] + s = 1-c**2 + matric = np.array([ [0, -vec[2], vec[1]], + [vec[2], 0, -vec[0]], + [-vec[1], vec[0], 0] ]) + R = np.identity(3) + matric + np.matmul(matric, matric)/(1+c) + disk_pos_new = np.matmul(R, np.transpose(disk_pos)) + disk_vel_new = np.matmul(R, np.transpose(disk_vel)) + disk_posx = disk_pos_new[0,:] + disk_posy = disk_pos_new[1,:] + disk_posz = disk_pos_new[2,:] + disk_velx = disk_vel_new[0,:] + disk_vely = disk_vel_new[1,:] + disk_velz = disk_vel_new[2,:] + + disk_r = (disk_posx**2 + disk_posy**2)**0.5 + disk_sortR = np.sort(disk_r)/3.08568e+21 + disk_indexR = np.argsort(disk_r) + + num_3 = SearchIndex( 3, disk_sortR, disk_size ) + num_5 = SearchIndex( 5, disk_sortR, disk_size ) + num_7 = SearchIndex( 7, disk_sortR, disk_size ) + num_9 = SearchIndex( 9, disk_sortR, disk_size ) + num_11 = SearchIndex( 11, disk_sortR, disk_size ) + + vel = [] + vel.append( np.sort( disk_velz[disk_indexR[num_3:num_5]] )/1e+5 ) + vel.append( np.sort( disk_velz[disk_indexR[num_5:num_7]] )/1e+5 ) + vel.append( np.sort( disk_velz[disk_indexR[num_7:num_9]] )/1e+5 ) + vel.append( np.sort( disk_velz[disk_indexR[num_9:num_11]] )/1e+5 ) + + vmean = np.zeros(4) + disp = np.zeros(4) + mean_new = np.zeros(4) + disp_new = np.zeros(4) + vmax = np.zeros(4) + vmin = np.zeros(4) + dv = np.zeros(4) + index_max = np.zeros(4, dtype = np.int32) + index_min = np.zeros(4, dtype = np.int32) + + for i in range(4): + vmean[i] = np.mean(vel[i]) + disp[i] = np.mean( (vel[i] - vmean[i])**2. )**0.5 + vmax[i] = vmean[i] + 3.*disp[i] + vmin[i] = vmean[i] - 3.*disp[i] + index_max[i] = SearchIndex( vmax[i], vel[i], len(vel[i]) ) + index_min[i] = SearchIndex( vmin[i], vel[i], len(vel[i]) ) + mean_new[i] = np.mean(vel[i][index_min[i]:index_max[i]]) + disp_new[i] = np.mean( (vel[i][index_min[i]:index_max[i]] - mean_new[i])**2. )**0.5 + dv[i] = (vmax[i]-vmin[i])/Div_disk + + Data = np.zeros((4, 2, Div_disk)) + for i in range(4): + num_v = index_min[i] + v = vmin[i] + for j in range(Div_disk): + num_v_pre = num_v + Data[i,0,j] = v + 0.5*dv[i] + num_v = SearchIndex( v + dv[i], vel[i], len(vel[i]) ) + Data[i,1,j] = (num_v - num_v_pre)/dv[i]/(index_max[i]- index_min[i]) + v = v + dv[i] + + np.savez('Vel_data_%06d'%idx, data = Data, mean = mean_new, disp = disp_new, time_Myr = time) + + print('plotting vel_%06d.png ...'%idx) + fig, axs = plt.subplots(2, 2) + fig_num = 0 + for i in range(2): + for j in range(2): + axs[i, j].plot(Data[fig_num,0,:], Data[fig_num,1,:], color = 'blue') + axs[i, j].plot(Data[fig_num,0,:], Gauss(Data[fig_num,0,:],mean_new[fig_num],disp_new[fig_num]),'--', lw=0.8, color = 'red', label ='Gaussian') + axs[i, j].set_title(label[fig_num], fontsize=6) + axs[i, j].set_xlim((vmean[fig_num]-3*disp[fig_num],vmean[fig_num]+3*disp[fig_num])) + axs[i, j].tick_params(axis='both', labelsize=6) + axs[i, j].grid(ls='--',lw=0.3) + axs[i, j].legend(loc='best', shadow=True, fontsize=6) + + fig_num += 1 + axs[0, 0].set_ylabel('number density', fontsize=6) + axs[1, 0].set_ylabel('number density', fontsize=6) + axs[1, 0].set_xlabel('$v(km/s)$', fontsize=6 ) + axs[1, 1].set_xlabel('$v(km/s)$', fontsize=6 ) + + fig.tight_layout() + fig.suptitle('t = %4.1f Myr' %time, fontsize=8) + plt.savefig("vel_%06d.png"%idx, dpi = 140) + plt.close() + + print(' ') diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/data_halo.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/data_halo.py new file mode 100644 index 0000000000..cee42bb490 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/data_halo.py @@ -0,0 +1,204 @@ +import numpy as np +import argparse +import sys +import warnings +from scipy.integrate import quad + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 500 # total number of evenly divided radial bins +r_min = 5.0e-2 +r_max = 1.2e+2 +ratio = (r_max/r_min)**(1.0/(Div_disk-1.0)) +r_log = np.zeros(Div_disk) + +for i in range(Div_disk): + r_log[i] = r_min * (ratio)**(i) + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Convert halo properties for later heating rate calculation' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +# ------------------------------------------------------------------------------------------------------------------------- +# specify script unit setting +Unit_L = 3.08568e+21 +Unit_M = 1.989e+43 +Unit_V = 100000 +G = 6.6743e-8*Unit_M/Unit_L + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def interpolation(r, array_radius, array_value): + count = 0 + log_radius = np.log(array_radius) + while(r > array_radius[count]): + count = count+1 + output = array_value[count-1]+(array_value[count]-array_value[count-1])*(np.log(r)-log_radius[count-1])/(log_radius[count]-log_radius[count-1]) + return output + +def log_interpolation(r, array_radius, array_value): + count = 0 + log_value = np.log(array_value) + log_radius = np.log(array_radius) + while(r > array_radius[count]): + count = count+1 + output = log_value[count-1]+(log_value[count]-log_value[count-1])*(np.log(r)-log_radius[count-1])/(log_radius[count]-log_radius[count-1]) + return np.exp(output) +def density(r): + if (isinstance(r, int) or isinstance(r, float)): + if r < r_d[0]: + return dens[0] + elif r > r_d[-1]: + return 0 + else: + return log_interpolation(r, r_d, dens) + else: + out = np.zeros(len(r)) + for i in range(len(r)): + if r[i] < r_d[0]: + out[i] = dens[0] + elif r[i] > r_d[-1]: + out[i] = dens[-1] + else: + out[i] = log_interpolation(r[i], r_d, dens) + return out + +def enclosed_mass_int(r): + warnings.filterwarnings("ignore") + if (isinstance(r, int) or isinstance(r, float)): + return quad(lambda x: 4*np.pi*x*x*density(x), 0, r)[0] + else: + temp = [] + for i in r: + temp.append( quad(lambda x: 4*np.pi*x*x*density(x), 0, i)[0] ) + result = np.array(temp) + return result + +def potential(r): + if (isinstance(r, int) or isinstance(r, float)): + if r < r_p[0]: + return pote[0] + elif r > r_p[-1]: + return pote[-1] + else: + return log_interpolation(r, r_p, pote) + else: + out = np.zeros(len(r)) + for i in range(len(r)): + if r[i] < r_p[0]: + out[i] = pote[0] + elif r[i] > r_p[-1]: + out[i] = pote[-1] + else: + out[i] = interpolation(r[i], r_p, pote) + return out + +def jeans_int(r): + if (isinstance(r, int) or isinstance(r, float)): + return quad(lambda x: density(x)*jeans_interpolation(x), r, np.inf)[0] + else: + temp = [] + for i in r: + temp.append( quad(lambda x: density(x)*jeans_interpolation(x), i, np.inf)[0] ) + result = np.array(temp) + return result + +def jeans_interpolation(r): + if (isinstance(r, int) or isinstance(r, float)): + if r < r_p[0]: + return DPotDr[0] + elif r > r_p[-1]: + return DPotDr[-1] + else: + return interpolation(r, r_p, DPotDr) + else: + out = np.zeros(len(r)) + for i in range(len(r)): + if r[i] < r_p[0]: + out[i] = DPotDr[0] + elif r[i] > r_p[-1]: + out[i] = DPotDr[-1] + else: + out[i] = interpolation(r[i], r_p, DPotDr) + return out + + +#----------------------------------------------------------------------------------------------------------------- +# analyze simulation data and generate processed output files +for idx in range(idx_start, idx_end+1, didx): + Dens = np.load('Halo_Dens_Data_%06d.npy'%idx) + r_d = Dens[0] + dens_cgs = Dens[1] + dens = dens_cgs*Unit_L**3/Unit_M + enmass = enclosed_mass_int(r_log) + mass_cgs = enmass*Unit_M + + Pote = np.load('Halo_Pote_Data_%06d.npy'%idx) + r_p = Pote[0] + pote = Pote[1] + pote_cgs = potential(r_log) + + DPotDr = np.gradient(pote, r_p) + sigma = ((jeans_int(r_p)/density(r_p))*0.5)**0.5 + + r_in_kpc = np.zeros(Div_disk) + rho = np.zeros(Div_disk) + M_h = np.zeros(Div_disk) + disp = np.zeros(Div_disk) + phi_h = np.zeros(Div_disk) + + r = 0 + num_dens = 0 + num_log = 0 + num_disp = 0 + + dr = 15.0/float(Div_disk) + for j in range(Div_disk): + r_in_kpc[j] = r + 0.5 * dr + while r_d[num_dens] < (r + 0.5*dr): + num_dens += 1 + rho[j] = dens_cgs[num_dens-1] + (dens_cgs[num_dens] - dens_cgs[num_dens-1])*(r_in_kpc[j]-r_d[num_dens-1])/(r_d[num_dens]-r_d[num_dens-1]) + while r_log[num_log] < (r + 0.5*dr): + num_log += 1 + M_h[j] = mass_cgs[num_log-1] + (mass_cgs[num_log] - mass_cgs[num_log-1])*(r_in_kpc[j]-r_log[num_log-1])/(r_log[num_log]-r_log[num_log-1]) + while r_p[num_disp] < (r + 0.5*dr): + num_disp += 1 + disp[j] = sigma[num_disp-1] + (sigma[num_disp] - sigma[num_disp-1])*(r_in_kpc[j]-r_p[num_disp-1])/(r_p[num_disp]-r_p[num_disp-1]) + + r = r + dr + + Data = np.zeros((6,Div_disk)) + Data[0] = r_in_kpc + Data[1] = rho # halo density + Data[2] = M_h # enclosed mass + Data[3] = disp # Jeans velocity dispersion/sqrt(2) + Data[4] = r_log # r in log bins (used for computation) + Data[5] = pote_cgs # Phi(halo potential) in log bins (used for computation) + + np.save('Data_Halo_%06d'%idx, Data) + print('Data_Halo_%06d.npy completed'%idx) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_density.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_density.py new file mode 100644 index 0000000000..314105a8ec --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_density.py @@ -0,0 +1,70 @@ +import argparse +import sys +import yt +import numpy as np + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +dpi = 150 +nbin = 256*2 + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Get the halo density profile' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])), +print( '' ) +print( '-------------------------------------------------------------------\n' ) + +# ------------------------------------------------------------------------------------------------------------------------- +# output figures +field = ('gamer', 'Dens') +center_mode = 'max' +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + +for idx in range(idx_start, idx_end+1, didx): + ds = yt.load( '../../Data_%06d'%idx ) + current_step = ds.parameters["Step"] + CellSize = ds.parameters["CellSize"] + MaxLevel = ds.parameters["MaxLevel"] + Resolution = CellSize[MaxLevel]*ds.parameters["Unit_L"]/3.08568e+21 # spatial resolution in kpc + print("Current Simulation Time = %.5e [code units]"%ds.parameters["Time"][0]) + print("Current Simulation Step = %i"%current_step) + + sp = ds.sphere( Center[current_step,3:6], 0.5*ds.domain_width.to_value().max() ) + prof = yt.create_profile( sp, 'radius', field, + weight_field='cell_volume', + n_bins=nbin, + units={"radius":"kpc", field:"g/cm**3"}, + logs={"radius":True, field:True}, + extrema={"radius":((Resolution/2., "kpc"), (0.5*ds.domain_width.to_value().max(), "code_length"))}) + plot = yt.ProfilePlot.from_profiles(prof) + plot.save( mpl_kwargs={"dpi":dpi} ) + + A = prof.x + B = prof[field] + C = [A,B] + Data = np.asarray(C) + Data = np.delete(Data, np.nonzero(Data[1]==0)[0] , axis = 1) + np.save("Halo_Dens_"+str(ds),Data) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_potential.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_potential.py new file mode 100644 index 0000000000..684e9c2608 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_potential.py @@ -0,0 +1,70 @@ +import argparse +import sys +import yt +import numpy as np + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +dpi = 150 +nbin = 256*2 + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Get the halo potential profile' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])), +print( '' ) +print( '-------------------------------------------------------------------\n' ) + +# ------------------------------------------------------------------------------------------------------------------------- +# output figures +field = ('gamer', 'Pote') +center_mode = 'max' +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + +for idx in range(idx_start, idx_end+1, didx): + ds = yt.load( '../../Data_%06d'%idx ) + current_step = ds.parameters["Step"] + CellSize = ds.parameters["CellSize"] + MaxLevel = ds.parameters["MaxLevel"] + Resolution = CellSize[MaxLevel]*ds.parameters["Unit_L"]/3.08568e+21 # spatial resolution in kpc + print("Current Simulation Time = %.5e [code units]"%ds.parameters["Time"][0]) + print("Current Simulation Step = %i"%current_step) + + sp = ds.sphere( Center[current_step,3:6], 0.5*ds.domain_width.to_value().max() ) + prof = yt.create_profile( sp, 'radius', field, + weight_field='cell_volume', + n_bins=nbin, + units={"radius":"kpc", field:"cm**2/s**2"}, + logs={"radius":True, field:False}, + extrema={"radius":((Resolution/2., "kpc"), (0.5*ds.domain_width.to_value().max(), "code_length"))}) + plot = yt.ProfilePlot.from_profiles(prof) + plot.set_log(field, False) + plot.save( mpl_kwargs={"dpi":dpi} ) + A = prof.x + B = prof[field] + C = [A,B] + Data = np.asarray(C) + Data = np.delete(Data, np.nonzero(Data[1]==0)[0] , axis = 1) + np.save("Halo_Pote_"+str(ds),Data) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_slice.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_slice.py new file mode 100644 index 0000000000..6273f14820 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/halo/plot_halo_slice.py @@ -0,0 +1,66 @@ +import argparse +import sys +import yt +import numpy as np + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +colormap = 'algae' +dpi = 150 + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Plot the halo slices' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +# ------------------------------------------------------------------------------------------------------------------------- +# output figures +yt.enable_parallelism() + +field = ('Dens') +center_mode = 'max' +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + +for idx in range(idx_start, idx_end+1, didx): + ds = yt.load( '../../Data_%06d'%idx ) + if sys.version_info[0] == 2: + ds.periodicity = (True, True, True) + current_step = ds.parameters["Step"] + print("Current Simulation Time = %.5e [code units]"%ds.parameters["Time"][0]) + print("Current Simulation Step = %i"%current_step) + + dens = yt.SlicePlot( ds, 0, fields = field, center = Center[current_step,3:6], width = (60 ,'kpc')) + dens.set_background_color( field ) + dens.set_zlim( field, 1.0e+1, 1.0e+9, dynamic_range=None) + dens.set_cmap( field, colormap ) + dens.set_font( {'size':16} ) + dens.set_axes_unit( 'kpc' ) + dens.annotate_grids() + dens.annotate_sphere( center = Center[current_step,3:6], radius=(0.05, 'kpc'),circle_args={'color':'red'}) + dens.annotate_timestamp( time_unit='Myr', corner='upper_right', text_args={'color':'k'} ) + dens.save( mpl_kwargs={"dpi":dpi} ) diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/plot_data_example.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/plot_data_example.py new file mode 100644 index 0000000000..cdaa9f00a0 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/plot_data_example.py @@ -0,0 +1,87 @@ +import matplotlib +matplotlib.use('Agg') +import numpy as np +import argparse +import sys +import matplotlib.pyplot as plt +from matplotlib.pyplot import cm + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +output_mode = int(3) # [1] angle-averaged rotation curve; [2] shell-averaged halo density profile; [3] both +figure_dpi = 140 + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='An example script to plot disk/halo data from existing .npy files' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#------------------------------------------------------------------------------------------------------------------------- +# analyze and plot input simulation data +# [1] angle-averaged rotation curve +if (output_mode == 1 or output_mode == 3): + plt.figure(dpi = figure_dpi) + for idx in range(idx_start, idx_end+1, didx): + + data = np.load('disk/Data_Disk_%06d.npy'%idx) + r = data[0] + r_in_kpc = r/3.08568e+21 + v_c = data[1]/1.0e5 # rotaion speed + sigma_r = data[2]/1.0e5 # sigma_r + sigma_z = data[3]/1.0e5 # sigma_z + Sigma = data[4] # surface density + height = data[7] # scale height + plt.plot(r_in_kpc, v_c, color = cm.Blues(0.3+0.6*(idx-idx_start)/(idx_end+1-idx_start)), label ='id=%d'%idx) + + plt.xlim((0, 15)) + plt.ylim((0, 100)) + plt.grid(ls='--') + plt.legend(loc='lower right', shadow=True, prop={'size':8}) ## loc='best', 'upper left', 'upper right', 'lower left', 'lower right' + plt.xlabel(r"$R$ (kpc)") + plt.ylabel(r"$v_{\rm cir}$ (km/s)") + plt.savefig("Rotation_Curve.png", dpi = figure_dpi) + plt.close() + +# [2] plot halo density profile +if (output_mode == 2 or output_mode == 3): + plt.figure(dpi = figure_dpi) + for idx in range(idx_start, idx_end+1, didx): + + Dens = np.load('halo/Halo_Dens_Data_%06d.npy'%idx) + Dens = np.delete(Dens, np.nonzero(Dens[1]==0)[0] , axis = 1) + r_d = Dens[0] + dens_cgs = Dens[1] + plt.plot(r_d, dens_cgs, color = cm.Blues(0.3+0.6*(idx-idx_start)/(idx_end+1-idx_start)), label ='id=%d'%idx) + + plt.xlim((1e-1, 1e2)) + plt.grid(ls='--') + plt.legend(loc='lower left', shadow=True, prop={'size':8}) ## loc='best', 'upper left', 'upper right', 'lower left', 'lower right' + plt.xlabel(r"$r$ (kpc)") + plt.ylabel(r"$\rho_{\rm h}$ (g/cm$^3$)") + plt.xscale('log') + plt.yscale('log') + plt.savefig("Halo_Density_Profile.png", dpi = figure_dpi) + plt.close() diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/data_thin_disk.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/data_thin_disk.py new file mode 100644 index 0000000000..bcebd1751a --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/data_thin_disk.py @@ -0,0 +1,172 @@ +import h5py +import math +import numpy as np +import argparse +import sys + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +Div_disk = 500 # total number of data points sampled per disk rotation curve +ratio = 0.76159415595 # tanh(1) + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Get disk properties' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +#----------------------------------------------------------------------------------------------------------------- +# predefined functions +def SearchIndex(x, A, N): + i = 0 + j = N - 1 + while(i <= j): + mid = int(i + (j - i)/2) + if(A[mid] == x): + i = mid + break + elif(A[mid] > x): + j = mid - 1 + else: i = mid + 1 + return i + +f = h5py.File('../../Data_%06d'%idx_start, 'r') +Unit_L = f['Info']['InputPara']['Unit_L'] +Unit_T = f['Info']['InputPara']['Unit_T']*(3.16887646e-14) +Unit_V = f['Info']['InputPara']['Unit_V'] +Unit_M = f['Info']['InputPara']['Unit_M'] +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +Center = Center * Unit_L +f.close() +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + + +#----------------------------------------------------------------------------------------------------------------- +# analyze simulation data and generate processed output files +for idx in range(idx_start, idx_end+1, didx): + print('loading file Data_%06d ...'%idx) + with h5py.File('../../Data_%06d'%idx, 'r') as f: + disk_mass = np.array(f['Particle/ParMass'], dtype = np.float64) * Unit_M + disk_posx = np.array(f['Particle/ParPosX'], dtype = np.float64) * Unit_L + disk_posy = np.array(f['Particle/ParPosY'], dtype = np.float64) * Unit_L + disk_posz = np.array(f['Particle/ParPosZ'], dtype = np.float64) * Unit_L + disk_velx = np.array(f['Particle/ParVelX'], dtype = np.float64) * Unit_V + disk_vely = np.array(f['Particle/ParVelY'], dtype = np.float64) * Unit_V + disk_velz = np.array(f['Particle/ParVelZ'], dtype = np.float64) * Unit_V + disk_type = np.array(f['Particle/ParType'], dtype = np.int32) + current_step = f['Info/KeyInfo']['Step'] + time = f['Info/KeyInfo']['Time'][0]*Unit_T + # particle filter: thin disk particles have ParType=3 + disk_index = (disk_type==3) + disk_mass = disk_mass[disk_index] + disk_posx = disk_posx[disk_index] + disk_posy = disk_posy[disk_index] + disk_posz = disk_posz[disk_index] + disk_velx = disk_velx[disk_index] + disk_vely = disk_vely[disk_index] + disk_velz = disk_velz[disk_index] + disk_size = np.size(disk_mass) + print("number of thin disk particles = %d"%disk_size) + VCM = [ np.sum(disk_mass*disk_velx)/ np.sum(disk_mass), + np.sum(disk_mass*disk_vely)/ np.sum(disk_mass), + np.sum(disk_mass*disk_velz)/ np.sum(disk_mass) ] + center = Center[current_step,3:6] + disk_posx = disk_posx - center[0] + disk_posy = disk_posy - center[1] + disk_posz = disk_posz - center[2] + disk_velx = disk_velx - VCM[0] + disk_vely = disk_vely - VCM[1] + disk_velz = disk_velz - VCM[2] + # compute angular momentum + disk_pos = np.zeros((disk_size, 3)) + disk_vel = np.zeros((disk_size, 3)) + disk_pos[:,0] = disk_posx + disk_pos[:,1] = disk_posy + disk_pos[:,2] = disk_posz + disk_vel[:,0] = disk_velx + disk_vel[:,1] = disk_vely + disk_vel[:,2] = disk_velz + disk_L = np.cross(disk_pos, disk_vel) + tot_L = np.array([np.sum(disk_L[:,0]),np.sum(disk_L[:,1]),np.sum(disk_L[:,2])]) + tot_L = tot_L/(tot_L[0]**2+tot_L[1]**2+tot_L[2]**2)**0.5 + + vec = np.cross(tot_L, np.array([0,0,1])) + c = tot_L[2] + s = 1-c**2 + matric = np.array([ [0, -vec[2], vec[1]], + [vec[2], 0, -vec[0]], + [-vec[1], vec[0], 0] ]) + R = np.identity(3) + matric + np.matmul(matric, matric)/(1+c) + disk_pos_new = np.matmul(R, np.transpose(disk_pos)) + disk_vel_new = np.matmul(R, np.transpose(disk_vel)) + disk_posx = disk_pos_new[0,:] + disk_posy = disk_pos_new[1,:] + disk_posz = disk_pos_new[2,:] + disk_velx = disk_vel_new[0,:] + disk_vely = disk_vel_new[1,:] + disk_velz = disk_vel_new[2,:] + + disk_r = (disk_posx**2 + disk_posy**2)**0.5 + disk_velr = (disk_posx*disk_velx + disk_posy*disk_vely)/disk_r + disk_velp = (disk_posx*disk_vely - disk_posy*disk_velx)/disk_r + disk_sortR = np.sort(disk_r) + disk_indexR = np.argsort(disk_r) + Data = np.zeros((8,Div_disk)) + + r = 0 + disk_num = 0 + dr = 15.0*3.08568e+21/Div_disk + for j in range(Div_disk): + disk_num_pre = disk_num + Data[0,j] = r + 0.5 * dr + disk_num = SearchIndex( r+dr, disk_sortR, disk_size ) + mean_vr = np.mean( disk_velr[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vp = np.mean( disk_velp[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vz = np.mean( disk_velz[ disk_indexR[ disk_num_pre:disk_num ] ] ) + mean_vp2 = np.mean( disk_velp[ disk_indexR[ disk_num_pre:disk_num ] ]**2 ) + Data[1,j] = mean_vp # rotation curve + Data[2,j] = (np.mean(( disk_velr[ disk_indexR[ disk_num_pre:disk_num ] ] - mean_vr )**2))**0.5 # sigma_r + Data[3,j] = (np.mean(( disk_velz[ disk_indexR[ disk_num_pre:disk_num ] ] - mean_vz )**2))**0.5 # sigma_z + mass = disk_mass[0]*( disk_num-disk_num_pre ) + mass_total = disk_mass[0]*disk_num + area = (math.pi *((r + dr)**2 - r**2)) + Data[4,j] = mass/area # surface density Sigma + Data[5,j] = mean_vp2 # sigma_phi^2 + Data[6,j] = mass_total # enclosed mass + # get sacle height + CMZ = np.average(disk_posz[ disk_indexR[ disk_num_pre:disk_num ] ]) + disk_z = np.abs( disk_posz[ disk_indexR[ disk_num_pre:disk_num ] ] - CMZ ) + sortZ = np.sort(disk_z) + target_index = ratio*len(disk_z) + index_plus = int(target_index) + index_minus = int(target_index - 1) + target_z = sortZ[index_minus] + (target_index - index_minus)*(sortZ[index_plus] - sortZ[index_minus]) + Data[7,j] = target_z # scale height + + r = r + dr + + np.save('Data_Thin_Disk_%06d'%idx, Data) + print(' ') diff --git a/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/particle_proj_thin_disk.py b/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/particle_proj_thin_disk.py new file mode 100644 index 0000000000..fe3251f135 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/analysis_script/thin_disk/particle_proj_thin_disk.py @@ -0,0 +1,139 @@ +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import numpy as np +import yt +from mpl_toolkits.axes_grid1 import AxesGrid +import argparse +import sys + + +# ------------------------------------------------------------------------------------------------------------------------- +# user-specified parameters +FONT_SIZE = 24.0 +LINE_WIDTH = 0.5 +colormap = 'arbre' + + +#------------------------------------------------------------------------------------------------------------------------- +# load the command-line parameters +parser = argparse.ArgumentParser( description='Plot the face-on and edge-on projection of stellar disk' ) + +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx + +# print command-line parameters +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +for t in range( len(sys.argv) ): + print( str(sys.argv[t])) +print( '' ) +print( '-------------------------------------------------------------------\n' ) + + +# ------------------------------------------------------------------------------------------------------------------------- +# specify script unit and output figure settings +cm = 1/2.54 # centimeters in inches +field = ('ThinDisk', 'ParMass') +Center = np.loadtxt('../../Record__Center', skiprows=1, dtype=float) +if Center.ndim == 1: + Center = Center.reshape(1,len(Center)) # reshape the array if there is only one row + + +plt.rcParams['font.size'] = FONT_SIZE +plt.rcParams['figure.titlesize'] = 2*FONT_SIZE + +plt.rcParams['axes.titlesize'] = 2*FONT_SIZE +plt.rcParams['axes.labelsize'] = FONT_SIZE +plt.rcParams['axes.labelpad'] = 0.05 +plt.rcParams['axes.linewidth'] = LINE_WIDTH + +plt.rcParams['legend.fontsize'] = 6.0 +plt.rcParams['lines.linewidth'] = LINE_WIDTH + +plt.rcParams['xtick.major.size'] = 2 +plt.rcParams['xtick.major.width'] = 0.5 +plt.rcParams['xtick.minor.size'] = 1 +plt.rcParams['xtick.minor.width'] = 0.25 +plt.rcParams['xtick.labelsize'] = FONT_SIZE + +plt.rcParams['ytick.major.size'] = 2 +plt.rcParams['ytick.major.width'] = 0.5 +plt.rcParams['ytick.minor.size'] = 1 +plt.rcParams['ytick.minor.width'] = 0.25 +plt.rcParams['ytick.labelsize'] = FONT_SIZE + +plt.rcParams['font.family'] = 'STIXGeneral' +plt.rcParams['mathtext.fontset'] = 'custom' +plt.rcParams['mathtext.rm'] = 'STIXGeneral:regular' +plt.rcParams['mathtext.it'] = 'STIXGeneral:italic' +plt.rcParams['mathtext.bf'] = 'STIXGeneral:italic:bold' + +# add yt particle filter +def Disk(pfilter, data): + filter = data['all', 'ParType'] == 3 + return filter +yt.add_particle_filter('ThinDisk', function = Disk, filtered_type = 'all', requires = ['ParType']) + +# ------------------------------------------------------------------------------------------------------------------------- +# output figures +fig = plt.figure() +fig.dpi = 150 + +grid = AxesGrid( fig, (0.1, 0.05, 3.2, 2.7), nrows_ncols=(1, 2), axes_pad=(1.2,0.5), label_mode="all", share_all=True, cbar_location="right", cbar_mode="single", cbar_size="2%", cbar_pad="2%") + +for idx in range(idx_start, idx_end+1, didx): + ds = yt.load( '../../Data_%06d'%idx ) + ds.add_particle_filter('ThinDisk') + if sys.version_info[0] == 2: + ds.periodicity=(True,True,True) + current_step = ds.parameters["Step"] + print("Current Simulation Time = %.5e [code units]"%ds.parameters["Time"][0]) + print("Current Simulation Step = %i"%current_step) + + parx = yt.ParticleProjectionPlot( ds, 'x', fields = field, center = Center[current_step,3:6], width =( (60, 'kpc'),(60,'kpc'))) + parx.set_background_color( field ) + parx.set_cmap( field, colormap ) + parx.set_colorbar_label(field, "Projected stellar mass (M$_\odot$)") + parx.set_font( {'size':FONT_SIZE} ) + parx.set_axes_unit( 'kpc' ) + parx.set_unit( field, 'Msun' ) + parx.set_zlim( field, 5.5e+5, 5.5e+2, dynamic_range=None) + parx.annotate_text([ 0.05, 0.92], 'edge-on', coord_system="axis",text_args={"size":FONT_SIZE,"color":"black","weight":"normal","bbox":dict(boxstyle="round",ec='white',fc='white',alpha=0.7)}) + plot = parx.plots[field] + plot.figure = fig + plot.axes = grid[0].axes + plot.cax = grid.cbar_axes[0] + parx._setup_plots() + + parz = yt.ParticleProjectionPlot( ds, 'z', fields = field, center = Center[current_step,3:6], width =( (60, 'kpc'),(60,'kpc'))) + parz.set_background_color( field ) + parz.set_cmap( field, colormap ) + parz.set_colorbar_label(field, "Projected stellar mass (M$_\odot$)") + parz.set_font( {'size':FONT_SIZE} ) + parz.set_axes_unit( 'kpc' ) + parz.set_unit( field, 'Msun' ) + parz.set_zlim( field, 5.5e+5, 5.5e+2, dynamic_range=None) + parz.annotate_text([ 0.05, 0.92], 'face-on', coord_system="axis",text_args={"size":FONT_SIZE,"color":"black","weight":"normal","bbox":dict(boxstyle="round",ec='white',fc='white',alpha=0.7)}) + parz.annotate_text([ 0.75, 0.92], r'$t_{\rm rel}$ = %2.1f Gyr'%(140.59*idx/1000.), coord_system="axis", text_args={"size":FONT_SIZE,"color":"white"}) + plot = parz.plots[field] + plot.figure = fig + plot.axes = grid[1].axes + plot.cax = grid.cbar_axes[1] + parz._setup_plots() + + fig.set_size_inches(18*cm, 8*cm) + fig.savefig("particle_proj_%06d.png"%idx, bbox_inches='tight',pad_inches=0.02) + #fig.savefig("particle_proj_%06d.pdf"%idx, bbox_inches='tight',pad_inches=0.02) + + print('\nparticle_proj_%06d.png completed\n'%idx) diff --git a/example/test_problem/ELBDM/DiskHeating/clean.sh b/example/test_problem/ELBDM/DiskHeating/clean.sh new file mode 100644 index 0000000000..dd2cc95324 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/clean.sh @@ -0,0 +1,6 @@ +rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Dump Record__MemInfo Record__L1Err \ + Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \ + Diag* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \ + PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \ + Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance Record__Center \ + GRACKLE_INFO Record__DivB Record__Hybrid diff --git a/example/test_problem/ELBDM/DiskHeating/download_ic.sh b/example/test_problem/ELBDM/DiskHeating/download_ic.sh new file mode 100644 index 0000000000..34df8e4e4d --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/download_ic.sh @@ -0,0 +1,9 @@ +filename=disk-heating-ic + +curl https://girder.hub.yt/api/v1/item/6645cffcff473673ea91b24d/download -o ${filename}.tgz +tar -zxvf ${filename}.tgz +rm ${filename}.tgz +ln -s ${filename}/UM_IC_0.4_M7 UM_IC +ln -s ${filename}/PAR_IC_0.4_M7_low_res PAR_IC + + diff --git a/example/test_problem/ELBDM/DiskHeating/generate_make.sh b/example/test_problem/ELBDM/DiskHeating/generate_make.sh new file mode 100644 index 0000000000..47257f54a3 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/generate_make.sh @@ -0,0 +1,7 @@ +# This script should run in the same directory as configure.py + +PYTHON=python3 + +${PYTHON} configure.py --machine=eureka_intel --mpi=true --hdf5=true --fftw=FFTW3 --gpu=true \ + --model=ELBDM --gravity=true --particle=true --store_par_acc=true \ + --gsl=true --max_patch=20000000 diff --git a/example/test_problem/ELBDM/DiskHeating/get_par_ic.py b/example/test_problem/ELBDM/DiskHeating/get_par_ic.py new file mode 100644 index 0000000000..471d8ce186 --- /dev/null +++ b/example/test_problem/ELBDM/DiskHeating/get_par_ic.py @@ -0,0 +1,85 @@ +import h5py +import numpy as np + +filename = "snap_001.hdf5" +UNIT_L_GALIC = 3.08568e+21 +UNIT_M_GALIC = 1.989e+43 +UNIT_D_GALIC = UNIT_M_GALIC/UNIT_L_GALIC**3 +UNIT_V_GALIC = 1.0e+5 + +UNIT_L = 4.436632034507548e+24 +UNIT_D = 2.5758579703106658e-30 +UNIT_M = UNIT_D*UNIT_L**3 +UNIT_V = 1.0e+7 + +def SearchIndex(x, A, N): + i = 0 + j = N - 1 + while(i <= j): + mid = int(i + (j - i)/2) + if(A[mid] == x): + i = mid + break + elif(A[mid] > x): + j = mid - 1 + else: i = mid + 1 + return i + +with h5py.File(filename, "r") as f: + N_Disk = f['PartType2']['Masses'].size + mass = f['PartType2']['Masses'][0] *UNIT_M_GALIC/UNIT_M + disk_posx = np.array(f['PartType2']['Coordinates'][:,0], dtype = np.float64) *UNIT_L_GALIC/UNIT_L + disk_posy = np.array(f['PartType2']['Coordinates'][:,1], dtype = np.float64) *UNIT_L_GALIC/UNIT_L + disk_posz = np.array(f['PartType2']['Coordinates'][:,2], dtype = np.float64) *UNIT_L_GALIC/UNIT_L + disk_velx = np.array(f['PartType2']['Velocities'][:,0], dtype = np.float64) *UNIT_V_GALIC/UNIT_V + disk_vely = np.array(f['PartType2']['Velocities'][:,1], dtype = np.float64) *UNIT_V_GALIC/UNIT_V + disk_velz = np.array(f['PartType2']['Velocities'][:,2], dtype = np.float64) *UNIT_V_GALIC/UNIT_V + +CM = [ np.sum(disk_posx)/ N_Disk, + np.sum(disk_posy)/ N_Disk, + np.sum(disk_posz)/ N_Disk ] +VCM = [ np.sum(disk_velx)/ N_Disk, + np.sum(disk_vely)/ N_Disk, + np.sum(disk_velz)/ N_Disk ] +#print(np.array(VCM)*UNIT_V/UNIT_V_GALIC) + +center = [7.9770907e-02, 7.9854590e-02, 8.0858787e-02] +disk_posx = disk_posx - CM[0] +disk_posy = disk_posy - CM[1] +disk_posz = disk_posz - CM[2] +disk_velx = disk_velx - VCM[0] +disk_vely = disk_vely - VCM[1] +disk_velz = disk_velz - VCM[2] +disk_r = (disk_posx**2+disk_posy**2+disk_posz**2)**0.5 +sortR_kpc = np.sort(disk_r)*UNIT_L/UNIT_L_GALIC + +#print(sortR_kpc[0]) +#print(sortR_kpc[-1]) + +indexR = np.argsort(disk_r) + +# exclude particles with r > 70 kpc +num = SearchIndex( 70, sortR_kpc, N_Disk ) + +disk_posx = disk_posx[ indexR[:num] ] + center[0] +disk_posy = disk_posy[ indexR[:num] ] + center[1] +disk_posz = disk_posz[ indexR[:num] ] + center[2] +disk_velx = disk_velx[ indexR[:num] ] +disk_vely = disk_vely[ indexR[:num] ] +disk_velz = disk_velz[ indexR[:num] ] + +disk_type = np.full(num, 2) +disk_mass = np.full(num, mass) +with open('PAR_IC', 'wb') as output: + output.write(disk_mass.astype('f').tobytes()) + output.write(disk_posx.astype('f').tobytes()) + output.write(disk_posy.astype('f').tobytes()) + output.write(disk_posz.astype('f').tobytes()) + output.write(disk_velx.astype('f').tobytes()) + output.write(disk_vely.astype('f').tobytes()) + output.write(disk_velz.astype('f').tobytes()) + output.write(disk_type.astype('f').tobytes()) +output.close() +#print('Disk particles = '+ str(len(disk_posx))) +print('Disk particles = '+ str(num)) +print('PAR_IC complete') diff --git a/example/test_problem/ELBDM/LSS_Hybrid/download_heavy_halo_ic.sh b/example/test_problem/ELBDM/LSS_Hybrid/download_heavy_halo_ic.sh index 3262c7e4eb..8d03ee68cf 100644 --- a/example/test_problem/ELBDM/LSS_Hybrid/download_heavy_halo_ic.sh +++ b/example/test_problem/ELBDM/LSS_Hybrid/download_heavy_halo_ic.sh @@ -3,5 +3,5 @@ link=https://use.yt/upload/4587e2e6 curl -L ${link} -o ${filename} ln -sf ${filename} UM_IC_wave_heavy -python3 elbdm_wave_to_hybrid.py -input UM_IC_wave_heavy -output UM_IC_hybrid_heavy -resolution 256 +python3 elbdm_wave_to_hybrid_IC.py -input UM_IC_wave_heavy -output UM_IC_hybrid_heavy -resolution 256 ln -sf UM_IC_hybrid_heavy UM_IC diff --git a/include/Macro.h b/include/Macro.h index 76a29ac415..408b5034a1 100644 --- a/include/Macro.h +++ b/include/Macro.h @@ -917,7 +917,7 @@ # define GRAMFE_FFT_FLOAT8 #endif -#if ( GRAMFE_SCHEME == GRAMFE_MATMUL ) +#if ( ( GRAMFE_SCHEME == GRAMFE_MATMUL ) && defined( FLOAT8 ) ) # define GRAMFE_MATMUL_FLOAT8 #endif diff --git a/include/Typedef.h b/include/Typedef.h index 6967585b2a..2ca69553e6 100644 --- a/include/Typedef.h +++ b/include/Typedef.h @@ -92,8 +92,8 @@ const TestProbID_t TESTPROB_ELBDM_LSS = 1009, TESTPROB_ELBDM_PLANE_WAVE = 1010, TESTPROB_ELBDM_PERTURBATION = 1011, - TESTPROB_ELBDM_HALO_MERGER = 1012; - + TESTPROB_ELBDM_HALO_MERGER = 1012, + TESTPROB_ELBDM_DISK_HEATING = 1013; // program initialization options typedef int OptInit_t; diff --git a/src/Init/Init_TestProb.cpp b/src/Init/Init_TestProb.cpp index b10f304b0f..819bc5cd4a 100644 --- a/src/Init/Init_TestProb.cpp +++ b/src/Init/Init_TestProb.cpp @@ -44,6 +44,7 @@ void Init_TestProb_ELBDM_LSS(); void Init_TestProb_ELBDM_PlaneWave(); void Init_TestProb_ELBDM_Perturbation(); void Init_TestProb_ELBDM_HaloMerger(); +void Init_TestProb_ELBDM_DiskHeating(); @@ -109,6 +110,7 @@ void Init_TestProb() case TESTPROB_ELBDM_PLANE_WAVE : Init_TestProb_ELBDM_PlaneWave(); break; case TESTPROB_ELBDM_PERTURBATION : Init_TestProb_ELBDM_Perturbation(); break; case TESTPROB_ELBDM_HALO_MERGER : Init_TestProb_ELBDM_HaloMerger(); break; + case TESTPROB_ELBDM_DISK_HEATING : Init_TestProb_ELBDM_DiskHeating(); break; default: Aux_Error( ERROR_INFO, "unsupported TESTPROB_ID (%d) !!\n", TESTPROB_ID ); } // switch( TESTPROB_ID ) diff --git a/src/Model_ELBDM/CPU_ELBDMGravity/CPU_ELBDMGravitySolver_HJ.cpp b/src/Model_ELBDM/CPU_ELBDMGravity/CPU_ELBDMGravitySolver_HJ.cpp index 7904b8c49c..fb50b961cf 100644 --- a/src/Model_ELBDM/CPU_ELBDMGravity/CPU_ELBDMGravitySolver_HJ.cpp +++ b/src/Model_ELBDM/CPU_ELBDMGravity/CPU_ELBDMGravitySolver_HJ.cpp @@ -13,7 +13,7 @@ //----------------------------------------------------------------------------------------- // Function : CPU/CUPOT_ELBDMGravitySolver_HamiltonJacobi -// Description : CPU/GPU ELBDM gravity solver for advancing the phase S by S = S - i*Eta*(Phi + Lambda*Rho)*dt +// Description : CPU/GPU ELBDM gravity solver for advancing the phase S by S = S - Eta*(Phi + Lambda*Rho)*dt // // Note : 1. ELBDM gravity solver requires NO potential and fluid ghost zone // --> Optimized performance can be achieved if GRA_GHOST_SIZE==0 (and thus GRA_NXT==PATCH_SIZE) diff --git a/src/Model_ELBDM/ELBDM_GramFE_EvolutionMatrix.cpp b/src/Model_ELBDM/ELBDM_GramFE_EvolutionMatrix.cpp index 6c3efb6437..b4f0477f05 100644 --- a/src/Model_ELBDM/ELBDM_GramFE_EvolutionMatrix.cpp +++ b/src/Model_ELBDM/ELBDM_GramFE_EvolutionMatrix.cpp @@ -95,8 +95,8 @@ void ELBDM_GramFE_ComputeTimeEvolutionMatrix( gramfe_matmul_float (*output)[2 * long double K, Filter, Coeff; gramfe_evo_complex_type ExpCoeff; - const long double filterDecay = (long double) 32.0 * (long double) 2.302585092994046; // decay of k-space filter ( 32 * log(10) ) - const long double filterDegree = (long double) 100; // degree of k-space filter + const long double filterDecay = (long double) 16.0 * (long double) 2.302585092994046; // decay of k-space filter ( 16 * log(10) ) + const long double filterDegree = (long double) 50; // degree of k-space filter const long double kmax = (long double) M_PI / dh; // maximum value of k const long double dk = (long double) + 2.0 * kmax / GRAMFE_FLU_NXT; // k steps in k-space const long double dT = (long double) - 0.5 * dt / Eta; // coefficient in time evolution operator diff --git a/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cpp b/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cpp new file mode 100644 index 0000000000..9f8605fc7f --- /dev/null +++ b/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cpp @@ -0,0 +1,213 @@ +#include "CUPOT.h" +#ifdef __CUDACC__ +#include "CUDA_CheckError.h" +#endif + +#ifdef GRAVITY + +extern double m_22; +extern double CoreRadius; +extern double Cen[3]; + + +// ================================= +// I. Set auxiliary arrays +// ================================= + +#ifndef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : SetExtPotAuxArray_Soliton +// Description : Set the auxiliary arrays ExtPot_AuxArray_Flt/Int[] used by ExtPot_Soliton() +// +// Note : 1. Invoked by Init_ExtPot_Soliton() +// 2. AuxArray_Flt/Int[] have the size of EXT_POT_NAUX_MAX defined in Macro.h (default = 20) +// 3. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : AuxArray_Flt/Int : Floating-point/Integer arrays to be filled up +// Time : Target physical time +// +// Return : AuxArray_Flt/Int[] +//------------------------------------------------------------------------------------------------------- +void SetExtPotAuxArray_Soliton( double AuxArray_Flt[], int AuxArray_Int[], const double Time ) +{ + + const double A = 0.0019/m_22/m_22/POW(CoreRadius,4)*1.0e10*Const_Msun/Const_kpc; // soliton central density in g/cm^3 + const double B = 9.1*0.01; + const double G = Const_NewtonG/UNIT_V/UNIT_V; // convert output potential to code unit + + AuxArray_Flt[0] = Cen[0]; // x coordinate of the external potential center + AuxArray_Flt[1] = Cen[1]; // y ... + AuxArray_Flt[2] = Cen[2]; // z ... + AuxArray_Flt[3] = CoreRadius*Const_kpc; + AuxArray_Flt[4] = A; + AuxArray_Flt[5] = B; + AuxArray_Flt[6] = G; + AuxArray_Flt[7] = UNIT_L; + + if ( MPI_Rank == 0 ) + { + Aux_Message( stdout, "EXT_POT_AUX_ARRAY:\n" ); + Aux_Message( stdout, "=============================================================================\n" ); + Aux_Message( stdout, " CenX = %13.7e\n", Cen[0] ); + Aux_Message( stdout, " CenY = %13.7e\n", Cen[1] ); + Aux_Message( stdout, " Cenz = %13.7e\n", Cen[2] ); + Aux_Message( stdout, " CoreRadius = %13.7e\n", CoreRadius ); + Aux_Message( stdout, " UNIT_L = %13.7e\n", UNIT_L ); + Aux_Message( stdout, " UNIT_V = %13.7e\n", UNIT_V ); + Aux_Message( stdout, " m_22 = %13.7e\n", m_22 ); + Aux_Message( stdout, "=============================================================================\n" ); + } + +} // FUNCTION : SetExtPotAuxArray_Soliton +#endif // #ifndef __CUDACC__ + + + +// ================================= +// II. Specify external potential +// ================================= + +//----------------------------------------------------------------------------------------- +// Function : ExtPot_Soliton +// Description : Calculate the external potential at the given coordinates and time +// +// Note : 1. This function is shared by CPU and GPU +// 2. Auxiliary arrays UserArray_Flt/Int[] are set by SetExtPotAuxArray_SOliton(), where +// UserArray_Flt[0] = x coordinate of the external potential center +// UserArray_Flt[1] = y ... +// UserArray_Flt[2] = z .. +// UserArray_Flt[3] = gravitational_constant*point_source_mass +// 3. Currently it does not support the soften length +// 4. GenePtr has the size of EXT_POT_NGENE_MAX defined in Macro.h (default = 6) +// +// Parameter : x/y/z : Target spatial coordinates +// Time : Target physical time +// UserArray_Flt/Int : User-provided floating-point/integer auxiliary arrays +// Usage : Different usages of external potential when computing total potential on level Lv +// --> EXT_POT_USAGE_ADD : add external potential on Lv +// EXT_POT_USAGE_SUB : subtract external potential for preparing self-gravity potential on Lv-1 +// EXT_POT_USAGE_SUB_TINT: like SUB but for temporal interpolation +// --> This parameter is useless in most cases +// PotTable : 3D potential table used by EXT_POT_TABLE +// GenePtr : Array of pointers for general potential tables +// +// Return : External potential at (x,y,z,Time) +//----------------------------------------------------------------------------------------- +GPU_DEVICE_NOINLINE +static real ExtPot_Soliton( const double x, const double y, const double z, const double Time, + const double UserArray_Flt[], const int UserArray_Int[], + const ExtPotUsage_t Usage, const real PotTable[], void **GenePtr ) +{ + + const double Center[3] = { UserArray_Flt[0], UserArray_Flt[1], UserArray_Flt[2] }; + const double r_sol = UserArray_Flt[3]; + const double dx = (x - Center[0]); + const double dy = (y - Center[1]); + const double dz = (z - Center[2]); + const double A = UserArray_Flt[4]; + const double B = UserArray_Flt[5]; + const double G = UserArray_Flt[6]; + const double unit_l = UserArray_Flt[7]; + const double r = SQRT( dx*dx + dy*dy + dz*dz )*unit_l; + const double Y = r/r_sol; + +// soliton potential, consistent with Eq. [6] in Chiang et al., PRD 103, 103019 (2021) +// --> this form may suffer from large numerical errors at r->0 as both the numerator and denominator approach 0 + double soliton_potential = -M_PI*pow(Y,2.)*A*G/53760./pow(B,1.5) + *( pow(B,0.5)/pow((1.+B*pow(Y,2.)),6.) + *( 11895. + 36685.*B*pow(Y,2.) + + 55638.*pow(B,2.)*pow(Y,4.) + 45738.*pow(B,3.)*pow(Y,6.) + + 19635.*pow(B,4.)*pow(Y,8.) + 3465.*pow(B,5.)*pow(Y,10.) ) + + 3465.*atan(pow(B,0.5)*Y)/Y ); + + return (real)soliton_potential; + +} // FUNCTION : ExtPot_Soliton + + + +// ================================= +// III. Set initialization functions +// ================================= + +#ifdef __CUDACC__ +# define FUNC_SPACE __device__ static +#else +# define FUNC_SPACE static +#endif + +FUNC_SPACE ExtPot_t ExtPot_Ptr = ExtPot_Soliton; + +//----------------------------------------------------------------------------------------- +// Function : SetCPU/GPUExtPot_Soliton +// Description : Return the function pointers of the CPU/GPU external potential routines +// +// Note : 1. Invoked by Init_ExtPot_Soliton() +// 2. Must obtain the CPU and GPU function pointers by **separate** routines +// since CPU and GPU functions are compiled completely separately in GAMER +// --> In other words, a unified routine like the following won't work +// +// SetExtPot_Soliton( ExtPot_t &CPUExtPot_Ptr, ExtPot_t &GPUExtPot_Ptr ) +// +// Parameter : CPU/GPUExtPot_Ptr (call-by-reference) +// +// Return : CPU/GPUExtPot_Ptr +//----------------------------------------------------------------------------------------- +#ifdef __CUDACC__ +__host__ +void SetGPUExtPot_Soliton( ExtPot_t &GPUExtPot_Ptr ) +{ + CUDA_CHECK_ERROR( cudaMemcpyFromSymbol( &GPUExtPot_Ptr, ExtPot_Ptr, sizeof(ExtPot_t) ) ); +} + +#else // #ifdef __CUDACC__ + +void SetCPUExtPot_Soliton( ExtPot_t &CPUExtPot_Ptr ) +{ + CPUExtPot_Ptr = ExtPot_Ptr; +} + +#endif // #ifdef __CUDACC__ ... else ... + + + +#ifndef __CUDACC__ + +// local function prototypes +void SetExtPotAuxArray_Soliton( double [], int [], const double ); +void SetCPUExtPot_Soliton( ExtPot_t & ); +#ifdef GPU +void SetGPUExtPot_Soliton( ExtPot_t & ); +#endif + +//----------------------------------------------------------------------------------------- +// Function : Init_ExtPot_Soliton +// Description : Initialize external potential +// +// Note : 1. Set auxiliary arrays by invoking SetExtPotAuxArray_*() +// --> They will be copied to GPU automatically in CUAPI_SetConstMemory() +// 2. Set the CPU/GPU external potential major routines by invoking SetCPU/GPUExtPot_*() +// 3. Invoked by Init_ExtAccPot() +// --> Enable it by linking to the function pointer "Init_ExtPot_Ptr" +// 4. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : None +// +// Return : None +//----------------------------------------------------------------------------------------- +void Init_ExtPot_Soliton() +{ + + SetExtPotAuxArray_Soliton( ExtPot_AuxArray_Flt, ExtPot_AuxArray_Int, Time[0] ); + SetCPUExtPot_Soliton( CPUExtPot_Ptr ); +# ifdef GPU + SetGPUExtPot_Soliton( GPUExtPot_Ptr ); +# endif + +} // FUNCTION : Init_ExtPot_Soliton + +#endif // #ifndef __CUDACC__ + + + +#endif // #ifdef GRAVITY diff --git a/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cu b/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cu new file mode 120000 index 0000000000..df762c8fb8 --- /dev/null +++ b/src/TestProblem/ELBDM/DiskHeating/ExtPot_Soliton.cu @@ -0,0 +1 @@ +ExtPot_Soliton.cpp \ No newline at end of file diff --git a/src/TestProblem/ELBDM/DiskHeating/Init_TestProb_ELBDM_DiskHeating.cpp b/src/TestProblem/ELBDM/DiskHeating/Init_TestProb_ELBDM_DiskHeating.cpp new file mode 100644 index 0000000000..502e6db237 --- /dev/null +++ b/src/TestProblem/ELBDM/DiskHeating/Init_TestProb_ELBDM_DiskHeating.cpp @@ -0,0 +1,733 @@ +#include "GAMER.h" +#include "TestProb.h" + +#ifdef SUPPORT_GSL +#include +#include +#endif + + +static void SetGridIC( real fluid[], const double x, const double y, const double z, const double Time, + const int lv, double AuxArray[] ); +static void BC( real Array[], const int ArraySize[], real fluid[], const int NVar_Flu, + const int GhostSize, const int idx[], const double pos[], const double Time, + const int lv, const int TFluVarIdxList[], double AuxArray[] ); +static void End_DiskHeating(); +static double Get_Dispersion( double r ); +static double Halo_Density( double r ); + +// problem-specific global variables +// ======================================================================================= +static RandomNumber_t *RNG = NULL; + double Cen[3]; // center +static bool AddFixedHalo; // add a fixed halo, must enable OPT__FREEZE_FLUID +static bool HaloUseTable; // 0 = from analytical profile, 1 = from table + double m_22; // ELBDM particle mass, used for soliton profile of the fixed halo + double CoreRadius; // soliton radius of the fixed halo (in kpc) +static double Rho_0; // halo rho_0 (in 1.0e+10 Msun*kpc^-3) +static double Rs; // halo Rs (in kpc) +static double Alpha; // dimensionless, used for alpha-beta-gamma density profile of the fixed halo +static double Beta; // dimensionless, used for alpha-beta-gamma density profile of the fixed halo +static double Gamma; // dimensionless, used for alpha-beta-gamma density profile of the fixed halo +static char DensTableFile[MAX_STRING]; // fixed halo density profile filename +static double *DensTable = NULL; // density table, radius should be in kpc and density should be in g/cm^3 +static int DensTable_Nbin; // number of bins of density table +static bool AddParWhenRestart; // add a new disk to an existing snapshot, must enable OPT__RESTART_RESET +static bool AddParWhenRestartByFile; // add a new disk via PAR_IC +static long AddParWhenRestartNPar; // particle number of the new disk +static int NewDisk_RSeed; // random seed for setting new disk particle position and velocity +static double Disk_Mass; // thin disk total mass (code unit) +static double Disk_R; // thin disk scale radius (code unit) +static char DispTableFile[MAX_STRING]; // velocity dispersion table filename used for thin disk +static double *DispTable = NULL; // velocity dispersion table, radius should be in kpc and dispersion should be in km/s +static int DispTable_Nbin; // number of bins of velocity dispersion table +// ======================================================================================= + +#ifdef PARTICLE +void Par_Init_ByFunction_DiskHeating( const long NPar_ThisRank, const long NPar_AllRank, + real *ParMass, real *ParPosX, real *ParPosY, real *ParPosZ, + real *ParVelX, real *ParVelY, real *ParVelZ, real *ParTime, + real *ParType, real *AllAttribute[PAR_NATT_TOTAL]); +static void Init_NewDiskRestart(); +static void Init_NewDiskVelocity(); +#endif +#ifdef GRAVITY +void Init_ExtPot_Soliton(); +#endif + + + +//------------------------------------------------------------------------------------------------------- +// Function : Validate +// Description : Validate the compilation flags and runtime parameters for this test problem +// +// Note : None +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Validate() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Validating test problem %d ...\n", TESTPROB_ID ); + + +// errors +# if ( MODEL != ELBDM ) + Aux_Error( ERROR_INFO, "MODEL != ELBDM !!\n" ); +# endif + +# ifndef GRAVITY + Aux_Error( ERROR_INFO, "GRAVITY must be enabled !!\n" ); +# endif + +# ifndef PARTICLE + Aux_Error( ERROR_INFO, "PARTICLE must be enabled !!\n" ); +# endif + +# ifdef COMOVING + Aux_Error( ERROR_INFO, "COMOVING must be disabled !!\n" ); +# endif + + +// warnings + if ( MPI_Rank == 0 ) + { + if ( FLAG_BUFFER_SIZE < 5 ) + Aux_Message( stderr, "WARNING : it's recommended to set FLAG_BUFFER_SIZE >= 5 for this test !!\n" ); + + if ( !OPT__RECORD_CENTER ) + Aux_Message( stderr, "WARNING : it's recommended to set OPT__RECORD_CENTER = 1 for this test !!\n" ); + } // if ( MPI_Rank == 0 ) + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Validating test problem %d ... done\n", TESTPROB_ID ); + +} // FUNCTION : Validate + + + +#if ( MODEL == ELBDM ) +//------------------------------------------------------------------------------------------------------- +// Function : SetParameter +// Description : Load and set the problem-specific runtime parameters +// +// Note : 1. Filename is set to "Input__TestProb" by default +// 2. Major tasks in this function: +// (1) load the problem-specific runtime parameters +// (2) set the problem-specific derived parameters +// (3) reset other general-purpose parameters if necessary +// (4) make a note of the problem-specific parameters +// 3. Must NOT call any EoS routine here since it hasn't been initialized at this point +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void SetParameter() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Setting runtime parameters ...\n" ); + + +// (1) load the problem-specific runtime parameters + const char FileName[] = "Input__TestProb"; + ReadPara_t *ReadPara = new ReadPara_t; + +// (1-1) add parameters in the following format: +// --> note that VARIABLE, DEFAULT, MIN, and MAX must have the same data type +// --> some handy constants (e.g., Useless_bool, Eps_double, NoMin_int, ...) are defined in "include/ReadPara.h" +// ******************************************************************************************************************************** +// ReadPara->Add( "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); +// ******************************************************************************************************************************** + ReadPara->Add( "CenX", &Cen[0], NoDef_double, NoMin_double, NoMax_double ); + ReadPara->Add( "CenY", &Cen[1], NoDef_double, NoMin_double, NoMax_double ); + ReadPara->Add( "CenZ", &Cen[2], NoDef_double, NoMin_double, NoMax_double ); + ReadPara->Add( "AddFixedHalo", &AddFixedHalo, false, Useless_bool, Useless_bool ); + ReadPara->Add( "HaloUseTable", &HaloUseTable, false, Useless_bool, Useless_bool ); + ReadPara->Add( "m_22", &m_22, 0.4, Eps_double, NoMax_double ); + ReadPara->Add( "CoreRadius", &CoreRadius, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Rho_0", &Rho_0, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Rs", &Rs, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Alpha", &Alpha, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Beta", &Beta, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Gamma", &Gamma, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "DensTableFile", DensTableFile, NoDef_str, Useless_str, Useless_str ); + ReadPara->Add( "AddParWhenRestart", &AddParWhenRestart, false, Useless_bool, Useless_bool ); + ReadPara->Add( "AddParWhenRestartByFile", &AddParWhenRestartByFile, true, Useless_bool, Useless_bool ); + ReadPara->Add( "AddParWhenRestartNPar", &AddParWhenRestartNPar, (long)0, (long)0, NoMax_long ); + ReadPara->Add( "NewDisk_RSeed", &NewDisk_RSeed, 1002, 0, NoMax_int ); + ReadPara->Add( "Disk_Mass", &Disk_Mass, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Disk_R", &Disk_R, 1.0, Eps_double, NoMax_double ); + ReadPara->Add( "DispTableFile", DispTableFile, NoDef_str, Useless_str, Useless_str ); + + + ReadPara->Read( FileName ); + + delete ReadPara; + +// (1-2) set the default values + +// (1-3) check the runtime parameters + + +// (2) set the problem-specific derived parameters +// use density table as background fixed halo profile + if ( HaloUseTable == 1 ) { + +// read the density table + const bool RowMajor_No = false; // load data into the column-major order + const bool AllocMem_Yes = true; // allocate memory for DensTable + const int NCol = 2; // total number of columns to load + const int Col[NCol] = {0, 1}; // target columns: (radius, density) + + DensTable_Nbin = Aux_LoadTable( DensTable, DensTableFile, NCol, Col, RowMajor_No, AllocMem_Yes ); + + double *DensTable_r = DensTable + 0*DensTable_Nbin; + double *DensTable_d = DensTable + 1*DensTable_Nbin; + +// convert to code unit and log-log scale (input units of radius and density should be fixed to kpc and g/cm^3 respectively) + for (int b=0; b a helper macro PRINT_WARNING is defined in TestProb.h + const long End_Step_Default = 2500; + const double End_T_Default = 2.5e-1; + + if ( END_STEP < 0 ) { + END_STEP = End_Step_Default; + PRINT_RESET_PARA( END_STEP, FORMAT_LONG, "" ); + } + + if ( END_T < 0.0 ) { + END_T = End_T_Default; + PRINT_RESET_PARA( END_T, FORMAT_REAL, "" ); + } + + +// (4) make a note + if ( MPI_Rank == 0 ) + { + Aux_Message( stdout, "=============================================================================\n" ); + Aux_Message( stdout, " test problem ID = %d\n", TESTPROB_ID ); + Aux_Message( stdout, " CenX = %13.7e\n", Cen[0] ); + Aux_Message( stdout, " CenY = %13.7e\n", Cen[1] ); + Aux_Message( stdout, " Cenz = %13.7e\n", Cen[2] ); + Aux_Message( stdout, " AddFixedHalo = %d\n", AddFixedHalo ); + Aux_Message( stdout, " HaloUseTable = %d\n", HaloUseTable ); + Aux_Message( stdout, " m_22 = %13.7e\n", m_22 ); + Aux_Message( stdout, " CoreRadius = %13.7e\n", CoreRadius ); + Aux_Message( stdout, " Rho_0 = %13.7e\n", Rho_0 ); + Aux_Message( stdout, " Rs = %13.7e\n", Rs ); + Aux_Message( stdout, " Alpha = %13.7e\n", Alpha ); + Aux_Message( stdout, " Beta = %13.7e\n", Beta ); + Aux_Message( stdout, " Gamma = %13.7e\n", Gamma ); + Aux_Message( stdout, " DensTableFile = %s\n", DensTableFile ); + Aux_Message( stdout, " AddParWhenRestart = %d\n", AddParWhenRestart ); + Aux_Message( stdout, " AddParWhenRestartByFile = %d\n", AddParWhenRestartByFile ); + Aux_Message( stdout, " AddParWhenRestartNPar = %d\n", AddParWhenRestartNPar ); + Aux_Message( stdout, " NewDisk_RSeed = %d\n", NewDisk_RSeed ); + Aux_Message( stdout, " Disk_Mass = %13.7e\n", Disk_Mass ); + Aux_Message( stdout, " Disk_R = %13.7e\n", Disk_R ); + Aux_Message( stdout, " DispTableFile = %s\n", DispTableFile ); + Aux_Message( stdout, "=============================================================================\n" ); + } + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Setting runtime parameters ... done\n" ); + +} // FUNCTION : SetParameter + + + +//------------------------------------------------------------------------------------------------------- +// Function : Halo_Density +// Description : alpha-beta-gamma density profile for fixed background halo with soliton as function of radius +// +// Note : 1. r should have unit in kpc +// 2. Returned density is in 1.0e10*Msun/kpc^3 +//------------------------------------------------------------------------------------------------------- +double Halo_Density( double r ) +{ + + double rho_halo = Rho_0 / pow( r/Rs, Alpha ) / pow( 1.0 + pow(r/Rs,Beta), (Gamma-Alpha)/Beta ); + if ( fabs(rho_halo) < __FLT_MIN__ ) rho_halo = 0.0; + + double rho_soliton = 0.0019 / m_22 / m_22 * pow( 1.0/CoreRadius/pow(1.0 + 0.091*pow(r/CoreRadius,2.0), 2.0), 4.0 ); + if ( fabs(rho_soliton) < __FLT_MIN__ ) rho_soliton = 0.0; + + double rho_max = 0.0019 / m_22 / m_22 * pow( 1.0 / CoreRadius, 4.0 ); + + if ( (rho_halo + rho_soliton) > rho_max ) return rho_max; + else return (rho_halo + rho_soliton); + +} // FUNCTION : Halo_Density + + + +//------------------------------------------------------------------------------------------------------- +// Function : SetGridIC +// Description : Set the problem-specific initial condition on grids +// +// Note : 1. This function may also be used to estimate the numerical errors when OPT__OUTPUT_USER is enabled +// --> In this case, it should provide the analytical solution at the given "Time" +// 2. This function will be invoked by multiple OpenMP threads when OPENMP is enabled +// (unless OPT__INIT_GRID_WITH_OMP is disabled) +// --> Please ensure that everything here is thread-safe +// 3. Even when DUAL_ENERGY is adopted for HYDRO, one does NOT need to set the dual-energy variable here +// --> It will be calculated automatically +// 4. For MHD, do NOT add magnetic energy (i.e., 0.5*B^2) to fluid[ENGY] here +// --> It will be added automatically later +// +// Parameter : fluid : Fluid field to be initialized +// x/y/z : Physical coordinates +// Time : Physical time +// lv : Target refinement level +// AuxArray : Auxiliary array +// +// Return : fluid +//------------------------------------------------------------------------------------------------------- +void SetGridIC( real fluid[], const double x, const double y, const double z, const double Time, + const int lv, double AuxArray[] ) +{ + +// set the output array + double dens = 0.0; + + if ( AddFixedHalo && OPT__FREEZE_FLUID ) // add a fixed halo at the background + { + const double dx = x - Cen[0]; + const double dy = y - Cen[1]; + const double dz = z - Cen[2]; + const double r = SQRT( dx*dx + dy*dy + dz*dz ); + + if ( HaloUseTable ) + { + const double *DensTable_r = DensTable + 0*DensTable_Nbin; + const double *DensTable_d = DensTable + 1*DensTable_Nbin; + + if ( r < exp( DensTable_r[0] ) ) + dens = exp( DensTable_d[0] ); + + else if ( r > exp( DensTable_r[DensTable_Nbin-1] ) ) + dens = 0.0; + + else + dens = exp( Mis_InterpolateFromTable(DensTable_Nbin, DensTable_r, DensTable_d, log(r)) ); + } + else + { + const double Unit_D_GALIC = 1.0e10*Const_Msun/CUBE(Const_kpc); + const double r_in_kpc = r*UNIT_L/Const_kpc; + dens = Halo_Density( r_in_kpc )*Unit_D_GALIC/UNIT_D; + } + } + + fluid[DENS] = (real)dens; + fluid[REAL] = sqrt( fluid[DENS] ); + fluid[IMAG] = 0.0; + +} // FUNCTION : SetGridIC + + + +//------------------------------------------------------------------------------------------------------- +// Function : BC +// Description : Boundary conditions +// +// Note : +//------------------------------------------------------------------------------------------------------- +void BC( real Array[], const int ArraySize[], real fluid[], const int NVar_Flu, + const int GhostSize, const int idx[], const double pos[], const double Time, + const int lv, const int TFluVarIdxList[], double AuxArray[] ) +{ + +// simply call the IC function + SetGridIC( fluid, pos[0], pos[1], pos[2], Time, lv, AuxArray ); + +} // FUNCTION : BC + + + +//------------------------------------------------------------------------------------------------------- +// Function : End_DiskHeating +// Description : Free memory before terminating the program +// +// Note : 1. Linked to the function pointer "End_User_Ptr" to replace "End_User()" +//------------------------------------------------------------------------------------------------------- +void End_DiskHeating() +{ + + delete [] DensTable; + delete [] DispTable; + +} // FUNCTION : End_DiskHeating + + + +#ifdef PARTICLE +#ifdef GRAVITY +//------------------------------------------------------------------------------------------------------- +// Function : Init_NewDiskRestart() +// Description : Add a new disk from an existing snapshot +// +// Note : Must enable OPT__RESTART_RESET and AddParWhenRestart +//------------------------------------------------------------------------------------------------------- +void Init_NewDiskRestart() +{ + + if ( AddParWhenRestart && !OPT__RESTART_RESET ) + Aux_Error( ERROR_INFO, "OPT__RESTART_RESET should be turned on to enable AddParWhenRestart !!\n" ); + + if ( amr->Par->Init != PAR_INIT_BY_RESTART || !OPT__RESTART_RESET || !AddParWhenRestart ) return; + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + + + const long NNewPar = ( MPI_Rank == 0 ) ? AddParWhenRestartNPar : 0; + const long NPar_AllRank = NNewPar; + real *NewParAtt[PAR_NATT_TOTAL]; + + for (int v=0; v = %ld != expect = %ld !!\n", + FileName, FileSize, ExpectSize ); + fclose( FileTemp ); + + Aux_Message( stdout, " Loading data ... " ); + + real *ParData_AllRank = new real [ NPar_AllRank*NParAtt ]; + +// note that fread() may fail for large files if sizeof(size_t) == 4 instead of 8 + FILE *File = fopen( FileName, "rb" ); + + for (int v=0; vSetSeed( 0, NewDisk_RSeed ); + + for (long p=0; pGetValue( 0, 0.0, 1.0 ); + R = 1.0; + + do + { + f = (1.0 + R) * exp(-R) + Ran - 1.0; + f_ = -R * exp(-R); + Rold = R; + R = R - f / f_; + } + while( fabs(R - Rold) / R > 1e-7 ); + + RanR = Disk_R*R; + phi = 2*M_PI*RNG->GetValue( 0, 0.0, 1.0 ); + RanVec[0] = RanR*cos(phi); + RanVec[1] = RanR*sin(phi); + RanVec[2] = 0; + for (int d=0; d<3; d++) Pos_AllRank[d][p] = RanVec[d] + Cen[d]; + + } // for ( long p = 0; p < NPar_AllRank; p++ ) + + delete RNG; + } // if ( MPI_Rank == 0 ) + } // if ( !AddParWhenRestartByFile ) + + +// add particles here + Par_AddParticleAfterInit( NNewPar, NewParAtt ); + +// free memory + for (int v=0; vLB->Par_Weight; +# else + const double Par_Weight = 0.0; +# endif +# ifdef LOAD_BALANCE + const UseLBFunc_t UseLB = USELB_YES; +# else + const UseLBFunc_t UseLB = USELB_NO; +# endif + + for (int lv=0; lvPar->Init != PAR_INIT_BY_RESTART || !OPT__RESTART_RESET || + !AddParWhenRestart || AddParWhenRestartByFile ) return; + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + +# ifndef SUPPORT_GSL + Aux_Error( ERROR_INFO, "SUPPORT_GSL must be enabled when AddParWhenRestart=1 and AddParWhenRestartByFile=0 !!\n" ); +# endif + + const long NPar_ThisRank = amr->Par->NPar_AcPlusInac; + real *ParType = amr->Par->Type; + real *ParPos[3] = { amr->Par->PosX, amr->Par->PosY, amr->Par->PosZ }; + real *ParVel[3] = { amr->Par->VelX, amr->Par->VelY, amr->Par->VelZ }; + real *ParAcc[3] = { amr->Par->AccX, amr->Par->AccY, amr->Par->AccZ }; + + real ParRadius[2]; + real NormParRadius[2]; + double V_acc, RanV[3], sigma; + double ParR; + +# ifdef SUPPORT_GSL +// initialize the RNG + gsl_rng *random_generator; + random_generator = gsl_rng_alloc(gsl_rng_ranlxd1); + gsl_rng_set(random_generator, NewDisk_RSeed + MPI_Rank); + + for (long p=0; p DispTable_r[DispTable_Nbin-1] ) + disp = DispTable_d[DispTable_Nbin-1]; + + else + disp = Mis_InterpolateFromTable( DispTable_Nbin, DispTable_r, DispTable_d, r ); + + return disp; + +} // FUNCTION : Get_Dispersion +#endif // #ifdef GRAVITY +#endif // #ifdef PARTICLE +#endif // #if ( MODEL == ELBDM ) + + + +//------------------------------------------------------------------------------------------------------- +// Function : Init_TestProb_ELBDM_DiskHeating +// Description : Test problem initializer +// +// Note : None +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Init_TestProb_ELBDM_DiskHeating() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + + +// validate the compilation flags and runtime parameters + Validate(); + + +# if ( MODEL == ELBDM ) +// set the problem-specific runtime parameters + SetParameter(); + + +// set the function pointers + Init_Function_User_Ptr = SetGridIC; + BC_User_Ptr = BC; + End_User_Ptr = End_DiskHeating; +# ifdef GRAVITY + Init_ExtPot_Ptr = Init_ExtPot_Soliton; +# endif +# ifdef PARTICLE + Par_Init_ByFunction_Ptr = Par_Init_ByFunction_DiskHeating; + Init_User_Ptr = Init_NewDiskRestart; + Init_User_AfterPoisson_Ptr = Init_NewDiskVelocity; +# endif +# endif // #if ( MODEL == ELBDM ) + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ... done\n", __FUNCTION__ ); + +} // FUNCTION : Init_TestProb_ELBDM_DiskHeating diff --git a/src/TestProblem/ELBDM/DiskHeating/Par_Init_ByFunction_DiskHeating.cpp b/src/TestProblem/ELBDM/DiskHeating/Par_Init_ByFunction_DiskHeating.cpp new file mode 100644 index 0000000000..40b53e18ec --- /dev/null +++ b/src/TestProblem/ELBDM/DiskHeating/Par_Init_ByFunction_DiskHeating.cpp @@ -0,0 +1,142 @@ +#include "GAMER.h" +#ifdef PARTICLE + + + + +//------------------------------------------------------------------------------------------------------- +// Function : Par_Init_ByFile +// Description : Initialize particle attributes from a file +// +// Note : 1. Invoked by Init_GAMER() +// 2. Periodicity should be taken care of in this function +// --> No particles should lie outside the simulation box even for the periodic BC +// 3. Particles lying outside the active region will be removed later by Par_Aux_InitCheck() +// if non-periodic B.C. is adopted +// 4. Particles loaded here are only temporarily stored in this rank +// --> They will be redistributed when calling Par_FindHomePatch_UniformGrid() +// and LB_Init_LoadBalance() +// --> So there is no constraint on which particles should be set by this function +// 5. Currently the target file name is fixed to "PAR_IC" +// 6. The data format of the PAR_IC file is controlled by the runtime parameter "PAR_IC_FORMAT" +// --> PAR_IC_FORMAT_ATT_ID: [particle attribute][particle id] in a row-major order +// PAR_IC_FORMAT_ID_ATT: [particle id][particle attribute] in a row-major order +// 7 Currently it only loads particle mass, position x/y/z, and velocity x/y/z +// (and must be in the same order of PAR_MASS, PAR_POSX/Y/Z, and PAR_VELX/Y/Z) +// --> The mass of all particles can be set to PAR_IC_MASS instead (by having PAR_IC_MASS>=0.0) +// --> In this case, the PAR_IC file should exclude the partice mass data +// 8. For LOAD_BALANCE, the number of particles in each rank must be set in advance +// --> Currently it's set by Init_Parallelization() +// +// Parameter : None +// +// Return : ParMass, ParPosX/Y/Z, ParVelX/Y/Z, ParTime, AllAttribute +//------------------------------------------------------------------------------------------------------- +void Par_Init_ByFunction_DiskHeating( const long NPar_ThisRank, const long NPar_AllRank, + real *ParMass, real *ParPosX, real *ParPosY, real *ParPosZ, + real *ParVelX, real *ParVelY, real *ParVelZ, real *ParTime, + real *ParType, real *AllAttribute[PAR_NATT_TOTAL]) +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + + + const char FileName[] = "PAR_IC"; + const long NParAllRank = amr->Par->NPar_Active_AllRank; + long NParThisRank = amr->Par->NPar_AcPlusInac; // cannot be "const" due to MPI_Allgather() + const int NParAtt = 8; // mass, pos*3, vel*3, type + + +// check + if ( !Aux_CheckFileExist(FileName) ) + Aux_Error( ERROR_INFO, "file \"%s\" does not exist !!\n", FileName ); + + FILE *FileTemp = fopen( FileName, "rb" ); + + fseek( FileTemp, 0, SEEK_END ); + + const long ExpectSize = long(NParAtt)*NParAllRank*sizeof(real); + const long FileSize = ftell( FileTemp ); + if ( FileSize != ExpectSize ) + Aux_Error( ERROR_INFO, "size of the file <%s> = %ld != expect = %ld !!\n", + FileName, FileSize, ExpectSize ); + + fclose( FileTemp ); + + MPI_Barrier( MPI_COMM_WORLD ); + + +// set the file offset for this rank + long NPar_EachRank[MPI_NRank], NPar_Check=0, FileOffset=0; + + MPI_Allgather( &NParThisRank, 1, MPI_LONG, NPar_EachRank, 1, MPI_LONG, MPI_COMM_WORLD ); + +// check if the total number of particles is correct + for (int r=0; rAttribute[] are the same + ParMass[p] = ParData1[0]; + ParPosX[p] = ParData1[1]; + ParPosY[p] = ParData1[2]; + ParPosZ[p] = ParData1[3]; + ParVelX[p] = ParData1[4]; + ParVelY[p] = ParData1[5]; + ParVelZ[p] = ParData1[6]; + ParType[p] = ParData1[7]; // 1=CDM halo, 2=disk + +// synchronize all particles to the physical time at the base level + amr->Par->Time[p] = Time[0]; + } + + delete [] ParData_ThisRank; + delete [] ParData1; + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "done\n" ); + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ... done\n", __FUNCTION__ ); + +} // FUNCTION : Par_Init_ByFunction_DiskHeating + + + +#endif // #ifdef PARTICLE