diff --git a/extern/Doppler-spectrum/exact_flux_density.csv b/extern/Doppler-spectrum/exact_flux_density.csv new file mode 100644 index 000000000..b52217a02 --- /dev/null +++ b/extern/Doppler-spectrum/exact_flux_density.csv @@ -0,0 +1,1025 @@ +0.001,0.0039423068928400645 +0.0010113065323601315,0.0040319584961523835 +0.0010227409023942737,0.004123648857243841 +0.0010343045555032248,0.004217424339207206 +0.0010459989534302536,0.004313332359659339 +0.0010578255744458764,0.004411421414328767 +0.0010697859135347235,0.004511741101911226 +0.0010818814825845169,0.004614342148997213 +0.0010941138105771857,0.004719276435662793 +0.0011064844437821437,0.004826597021887879 +0.0011189949459517487,0.004936358174238318 +0.0011316468985189756,0.0050486153933554 +0.0011444419007973232,0.0051634254419988725 +0.0011573815701829784,0.005280846373764771 +0.0011704675423592722,0.0054009375624989465 +0.001183701471503441,0.005523759732106628 +0.00119708503049573,0.0056493749876159 +0.001210619911130859,0.005777846846178777 +0.0012243078243318797,0.005909240269538365 +0.0012381505003664502,0.006043621696610636 +0.0012521496890655567,0.006181059077304426 +0.001266307160044705,0.006321621906625406 +0.001280624702927617,0.0064653812600767935 +0.0012951041275724518,0.006612409829422612 +0.0013097472643005899,0.006762781959486367 +0.0013245559641279984,0.006916573685818445 +0.0013395320989992166,0.0070738627730657895 +0.0013546775620239863,0.007234728754207088 +0.0013699942677165547,0.0073992529710184786 +0.0013854841522376866,0.007567518615007686 +0.0014011491736394114,0.0077396107694860635 +0.0014169913121125372,0.007915616452631018 +0.001433012570236963,0.008095624661531695 +0.0014492149732348223,0.00827972641716773 +0.001465600569226489,0.008468014810340256 +0.0014821714294894757,0.008660585048859215 +0.0014989296487202609,0.00885753450559106 +0.001515877345299077,0.009058962767867818 +0.0015330166615576916,0.009264971687563095 +0.0015503497640502143,0.009475665432842674 +0.0015678788438269703,0.009691150540728092 +0.0015856061167114658,0.009911535971005773 +0.0016035338235804864,0.010136933161221861 +0.0016216642306473648,0.010367456083174114 +0.001639999629748447,0.010603221300451826 +0.0016585423386328015,0.010844348027429005 +0.0016772947012552017,0.01109095818941948 +0.0016962590880724208,0.011343176484549182 +0.001715437896342879,0.011601130446570112 +0.0017348335504296759,0.011864950509559066 +0.0017544485021070508,0.01213477007367672 +0.0017742852308703086,0.012410725572782836 +0.0017943462442492472,0.012692956543296179 +0.0018146340781251319,0.012981605694857668 +0.0018351512970512514,0.013276818982435011 +0.0018559004945770986,0.013578745680110929 +0.0018768842935762188,0.013887538456569901 +0.0018981053465777608,0.014203353452399637 +0.0019195663361017813,0.014526350358750597 +0.0019412699749983351,0.014856692498556664 +0.001963219006790406,0.015194546908567584 +0.001985416206020707,0.015540084424304637 +0.0020078643786024097,0.015893479766152675 +0.0020305663621738327,0.016254911627665434 +0.002053525026457146,0.016624562766224926 +0.002076743273621124,0.017002620095098397 +0.0021002240386480066,0.017389274778189818 +0.0021239702897045067,0.017784722326609416 +0.0021479850285170085,0.018189162697468023 +0.0021722712907510143,0.018602800395322923 +0.0021968321463948755,0.019025844574944023 +0.0022216707001478666,0.019458509147655765 +0.0022467900918126445,0.019901012889345496 +0.002272193496692147,0.020353579550925645 +0.0022978841259909772,0.020816437971623837 +0.002323865227221327,0.02128982219455946 +0.0023501400846134895,0.02177397158528409 +0.002376712019531014,0.02226913095275202 +0.0024035843908905554,0.022775550672805567 +0.0024307605955864666,0.023293486815292214 +0.0024582440689201973,0.02382320127310812 +0.002486038285034545,0.024364961894969218 +0.0025141467573528142,0.024919042620544753 +0.0025425730390229435,0.025485723618950778 +0.0025713207233666545,0.026065291430891223 +0.002600393444333677,0.026658039112926857 +0.002629794876961109,0.027264266386160925 +0.0026595287378379786,0.027884279787380693 +0.002689598785575043,0.028518392824287527 +0.002720008821279918,0.029166926134107717 +0.002750762689037563,0.029830207645431053 +0.0027818642763962086,0.030508572744232843 +0.0028133175148587763,0.03120236444341649 +0.002845126380379852,0.03191193355624197 +0.002877294893868281,0.032637638873746246 +0.0029098271216954435,0.03337984734591945 +0.0029427271762092816,0.03413893426767201 +0.0029759992162541305,0.034915283468387945 +0.003009647447696434,0.03570928750574236 +0.0030436761239564015,0.03652134786476948 +0.003078089546545674,0.03735187516016217 +0.0031128920656110755,0.03820128934463451 +0.0031480880804845043,0.03907001992060929 +0.0031836820402390466,0.039958506157947785 +0.003219678444251379,0.0408671973155199 +0.0032560818427705253,0.04179655286885343 +0.003292896837493047,0.04274704274202892 +0.003330128082144737,0.04371914754572802 +0.00336778028306889,0.044713358819712314 +0.0034058581998212207,0.045730179281954386 +0.003444366645771519,0.0467701230825304 +0.0034833104887120924,0.04783371606337225 +0.003522694651473102,0.04892149602449802 +0.003562524112544844,0.05003401299584514 +0.0036028039067070817,0.051171829515294115 +0.0036435391256654736,0.05233552091315337 +0.003684734918695216,0.05352567560316977 +0.00372639649329195,0.05474289537984292 +0.0037685291158300366,0.05598779572276814 +0.003811138112228267,0.0572610061079891 +0.0038542288686231065,0.05856317032618657 +0.003897806832049547,0.059894946807999914 +0.003941877511129657,0.06125700895733771 +0.00398644647676892,0.0626500454913798 +0.0040315193628604395,0.06407476078962107 +0.004077101866997118,0.06553187524890638 +0.004123199751191874,0.06702212564842888 +0.0041698188426060115,0.06854626552225317 +0.004216965034285822,0.07010506553984533 +0.004264644285907519,0.0716993138964723 +0.0043128626245305815,0.07332981671092179 +0.004361626145359639,0.07499739843358443 +0.004410941012514944,0.07670290226330109 +0.004460813459811576,0.07844719057352244 +0.004511249791547446,0.08023114534837242 +0.004562256383300213,0.08205566862882767 +0.004613839682733214,0.08392168296818268 +0.004666006210410497,0.08583013189947318 +0.004718762560621078,0.08778198041156736 +0.004772115402212517,0.08977821543765789 +0.004826071479433915,0.09181984635437566 +0.0048806376127884426,0.0939079054913396 +0.004935820699895511,0.09604344865411911 +0.004991627716362686,0.09822755565756541 +0.00504806571666747,0.10046133087148534 +0.005105141835049042,0.10274590377951442 +0.005162863286410086,0.10508242955022239 +0.005221237367228817,0.10747208962065333 +0.005280271456481317,0.10991609229423784 +0.005339973016574302,0.11241567335120199 +0.00540034959428843,0.11497209667382374 +0.005461408821732274,0.11758665488516121 +0.005523158417307099,0.12026067000299251 +0.005585606186682515,0.12299549410738239 +0.005648760023783192,0.12579251002536712 +0.005712627911786714,0.1286531320294065 +0.00577721792413272,0.13157880655247436 +0.00584253822554346,0.13457101291983978 +0.005908597073055873,0.13763126409640872 +0.0059754028170653575,0.14076110745183842 +0.006042963902381329,0.14396212554363128 +0.00611128886929471,0.14723593691553322 +0.0061803863546575025,0.15058419691777192 +0.006250265092974553,0.15400859854278245 +0.0063209339175076705,0.15751087328182017 +0.006392401761392223,0.16109279199923796 +0.006464677658766367,0.1647561658295196 +0.006537770745913028,0.16850284709158433 +0.0066116902624148155,0.17233473022599755 +0.0066864455523219755,0.17625375275200642 +0.006762046065333562,0.18026189624827668 +0.006838501357991957,0.1843611873535239 +0.006915821094890896,0.1885536987919748 +0.006994015049897161,0.19284155042074075 +0.0070730931073860704,0.19722691030187584 +0.007153065263490954,0.20171199579829233 +0.007233941627366748,0.20629907469500725 +0.0073157324224678725,0.2109904663451686 +0.00739844798784057,0.215788542843446 +0.007482098779429839,0.2206957302241365 +0.007566695371401163,0.2257145096890026 +0.0076522484574771685,0.23084741886037446 +0.0077387688522894005,0.236097053065262 +0.00782626749274539,0.24146606664602208 +0.007914755439411161,0.24695717430357822 +0.00800424387790939,0.2525731524690524 +0.008094744120333359,0.2583168407078518 +0.008186267606676892,0.2641911431542273 +0.008278825906280481,0.2701990299812203 +0.008372430719293736,0.2763435389000743 +0.008467093878154391,0.2826277766976923 +0.008562827349084015,0.28905492080593403 +0.008659643233600654,0.2956282209086123 +0.008757553770048554,0.30235100058398034 +0.0088565713351452,0.30922665898522245 +0.008956708445545832,0.316258672558462 +0.009057977759425661,0.3234505967998562 +0.009160392078079959,0.3308060680539276 +0.009263964347542264,0.33832880535110454 +0.009368707660220854,0.34602261228726894 +0.009474635256553754,0.35389137894768374 +0.009581760526682423,0.36193908387242335 +0.00969009701214439,0.37016979606791384 +0.009799658407585014,0.3785876770635866 +0.009910458562488608,0.3871969830163247 +0.01002251148292913,0.3960020668604433 +0.010135831333340659,0.4050073805089316 +0.01025043243830791,0.41421747710418627 +0.01036632928437698,0.4236370133182967 +0.010483536521886569,0.4332707517084768 +0.010602068966819901,0.443123563122432 +0.010721941602677596,0.4532004291621567 +0.010843169582371711,0.463506444700801 +0.010965768230141192,0.4740468204588228 +0.011089753043488984,0.4848268856365776 +0.011215139695141062,0.4958520906086265 +0.01134194403502757,0.5071280096795925 +0.01147018209228641,0.5186603438993328 +0.011599870077289447,0.5304549239460079 +0.011731024383691642,0.5425177130730362 +0.011863661590503345,0.5548548101225182 +0.01199779846418602,0.567472452607973 +0.012133451960771677,0.5803770198666285 +0.012270639228006244,0.5935750362838643 +0.012409377607517195,0.6070731745909596 +0.012549684637005683,0.6208782592354627 +0.012691578052463431,0.6349972698337659 +0.012835075790414747,0.6494373446949061 +0.01298019599018381,0.6642057844302998 +0.013126956996187676,0.6793100556423722 +0.013275377360255126,0.6947577946984526 +0.013425475843971808,0.7105568115887709 +0.013577271421051841,0.7267150938745958 +0.013730783279736253,0.7432408107242438 +0.013886030825218543,0.7601423170405287 +0.014043033682097661,0.7774281576842125 +0.014201811696858717,0.7951070717904184 +0.014362384940381745,0.8131879971851571 +0.014524773710478837,0.8316800749024081 +0.014688998534459955,0.8505926538017683 +0.01485508017172775,0.8699352952939114 +0.01502303961640174,0.8897177781700165 +0.015192898099972124,0.9099501035446663 +0.01536467709398364,0.9306424999075279 +0.015538398312749738,0.9518054282916621 +0.01571408371609746,0.9734495875600667 +0.01589175551214333,0.9955859198086622 +0.016071436160100677,1.01822561590015 +0.016253148373118646,1.0413801211103029 +0.01643691512115333,1.0650611409147883 +0.016622759633871387,1.0892806469020149 +0.016810705403586444,1.114050882822309 +0.01700077618822873,1.1393843707703184 +0.017192996014348295,1.1652939175128318 +0.017387389180152136,1.1917926209584542 +0.017583980258575726,1.2188938767723148 +0.01778279410038923,1.2466113851454457 +0.017983855837338837,1.2749591577117834 +0.01818719088532365,1.3039515246275923 +0.01839282494760845,1.3336031418097907 +0.01860078401807282,1.3639289983361444 +0.01881109438449698,1.3949444240231388 +0.019023782631884784,1.4266650971573398 +0.0192388756458243,1.4591070524218372 +0.019456400615886358,1.492286688992136 +0.01967638503906156,1.526220778816483 +0.019898856723236116,1.5609264750890501 +0.02012384379070701,1.5964213209095288 +0.02035137468173687,1.632723258143621 +0.02058147815814909,1.6698506364809234 +0.020814183306963545,1.7078222227030473 +0.02104951954407344,1.7466572101563391 +0.021287516617963728,1.786375228450977 +0.021528204613471574,1.8269963533602427 +0.021771613955589325,1.8685411169626884 +0.022017775413310486,1.9110305180091196 +0.02226672010351919,1.954486032520061 +0.022518479494923624,1.9989296246318653 +0.022773085412033933,2.0443837576812114 +0.023030570039185138,2.0908714055427713 +0.02329096592460546,2.138416064230942 +0.023554305984530732,2.1870417637567545 +0.02382062350736527,2.236773080250751 +0.024089952157889807,2.2876351483826487 +0.02436232598151701,2.3396536740286904 +0.024637779408595104,2.392854947252474 +0.02491634725876017,2.447265855574054 +0.025198064745337617,2.50291389753331 +0.025482967479793468,2.559827196569346 +0.025771091476235932,2.6180345152088003 +0.026062473155967904,2.677565269578135 +0.026357149352090915,2.7384495442461136 +0.02665515731416115,2.800718107397186 +0.026956534712898113,2.864402426359779 +0.027261319644946505,2.929534683475606 +0.02756955063769198,2.9961477923316915 +0.027881266654131334,3.0642754143656807 +0.028196507097797724,3.133951975839494 +0.028515311817741657,3.2052126852048968 +0.028837721113568193,3.2780935508561635 +0.029163775740531202,3.352631399292221 +0.029493516914685138,3.4288638936792037 +0.029826986318095113,3.5068295528485196 +0.030164226104105858,3.586567770709187 +0.030505278902670258,3.6681188361143233 +0.030850187825738135,3.751523953167715 +0.03119899647270598,3.8368252619936496 +0.03155174893592826,3.9240658599801552 +0.03190848980629108,4.013289823495792 +0.03226926417884884,4.1045422301028145 +0.032634117658524635,4.197869181276309 +0.03300309636587508,4.293317825630992 +0.033376246942920386,4.390936382659824 +0.03375361655904026,4.490774167080356 +0.034135252916936525,4.592881613622436 +0.03452120425866314,4.697310302456802 +0.034911519371724424,4.804112985182203 +0.03530624759524219,4.913343611349592 +0.03570543882619261,5.025057355680879 +0.036109143525713656,5.139310645823165 +0.036517412725483776,5.256161190777282 +0.036930298034172734,5.375668009935347 +0.03734785164396542,5.497891462725636 +0.0377701263371593,5.622893279279642 +0.03819717549283665,5.750736591070081 +0.03862905309361203,5.8814859629653755 +0.0390658137324562,6.0152074256474775 +0.03950751261959709,6.151968508832369 +0.039954205589498866,6.291838275158906 +0.040405949107919885,6.43488735524778 +0.04086280027905041,6.581187982756794 +0.04132481685273109,6.730814031004234 +0.041792057231753,6.883841050049991 +0.04226458047924028,7.0403463045825925 +0.04274244632611621,7.200408813116686 +0.04322571517865362,7.364109387177419 +0.0437144481261109,7.531530672264724 +0.044208706948454066,7.702757189291013 +0.044708554124166335,7.877875376948491 +0.04521405283814592,8.056973635529499 +0.04572526698969311,8.2401423707297 +0.046242261200587734,8.427474039437895 +0.046765100823257834,8.619063196067165 +0.047293851949040816,8.81500654001879 +0.04782858141653792,9.015402964225796 +0.048369356820063186,9.220353604775314 +0.04891624651818798,9.429961891656903 +0.04946931964238205,9.644333600655864 +0.050028646105752334,9.863576906406548 +0.050594296611880585,10.08780243663343 +0.05116634266376091,10.31712332763559 +0.05174485657283831,10.551655280991394 +0.05232991146814947,10.791516621550155 +0.05292158130556693,11.03682835672753 +0.05351994087714765,11.28771423712981 +0.05412506582058745,11.544300818544468 +0.05473703262878217,11.806717525302698 +0.05535591865949707,12.075096715095247 +0.05598180214514549,12.349573745222264 +0.05661476220267806,12.630287040326593 +0.057254878843583795,12.917378161656448 +0.05790223298400419,13.210991877872107 +0.0585569064549617,13.51127623744665 +0.059218982012703925,13.818382642679962 +0.059888543349164616,14.1324659253786 +0.060565675102543085,14.453684424217697 +0.0612504628680032,14.78220006384537 +0.061942993208493315,15.118178435749028 +0.06264335366568856,15.46178888093403 +0.06335163277105683,15.813204574447987 +0.06406792005704996,16.1726026117863 +0.06479230606842132,16.540164097238314 +0.06552488237367148,16.916074234200735 +0.06626574157662321,17.300522417502027 +0.0670149773281274,17.693702327800775 +0.06777268433790136,18.095812028045493 +0.06853895838650083,18.507054062148736 +0.0693138963374275,18.92763555580589 +0.07009759614937344,19.357768319581595 +0.07089015688860377,19.79766895425759 +0.07169167874147957,20.24755895860635 +0.07250226302712226,20.707664839430716 +0.0733220122102212,21.178218224171545 +0.07415102991398602,21.659455975938137 +0.07498942093324559,22.151620311322173 +0.07583729124769485,22.6549589197953 +0.07669474803529164,23.169725087808786 +0.0775618996858048,23.696177823802316 +0.07843885581451564,24.234581985901887 +0.07932572727607415,24.785208414684018 +0.08022262617851206,25.348334065503863 +0.08112966589741413,25.924242147149023 +0.0820469610902499,26.51322226132492 +0.08297462771086728,27.11557054503422 +0.08391278302415006,27.731589818441403 +0.08486154562084132,28.36158973347134 +0.08582103543253415,29.005886926906104 +0.0867913737468321,29.664805176607384 +0.08777268322268092,30.338675561208884 +0.08876508790587373,31.02783662326029 +0.08976871324473143,31.732634535994144 +0.09078368610596035,32.4534232736926 +0.09181013479068942,33.190564785882785 +0.0928481890506884,33.94442917523257 +0.09389798010476962,34.71539487940308 +0.09495964065537518,35.503848856847206 +0.09603330490535164,36.31018677669893 +0.09711910857491439,37.134813212719244 +0.09821718891880378,37.97814184159784 +0.0993276847436354,38.840595645394146 +0.10045073642544625,39.72260711856362 +0.1015864859274396,40.624618479318926 +0.10273507681793025,41.54708188569896 +0.10389665428849279,42.49045965617278 +0.10507136517231502,43.455224495228975 +0.10625935796275901,44.4418597235999 +0.10746078283213174,45.450859513589 +0.10867579165066832,46.482729129537866 +0.10990453800572951,47.537985173221074 +0.1111471772212166,48.61715583478335 +0.11240386637720558,49.72078114903211 +0.11367476432980335,50.84941325686784 +0.1149600317312286,52.00361667291192 +0.11625983105011949,53.18396855821043 +0.1175743265920711,54.39105899940754 +0.11890368452040502,55.625491293170334 +0.12024807287717386,56.88788223740512 +0.12160766160440309,58.178862427949284 +0.1229826225655732,59.49907656206876 +0.1243731295673447,60.84918374814943 +0.1257793583815287,62.22985782197331 +0.12720148676730605,63.64178766967647 +0.12863969449369744,65.08567755746047 +0.13009416336228788,66.56224746830246 +0.13156507723020783,68.0722334456942 +0.1330526220333744,69.61638794450025 +0.1345569858099951,71.1954801892934 +0.13607835872433757,72.8102965399747 +0.13761693309076786,74.46164086520035 +0.1391729033980607,76.15033492333254 +0.14074646633398435,77.8772187514437 +0.1423378208101637,79.64315106229968 +0.14394716798722443,81.44900964947664 +0.1455747113002213,83.2956918008446 +0.14722065648435403,85.18411472051072 +0.1488852116009742,87.11521595938123 +0.15056858706388565,89.08995385447032 +0.15227099566594277,91.10930797723698 +0.1539926526059492,93.1742796536724 +0.15573377551586087,95.2858921809212 +0.1574945844882964,97.44519167626221 +0.1592753021043588,99.65324731531277 +0.16107615346177148,101.91115188867751 +0.1628973662033325,104.22002230720389 +0.1647391705456907,106.5810001169915 +0.1666017993084468,108.99525202483112 +0.16848548794358387,111.46397043391697 +0.17039047456523057,113.98837398993426 +0.17231699997976052,116.56970813821106 +0.17426530771623247,119.20924569133122 +0.17623564405717435,121.90828740792728 +0.1782282580697154,124.66816258273909 +0.18024340163707053,127.49022964791251 +0.18228132949038026,130.3758766993324 +0.18434229924091106,133.32652246514024 +0.18642657141261945,136.34361642906705 +0.18853440947508462,139.4286398190531 +0.19066607987681297,142.583106182686 +0.19282185207891958,145.80856206087782 +0.1950019985891904,149.106587673027 +0.1972067949965294,152.47879761401225 +0.19943652000579548,155.9268415667289 +0.20169145547303305,159.45240502431386 +0.20397188644110092,163.0572100280008 +0.20627810117570433,166.7430159165928 +0.20861039120183392,170.51162009178233 +0.21096905134061716,174.3648587915697 +0.2133543797465861,178.30460788424028 +0.2157666779453667,182.33278367184244 +0.21820625087179404,186.4513437083021 +0.220673406908459,190.66228763276624 +0.22316845792468995,194.96765801730965 +0.22569171931597612,199.36954122492534 +0.22824351004383595,203.8700682873857 +0.23082415267613657,208.4714157935875 +0.23343397342786926,213.17580679632738 +0.23607330220238557,217.985511729493 +0.23874247263309994,222.90284934224613 +0.24144182212566392,227.93018764820076 +0.24417169190061688,233.06994488852817 +0.24693242703651927,238.32459051047394 +0.24972437651357351,243.69664616184437 +0.2525478932577379,249.18868669585103 +0.25540333418533956,254.8033411976078 +0.2582910602481916,260.54329402093197 +0.26121143647922046,266.41128584036295 +0.26416483203860924,272.41011472214745 +0.2671516202604625,278.54263720365105 +0.270172178699999,284.81176939292516 +0.2732268891812778,291.22048807645683 +0.276316137845464,297.77183184954146 +0.2794403151996403,304.4689022557203 +0.2825998161661704,311.31486493971767 +0.2857950401326204,318.3129508212102 +0.28902639100224503,325.46645727247 +0.2922942772450439,332.7787493160619 +0.2955991119493963,340.2532608338695 +0.2989413128742783,347.89349578957535 +0.30232130250207156,355.7030294631179 +0.3057395080919683,363.68550969789214 +0.3091963617339809,371.844658159125 +0.3126923004035611,380.184271605166 +0.31622776601683794,388.7082231662653 +0.31980320548647945,397.4204636367971 +0.3234190707781861,406.32502277855457 +0.32707581896782334,415.426010629274 +0.33077391229919956,424.7276188224231 +0.3345138182424978,434.2341219057367 +0.3382960095533678,443.94987869308574 +0.3421209643326863,453.87933358429507 +0.34598916608699326,464.0270179223298 +0.34990110378961076,474.39755133504934 +0.35385727194245375,484.9956430971169 +0.357858170638539,495.82609347735433 +0.3619043056252011,506.8937950988382 +0.3659961883680234,518.203734290536 +0.37013433611549124,529.7609924465771 +0.37431927196437687,541.5707473728531 +0.37855152492586297,553.6382746408877 +0.38283162999241444,565.968948924745 +0.3871601282054056,578.5682453365533 +0.3915375667235127,591.4417407513289 +0.39596449889187924,604.5951151207493 +0.40044148431206356,618.0341527744071 +0.404969088912777,631.7647437054735 +0.4095478850214223,645.7928848390254 +0.4141784514364405,660.1246812812365 +0.4188613735004758,674.7663475469451 +0.4235972431743681,689.7242087619126 +0.4283866591119816,705.0047018381357 +0.43323022673587985,720.6143766203937 +0.4381285583138562,736.5598969977209 +0.44308227303632963,752.8480419805134 +0.44809199709461556,769.4857067377299 +0.4531583637600818,786.4799035920765 +0.45828201346419944,803.8377629671871 +0.46346359387949865,821.5665342867777 +0.4687037600014401,839.6735868183736 +0.4740031742312117,858.1664104597639 +0.479362506459462,877.052616461772 +0.48478243415097966,896.3399380846 +0.49026364243033105,916.0362311842742 +0.49580682416846555,936.149474719944 +0.5014126800703004,956.6877711781124 +0.5070819187632956,977.6593469134292 +0.5128152568870303,999.0725523907511 +0.5186134191837928,1020.93586233965 +0.5244771385901927,1043.2578757904319 +0.530407156329812,1066.0473160054921 +0.5364042220069004,1089.313030292476 +0.5424690937011326,1113.0639896934374 +0.5486025380634357,1137.3092885378846 +0.5548053304129003,1162.0581438597424 +0.5610782548347872,1187.3198946646455 +0.5674221042796429,1213.1040010427703 +0.5738376806635346,1239.4200431172262 +0.5803257949694197,1266.2777198215342 +0.5868872673496606,1293.6868474787902 +0.5935229272296987,1321.657358225559 +0.6002336134129013,1350.1992982098159 +0.6070201741865929,1379.3228255654042 +0.6138834674292863,1409.0382081870837 +0.6208243607191253,1439.355821258463 +0.6278437314435541,1470.2861445454557 +0.6349424669102264,1501.8397594579603 +0.6421214644591686,1534.027345784402 +0.6493816315762113,1566.859678192104 +0.6567238860077028,1600.3476224215365 +0.6641491558765202,1634.5021311801875 +0.6716583797993921,1669.3342397057952 +0.6792525070055475,1704.8550602216683 +0.6869324974567063,1741.0757791482915 +0.6946993219684265,1778.0076501484891 +0.702553962332824,1815.6619866309318 +0.7104974114426788,1854.050157889612 +0.7185306734169451,1893.1835816554 +0.7266547637276809,1933.073716930645 +0.7348707093284117,1973.7320566667806 +0.7431795487839462,2015.1701196282527 +0.7515823324016598,2057.3994434433257 +0.7600801223642624,2100.4315697998913 +0.7686739928640667,2144.2780437235942 +0.7773650302387758,2188.950397773631 +0.7861543331088052,2234.460141346315 +0.7950430125161575,2280.8187558951304 +0.804032192064868,2328.0376737376628 +0.813123008063037,2376.128272036117 +0.8223166096664691,2425.1018583422087 +0.8316141590239368,2474.9696568461977 +0.8410168314240846,2525.7427942335116 +0.8505258154439962,2577.4322852225164 +0.8601423130994411,2630.0490154469253 +0.8698675399968183,2683.60372653671 +0.8797027254868204,2738.106998435166 +0.889649112819833,2793.5692314199982 +0.8997079593030929,2850.0006274645334 +0.9098805364596212,2907.4111706135654 +0.9201681301889558,2965.810606698683 +0.930572040929699,3025.2084222338713 +0.9410935838239043,3085.613822165873 +0.9517340888833214,3147.0357069663896 +0.9624949011575211,3209.4826487460505 +0.9733773809039202,3272.962866461361 +0.9843829037597305,3337.4842001165666 +0.9955128609158502,3403.0540840880662 +1.0067686592927223,3469.6795194472247 +1.018151721718182,3537.3670427034526 +1.029663487107312,3606.1227063705933 +1.041305410644337,3675.952032968999 +1.0530789639665672,3746.859993501296 +1.064985635350429,3818.8509721808186 +1.0770269298995938,3891.9287326652125 +1.0892043697352367,3966.096383229993 +1.101519494188445,4041.356340928523 +1.1139738599948024,4117.7102946946625 +1.1265690414911742,4195.159167404749 +1.1393066308147166,4273.703076971314 +1.1521882381041357,4353.341296314289 +1.1652154917032231,4434.072212200436 +1.1783900383666923,4515.893283503443 +1.191713543468342,4598.800998227781 +1.2051876912115738,4682.79082931072 +1.2188141848422898,4767.8571899772005 +1.2325947468641965,4853.993387924751 +1.2465311192565447,4941.191578861604 +1.2606250636943297,5029.44271820729 +1.2748783617709827,5118.736514753308 +1.2892928152235779,5209.061379999979 +1.3038702461605882,5300.404378465624 +1.3186124972922157,5392.751178980767 +1.333521432163324,5486.086002788149 +1.3485989353890075,5580.391572675205 +1.3638469128928223,5675.6490613695605 +1.3792672921477107,5771.838039642387 +1.3948620224196497,5868.936424323177 +1.410633075014056,5966.920426152955 +1.4265824435249745,6065.76449766282 +1.4427121440870851,6165.4412812314395 +1.4590242156305606,6265.921557519435 +1.4755207201388032,6367.17419389336 +1.492203742909097,6469.166094568243 +1.5090753928162084,6571.8621489548905 +1.526137802578963,6675.22518565563 +1.5433931290298422,6779.215920806193 +1.5608435533876228,6883.7929135931445 +1.5784912815331027,6988.912520423452 +1.5963385442879423,7094.528851228047 +1.614387597696659,7200.5937276525165 +1.6326407233118114,7307.056643386439 +1.6511002284824052,7413.864726627873 +1.669768446645562,7520.962705331083 +1.6886477376214868,7628.29287500271 +1.707740487911767,7735.795069862322 +1.7270491110010484,7843.406637195775 +1.7465760476621182,7951.062415168409 +1.7663237662644407,8058.69471474884 +1.786294763086179,8166.2333056862 +1.8064915626297464,8273.605406971034 +1.8269167179409243,8380.735681796827 +1.8475728109315888,8487.546239759036 +1.868462452706086,8593.956639681877 +1.889588283891298,8699.883903196098 +1.9109529749704406,8805.242531725165 +1.9325592266206335,8909.94453038633 +1.9544097700542906,9013.899438774339 +1.9765073673643667,9117.01436833233 +1.9988548118735103,9219.194047340949 +2.021454928487163,9320.340873492582 +2.0443105740506504,9420.354974436368 +2.0674246377103134,9519.134276526182 +2.090800041278718,9616.574581714473 +2.1144397396040002,9712.569653921171 +2.1383467209433813,9807.011310868249 +2.1625240073409087,9899.789540788279 +2.1869746550094704,9990.792596559966 +2.2117017547171223,10079.907136903365 +2.2367084321777915,10167.0183488794 +2.2619978484463887,10252.010105359164 +2.2875732003183957,10334.76508603905 +2.3134377207339654,10415.164981306672 +2.339594679186593,10493.09065008253 +2.366047382136408,10568.422301370194 +2.3927991734281377,10641.039692244196 +2.419853434713799,10710.822333836339 +2.4472135858801662,10777.649705958323 +2.474883085481074,10841.401484215481 +2.502865431174608,10901.9577746117 +2.5311641601652384,10959.199359405222 +2.5597828496509516,11013.007952712001 +2.58872511727544,11063.26646158287 +2.6179946215854,11109.859260212315 +2.647595062493006,11152.672469674644 +2.6775301817436077,11191.594245727054 +2.707803763388721,11226.515073765755 +2.7384196342643614,11257.32806985329 +2.769381664474791,11283.929287287441 +2.80069376788173,11306.218028543923 +2.8323599025991038,11324.097160335106 +2.8643840714933795,11337.473432667332 +2.896770322689565,11346.25779859818 +2.9295227500829233,11350.365736850243 +2.9626454938564777,11349.717572694115 +2.996142741004364,11344.238798438244 +3.030018725861103,11333.86039151568 +3.064277730636856,11318.519128642434 +3.0989240859587324,11298.157895041708 +3.133962171418216,11272.725987474802 +3.169396416124784,11242.17940916977 +3.205231299265784,11206.481155657784 +3.241471350672639,11165.601492320204 +3.2781211513934587,11119.518214883577 +3.3151853342721207,11068.216899633375 +3.3526685845339017,11011.69114102927 +3.390575640377731,10949.942767884282 +3.4289112935751356,10882.982043902348 +3.4676803900759636,10810.827849264972 +3.506887830620951,10733.507841160606 +3.5465385713612183,10651.058589922128 +3.586637624484769,10563.525696192408 +3.627190058850071,10470.96387902024 +3.668201000626807,10373.437039365308 +3.7096756339438612,10271.018297236105 +3.751619201544639,10163.789999553941 +3.7940370054497943,10051.843699269126 +3.836934407627449,9935.280104590856 +3.880316830670991,9814.20899789238 +3.9241897584845358,9688.749121794459 +3.968558736976138,9559.028035430607 +4.013429374758842,9425.181937269035 +4.058807343859655,9287.355456581148 +4.1046983804365444,9145.701412264121 +4.151108285503529,9000.38053727824 +4.1980429256639855,8851.561181326855 +4.2455082338522265,8699.418963291864 +4.293510210083482,8544.136411721733 +4.342054922212346,8385.902566247487 +4.391148506699809,8224.912552111602 +4.440797169388953,8061.367125910973 +4.49100718628943,7895.4722010264095 +4.541784904370795,7727.438335554626 +4.593136742364821,7557.480212355696 +4.645069191576877,7385.816088796023 +4.697588816706491,7212.667230570525 +4.750702256677176,7038.2573249931565 +4.804416225475647,6862.811885108641 +4.858737513000528,6686.557642333076 +4.913672985920654,6509.7219199988685 +4.96922958854307,6332.5320104685115 +5.025414343690856,6155.214542930454 +5.082234353590866,5977.99485276087 +5.139696800771513,5801.096343915776 +5.1978089489707004,5624.739874008648 +5.25657814405402,5449.143125353289 +5.316011814943327,5274.51999775157 +5.376117474555825,5101.080015838645 +5.436902720753759,4929.0277463848615 +5.498375237304849,4758.562237822614 +5.560542794853582,4589.876480047782 +5.623413251903491,4423.156889787381 +5.68699455581053,4258.582819803359 +5.751294743787694,4096.326100931488 +5.816321943920984,3936.550611786819 +5.88208437619687,3779.4118837499545 +5.9485903535413645,3625.056740391685 +6.015848282870847,3473.6229729741876 +6.083866666154767,3325.239058685486 +6.152654101490373,3180.0239097173594 +6.22221928418957,3038.0866697206034 +6.292571007878093,2899.526556425063 +6.363718165607093,2764.4327286301336 +6.435669750977286,2632.8842139201947 +6.508434859275831,2504.9498692261604 +6.582022688626042,2380.6883852722012 +6.6564425411501125,2260.1483314755815 +6.731703824144982,2143.3682470090107 +6.807816051271499,2030.3767609997396 +6.884788843757024,1921.1927616786966 +6.962631931611635,1815.825598357301 +7.0413551548580875,1714.275319290492 +7.120968464775669,1616.5329433594188 +7.201481925158132,1522.580763615521 +7.282905713585835,1432.392680363651 +7.36525012271228,1345.934561467306 +7.44852556156519,1263.1646268705663 +7.532742556862294,1184.0338560869243 +7.617911754341997,1108.4864143476323 +7.704043920109092,1036.460095580982 +7.7911499419956805,967.8867792344595 +7.8792408309374915,902.6928980756699 +7.968327722365756,840.7999139667121 +8.058421877614817,782.1247989988125 +8.149534685345662,726.5805190019669 +8.241677662985538,674.076516750729 +8.334862458183856,624.5191921747465 +8.42910085028456,577.8123769801585 +8.524404751815114,533.8578012217761 +8.62078620999237,492.55554947738904 +8.718257408245425,453.8045044414813 +8.81683066775571,417.5027758831135 +8.91651844901449,383.5481130846003 +9.017333353397982,351.8382991108456 +9.11928812476027,322.2715253729263 +9.222395651044236,294.746745190977 +9.326668965910704,269.164005256279 +9.432121250386007,245.4247540020454 +9.538765834528181,223.43212632545635 +9.646616199111993,203.09120399622196 +9.755685977333021,184.30925152975448 +9.86598895653102,166.99592734509264 +9.977539079932738,151.06347031432614 +10.090350448414473,136.42686188390914 +10.20443732228454,123.0039641755934 +10.319814123085884,110.71563483537534 +10.436495435419097,99.48581871016837 +10.55449600878603,89.2416179577598 +10.673830759454248,79.91334118596603 +10.794514772342584,71.43453239957388 +10.916563302927996,63.74198132919564 +11.039991779173977,56.77571611166695 +11.164815803480792,50.47897963969098 +11.291051154657756,44.798190888263726 +11.418713789917796,39.682892551436645 +11.547819846894583,35.08568634597809 +11.678385645682464,30.96215733728492 +11.81042769089947,27.270788633045505 +11.943962673773619,23.972867770912927 +12.079007474252844,21.03238608737521 +12.215579163138754,18.41593232936213 +12.353695004244532,16.092581693507917 +12.49337245657722,14.033781445798072 +12.634629176544685,12.213234194519112 +12.77748302018755,10.606779826789804 +12.921952045436331,9.192277033583938 +13.068054514394127,7.949485285642717 +13.215808895645086,6.859948046204727 +13.365233866589012,5.906877894790536 +13.516348315802327,5.075044203371626 +13.669171345425756,4.350663885597188 +13.823722273578996,3.721295693335178 +13.980020636802688,3.175738431490957 +14.138086192528005,2.703933421992661 +14.297938921574154,2.2968714503407193 +14.459599030674115,1.9465043817053924 +14.62308695502896,1.6456615679795483 +14.788423360891013,1.3879711129183148 +14.955629148176254,1.167786013263033 +15.124725453106233,0.9801151468554572 +15.295733650879885,0.8205590419306958 +15.468675358375512,0.6852503239911124 +15.643572436883355,0.5707987082316195 +15.820446994869037,0.47424038023283893 +15.999321390768271,0.39299158680041 +16.180218235813136,0.32480624012393694 +16.36316039689035,0.2677373333132984 +16.548170999431814,0.22010194712151682 +16.73527343033788,0.18044962887088828 +16.924491340933645,0.14753392130480572 +17.11584864995868,0.12028681817363185 +17.309369546590553,0.09779592741374835 +17.505078493502552,0.07928412735350837 +17.703000229955983,0.06409150813797643 +17.903159774927396,0.05165939840260269 +18.10558243027122,0.0415162867658424 +18.31029378391811,0.033265457356700626 +18.517319713109497,0.02657416938173453 +18.72668638766867,0.021164221860646266 +18.938420273308882,0.01680375549762264 +19.152548134978822,0.01330015521162039 +19.369097040245936,0.010493927669105571 +19.588094362718003,0.008253439023184931 +19.80956778550339,0.0064704084617957566 +20.0335453047104,0.005056063124834915 +20.26005523298627,0.003937869388959245 +20.48912620309608,0.0030567644005791844 +20.720787171542206,0.0023648200143755417 +20.95506742222465,0.0018232789801529065 +21.191996570142773,0.0014009102718563194 +21.43160456513889,0.0010726369182727984 +21.673921695684175,0.0008183955536283347 +21.918978592707386,0.0006221921985326174 +22.166806233466865,0.0004713235334323864 +22.417435945466323,0.0003557371632839696 +22.67089941041491,0.0002675081297839018 +22.927228668232054,0.00020041224382060625 +23.186456121097557,0.0001495797185595647 +23.44861453754752,0.0001112151213547045 +23.71373705661655,0.00008237186460311056 +23.981857192026837,0.000060771356721383726 +24.25300883642454,0.000044658566744002184 +24.527226265664133,0.00003268715045555007 +24.804544143141136,0.000023828470964264034 +25.084997524173872,0.00001729984844243988 +25.368621860434764,0.000012508216381812409 +25.655453004431713,9.006066798974251e-6 +25.945527214040155,6.457153809183881e-6 +26.238881157086375,4.609911142604508e-6 +26.535551915982616,3.2769398102395563e-6 +26.835576992414623,2.319250491646976e-6 +27.13899431208216,1.634213096143189e-6 +27.445842229493145,1.1463833360533437e-6 +27.756159532811978,8.005514713121728e-7 +28.06998544876269,5.564993888414769e-7 +28.387359647587548,3.8506471851030537e-7 +28.708322248061688,2.652001006836486e-7 +29.032913822564485,1.8178649626713763e-7 +29.36117540220822,1.2401503240500314e-7 +29.69314848202479,8.419539824638761e-8 +30.028875026210994,5.6882668990478024e-8 +30.36839747343319,3.8240641251683425e-8 +30.7117587421919,2.5579949786854332e-8 +31.05900223624705,1.7024681165499017e-8 +31.410171850104575,1.1272967250602148e-8 +31.76531197456508,7.4259558211082825e-9 +32.12446750233517,4.866267185620201e-9 +32.48768383370232,3.172074231152613e-9 +32.8550068822738,2.0566887562004345e-9 +33.226483080780575,1.3263110597862835e-9 +33.60215938694678,8.506428224055657e-10 +33.982083289425596,5.425593770986373e-10 +34.36630281380217,3.441266694671701e-10 +34.7548665286645,2.1703658708078903e-10 +35.1478235517429,1.361016408422853e-10 +35.54522355611888,8.485578839091784e-11 +35.947116776504245,5.259674664408699e-11 +36.35355401559122,3.240909188544561e-11 +36.76458665047429,1.9850724540301276e-11 +37.18026663914475,1.2085315843091875e-11 +37.60064652705856,7.312772737494184e-12 +38.02577945377861,4.397628428654792e-12 +38.45571915969198,2.6280703194226e-12 +38.89051999280318,1.5606542393748956e-12 +39.33023691560415,9.208663390847817e-13 +39.774925512022065,5.398526749209212e-13 +40.22464199444556,3.144205258495473e-13 +40.67944321083047,1.8191618947044476e-13 +41.13938665188586,1.0455002290389658e-13 +41.60453045834138,5.968103693225428e-14 +42.07493342829669,3.383567053025836e-14 +42.550655024654105,1.9050474235942726e-14 +43.031755382635154,1.0651122250425509e-14 +43.51829531738218,5.91301708012776e-15 +44.01033633164593,3.259209733372893e-15 +44.50794062355996,1.783486085209242e-15 +45.011171094503055,9.68826070095787e-16 +45.520091357050475,5.224017388786075e-16 +46.03476574301511,2.7958234373240063e-16 +46.55525931157959,1.4849930189956814e-16 +47.08163785752028,7.827289995728395e-17 +47.61396791952432,4.09386157625223e-17 +48.15231678860069,2.1244746397215882e-17 +48.69675251658631,1.0937753661244456e-17 +49.24734392474841,5.58629672922666e-18 +49.80416061248411,2.830086989053673e-18 +50.36727296611836,1.4220527021548341e-18 +50.936752167801345,7.086494805974369e-19 +51.5126702045066,3.501921645231485e-19 +52.095099877130636,1.7159292693049547e-19 +52.6841148096957,8.336211342426702e-20 +53.27978945865641,4.0148777640431573e-20 +53.88219912231171,1.9167595640011264e-20 +54.49141995032318,9.070089937180992e-21 +55.10752895334102,4.253639172366972e-21 +55.73060401273886,1.9768365508658653e-21 +56.36072389045857,9.103269619070214e-22 +56.99796823896648,4.1533213928538814e-22 +57.64241761132211,1.877239202353527e-22 +58.29415347136074,8.40472623157247e-23 +58.953258203991155,3.727012300917363e-23 +59.619815125609776,1.6367572334417447e-23 +60.29390849463255,7.117807818983005e-24 +60.975623522145916,3.0647751121077248e-24 +61.665046382678256,1.3064494385720478e-24 +62.36226422509303,5.51288539881079e-25 +63.06736518360512,2.302543547013738e-25 +63.78043838892178,9.517631235742949e-26 +64.5015739795095,3.8930596832530986e-26 +65.23086311298825,1.5755887581505155e-26 +65.96839797765456,6.308588257447771e-27 +66.71427180413495,2.4986548899791622e-27 +67.46857887717101,9.78840300237391e-28 +68.23141454753784,3.792232584771802e-28 +69.00287524409713,1.4527842550307246e-28 +69.78305848598663,5.502691683806356e-29 +70.5720628949474,2.060446620092143e-29 +71.36998820779036,7.626114009857827e-30 +72.17693528900395,2.789612822484186e-30 +72.99300614350419,1.0083831988134831e-30 +73.81830392952901,3.6015543726947214e-31 +74.65293297167825,1.2708025220293212e-31 +75.49699877410127,4.4292503092625953e-32 +76.35060803383345,1.5247060212512635e-32 +77.21386865428371,5.1830420696277354e-33 +78.08688975887432,1.7396599830046824e-33 +78.96978170483506,5.764513218831596e-34 +79.8626560971533,1.8854561130293748e-34 +80.76562580268183,6.086420062020301e-35 +81.67880496440614,1.938807023703275e-35 +82.60230901587306,6.093532694328715e-36 +83.53625469578262,1.889294244826269e-36 +84.48076006274468,5.7777560259614265e-37 +85.43594451020262,1.74253067400805e-37 +86.40192878152561,5.181971609711881e-38 +87.37883498527171,1.5192664596057057e-38 +88.3667866106233,4.3906353736876636e-39 +89.36590854299715,1.2505604314756292e-39 +90.3763270798311,3.509897821970834e-40 +91.39816994654906,9.705647771884819e-41 +92.43156631270651,2.6437586583421683e-41 +93.47664680831878,7.09271413759933e-42 +94.53354354037361,1.8737937259146197e-42 +95.60239010953076,4.8738733507896606e-43 +96.68332162701009,1.2479392421317118e-43 +97.7764747316709,3.144865884629114e-44 +98.88198760728412,7.798699260546751e-45 +100.,1.902718581787403e-45 \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e6930ae0b..f51a84c87 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -166,6 +166,7 @@ add_subdirectory(RadForce) add_subdirectory(RadMarshak) add_subdirectory(RadMarshakCGS) add_subdirectory(RadMarshakAsymptotic) +add_subdirectory(RadMarshakVaytet) add_subdirectory(RadMatterCoupling) add_subdirectory(RadMatterCouplingRSLA) add_subdirectory(RadPulse) @@ -184,8 +185,9 @@ add_subdirectory(RadhydroUniformAdvecting) add_subdirectory(RadhydroPulse) add_subdirectory(RadhydroPulseDyn) add_subdirectory(RadhydroPulseGrey) -add_subdirectory(RadhydroPulseMG) +add_subdirectory(RadhydroPulseMGconst) add_subdirectory(RadhydroPulseMGint) +add_subdirectory(RadhydroBB) add_subdirectory(BinaryOrbitCIC) add_subdirectory(Cooling) diff --git a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index 4d21a3f6d..e4ad425c7 100644 --- a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -71,18 +71,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::Comp return ComputePlanckOpacity(rho, Tgas); } -template <> -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto -RadSystem::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) -> quokka::valarray -{ - quokka::valarray opacity_deriv{}; - auto sigma_dT = (-3.0 * kappa / Tgas) * std::pow(Tgas / T_hohlraum, -3); // cm^-1 - for (int i = 0; i < nGroups_; ++i) { - opacity_deriv[i] = sigma_dT / rho; - } - return opacity_deriv; -} - template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double { return (1. / 3.); // Eddington approximation diff --git a/src/RadMarshakVaytet/CMakeLists.txt b/src/RadMarshakVaytet/CMakeLists.txt new file mode 100644 index 000000000..f92e06204 --- /dev/null +++ b/src/RadMarshakVaytet/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_radiation_marshak_Vaytet test_radiation_marshak_Vaytet.cpp ../fextract.cpp ../interpolate.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radiation_marshak_Vaytet) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME MarshakWaveVaytet COMMAND test_radiation_marshak_Vaytet MarshakVaytet.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp b/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp new file mode 100644 index 000000000..9b82db9df --- /dev/null +++ b/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp @@ -0,0 +1,486 @@ +/// \file test_radiation_marshak_vaytet.cpp +/// \brief Defines a Marshak wave problem with variable opacity. +/// + +#include "AMReX_BLassert.H" + +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "radiation_system.hpp" +#include "test_radiation_marshak_Vaytet.hpp" + +// constexpr int n_groups_ = 2; // Be careful +constexpr int n_groups_ = 4; +// constexpr int n_groups_ = 8; +// constexpr int n_groups_ = 16; +// constexpr int n_groups_ = 64; +// constexpr int n_groups_ = 128; +// constexpr int n_groups_ = 256; +// constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity; +// constexpr OpacityModel opacity_model_ = OpacityModel::PPL_opacity_fixed_slope_spectrum; +constexpr OpacityModel opacity_model_ = OpacityModel::PPL_opacity_full_spectrum; + +constexpr double kappa0 = 2000.0; // cm^2 g^-1 (opacity). +constexpr double nu_pivot = 4.0e13; // Powerlaw, kappa = kappa0 (nu/nu_pivot)^{-2} +constexpr int n_coll = 4; // number of collections = 6, to align with Vaytet +constexpr int the_model = 10; // 0: constant opacity (Vaytet et al. Sec 3.2.1), 1: nu-dependent opacity (Vaytet et al. Sec 3.2.2), 2: nu-and-T-dependent opacity + // (Vaytet et al. Sec 3.2.3) +// 10: bin-centered method with opacities propto nu^-2 + +// OLD +// constexpr int the_model = 0; // 0: constant opacity (Vaytet et al. Sec 3.2.1), 1: nu-dependent opacity (Vaytet et al. Sec 3.2.2), 2: nu-and-T-dependent +// opacity (Vaytet et al. Sec 3.2.3) constexpr int n_groups_ = 6; constexpr amrex::GpuArray group_edges_ = {0.3e12, 0.3e14, 0.6e14, +// 0.9e14, 1.2e14, 1.5e14, 1.5e16}; constexpr amrex::GpuArray group_opacities_ = {1000., 750., 500., 250., 10., 10.}; + +// NEW +constexpr amrex::GpuArray group_edges_ = []() constexpr { + // bins cover four orders of magnitutde from 6.0e10 to 6.0e14, with N bins logarithmically spaced + // in the space of x = h nu / (k_B T) where T = 1000 K, this roughly corresponds to x = 3e-3 to 3e1 + // n_groups_ = 2, 4, 8, 16, 128, 256 + if constexpr (n_groups_ == 2) { + return amrex::GpuArray{6.0e10, 6.0e12, 6.0e14}; + } else if constexpr (n_groups_ == 4) { + return amrex::GpuArray{6.0e10, 6.0e11, 6.0e12, 6.0e13, 6.0e14}; + } else if constexpr (n_groups_ == 8) { + return amrex::GpuArray{6.0000000e+10, 1.8973666e+11, 6.0000000e+11, 1.8973666e+12, 6.0000000e+12, + 1.8973666e+13, 6.0000000e+13, 1.8973666e+14, 6.0000000e+14}; + } else if constexpr (n_groups_ == 16) { + return amrex::GpuArray{6.00000000e+10, 1.06696765e+11, 1.89736660e+11, 3.37404795e+11, 6.00000000e+11, 1.06696765e+12, + 1.89736660e+12, 3.37404795e+12, 6.00000000e+12, 1.06696765e+13, 1.89736660e+13, 3.37404795e+13, + 6.00000000e+13, 1.06696765e+14, 1.89736660e+14, 3.37404795e+14, 6.00000000e+14}; + } else if constexpr (n_groups_ == 64) { + return amrex::GpuArray{ + 6.00000000e+10, 6.92869191e+10, 8.00112859e+10, 9.23955916e+10, 1.06696765e+11, 1.23211502e+11, 1.42282422e+11, 1.64305178e+11, + 1.89736660e+11, 2.19104476e+11, 2.53017902e+11, 2.92180515e+11, 3.37404795e+11, 3.89628979e+11, 4.49936526e+11, 5.19578594e+11, + 6.00000000e+11, 6.92869191e+11, 8.00112859e+11, 9.23955916e+11, 1.06696765e+12, 1.23211502e+12, 1.42282422e+12, 1.64305178e+12, + 1.89736660e+12, 2.19104476e+12, 2.53017902e+12, 2.92180515e+12, 3.37404795e+12, 3.89628979e+12, 4.49936526e+12, 5.19578594e+12, + 6.00000000e+12, 6.92869191e+12, 8.00112859e+12, 9.23955916e+12, 1.06696765e+13, 1.23211502e+13, 1.42282422e+13, 1.64305178e+13, + 1.89736660e+13, 2.19104476e+13, 2.53017902e+13, 2.92180515e+13, 3.37404795e+13, 3.89628979e+13, 4.49936526e+13, 5.19578594e+13, + 6.00000000e+13, 6.92869191e+13, 8.00112859e+13, 9.23955916e+13, 1.06696765e+14, 1.23211502e+14, 1.42282422e+14, 1.64305178e+14, + 1.89736660e+14, 2.19104476e+14, 2.53017902e+14, 2.92180515e+14, 3.37404795e+14, 3.89628979e+14, 4.49936526e+14, 5.19578594e+14, + 6.00000000e+14}; + } else if constexpr (n_groups_ == 128) { + return amrex::GpuArray{ + 6.00000000e+10, 6.44764697e+10, 6.92869191e+10, 7.44562656e+10, 8.00112859e+10, 8.59807542e+10, 9.23955916e+10, 9.92890260e+10, + 1.06696765e+11, 1.14657178e+11, 1.23211502e+11, 1.32404044e+11, 1.42282422e+11, 1.52897805e+11, 1.64305178e+11, 1.76563631e+11, + 1.89736660e+11, 2.03892500e+11, 2.19104476e+11, 2.35451386e+11, 2.53017902e+11, 2.71895018e+11, 2.92180515e+11, 3.13979469e+11, + 3.37404795e+11, 3.62577834e+11, 3.89628979e+11, 4.18698351e+11, 4.49936526e+11, 4.83505313e+11, 5.19578594e+11, 5.58343225e+11, + 6.00000000e+11, 6.44764697e+11, 6.92869191e+11, 7.44562656e+11, 8.00112859e+11, 8.59807542e+11, 9.23955916e+11, 9.92890260e+11, + 1.06696765e+12, 1.14657178e+12, 1.23211502e+12, 1.32404044e+12, 1.42282422e+12, 1.52897805e+12, 1.64305178e+12, 1.76563631e+12, + 1.89736660e+12, 2.03892500e+12, 2.19104476e+12, 2.35451386e+12, 2.53017902e+12, 2.71895018e+12, 2.92180515e+12, 3.13979469e+12, + 3.37404795e+12, 3.62577834e+12, 3.89628979e+12, 4.18698351e+12, 4.49936526e+12, 4.83505313e+12, 5.19578594e+12, 5.58343225e+12, + 6.00000000e+12, 6.44764697e+12, 6.92869191e+12, 7.44562656e+12, 8.00112859e+12, 8.59807542e+12, 9.23955916e+12, 9.92890260e+12, + 1.06696765e+13, 1.14657178e+13, 1.23211502e+13, 1.32404044e+13, 1.42282422e+13, 1.52897805e+13, 1.64305178e+13, 1.76563631e+13, + 1.89736660e+13, 2.03892500e+13, 2.19104476e+13, 2.35451386e+13, 2.53017902e+13, 2.71895018e+13, 2.92180515e+13, 3.13979469e+13, + 3.37404795e+13, 3.62577834e+13, 3.89628979e+13, 4.18698351e+13, 4.49936526e+13, 4.83505313e+13, 5.19578594e+13, 5.58343225e+13, + 6.00000000e+13, 6.44764697e+13, 6.92869191e+13, 7.44562656e+13, 8.00112859e+13, 8.59807542e+13, 9.23955916e+13, 9.92890260e+13, + 1.06696765e+14, 1.14657178e+14, 1.23211502e+14, 1.32404044e+14, 1.42282422e+14, 1.52897805e+14, 1.64305178e+14, 1.76563631e+14, + 1.89736660e+14, 2.03892500e+14, 2.19104476e+14, 2.35451386e+14, 2.53017902e+14, 2.71895018e+14, 2.92180515e+14, 3.13979469e+14, + 3.37404795e+14, 3.62577834e+14, 3.89628979e+14, 4.18698351e+14, 4.49936526e+14, 4.83505313e+14, 5.19578594e+14, 5.58343225e+14, + 6.00000000e+14}; + } +}(); + +constexpr amrex::GpuArray group_opacities_{}; + +struct SuOlsonProblemCgs { +}; // dummy type to allow compile-type polymorphism via template specialization + +constexpr int max_step_ = 1e6; +constexpr double rho0 = 1.0e-3; // g cm^-3 +constexpr double T_initial = 300.0; // K +constexpr double T_L = 1000.0; // K +constexpr double T_R = 300.0; // K +constexpr double rho_C_V = 1.0e-3; // erg cm^-3 K^-1 +constexpr double c_v = rho_C_V / rho0; +constexpr double mu = 1.0 / (5. / 3. - 1.) * C::k_B / c_v; + +constexpr double a_rad = radiation_constant_cgs_; +constexpr double Erad_floor_ = a_rad * T_initial * T_initial * T_initial * T_initial * 1e-20; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = C::k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = false; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_groups_; // number of radiation groups +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c_light_cgs_; + static constexpr double c_hat = c_light_cgs_; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = Erad_floor_; + static constexpr int beta_order = 0; + static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz + static constexpr amrex::GpuArray radBoundaries = group_edges_; + static constexpr OpacityModel opacity_model = opacity_model_; +}; + +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray rad_boundaries, const double /*rho*/, + const double Tgas) -> amrex::GpuArray, 2> +{ + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_ + 1; ++i) { + if constexpr (the_model < 10) { + exponents_and_values[0][i] = 0.0; + } else { + exponents_and_values[0][i] = -2.0; + } + } + if constexpr (the_model == 0) { + for (int i = 0; i < nGroups_; ++i) { + exponents_and_values[1][i] = kappa0; + } + } else if constexpr (the_model == 1) { + for (int i = 0; i < nGroups_; ++i) { + exponents_and_values[1][i] = group_opacities_[i]; + } + } else if constexpr (the_model == 2) { + for (int i = 0; i < nGroups_; ++i) { + exponents_and_values[1][i] = group_opacities_[i] * std::pow(Tgas / T_initial, 3. / 2.); + } + } else if constexpr (the_model == 10) { + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + for (int i = 0; i < nGroups_; ++i) { + auto const bin_center = std::sqrt(rad_boundaries[i] * rad_boundaries[i + 1]); + exponents_and_values[1][i] = kappa0 * std::pow(bin_center / nu_pivot, -2.); + } + } else { + for (int i = 0; i < nGroups_ + 1; ++i) { + exponents_and_values[1][i] = kappa0 * std::pow(rad_boundaries[i] / nu_pivot, -2.); + } + } + } + return exponents_and_values; +} + +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, int /*dcomp*/, + int /*numcomp*/, amrex::GeometryData const &geom, const amrex::Real /*time*/, + const amrex::BCRec * /*bcr*/, int /*bcomp*/, int /*orig_comp*/) +{ +#if (AMREX_SPACEDIM == 1) + auto i = iv.toArray()[0]; + int j = 0; + int k = 0; +#endif +#if (AMREX_SPACEDIM == 2) + auto [i, j] = iv.toArray(); + int k = 0; +#endif +#if (AMREX_SPACEDIM == 3) + auto [i, j, k] = iv.toArray(); +#endif + + amrex::Box const &box = geom.Domain(); + amrex::GpuArray lo = box.loVect3d(); + amrex::GpuArray hi = box.hiVect3d(); + + auto const radBoundaries_g = RadSystem::radBoundaries_; + + if (i < lo[0] || i >= hi[0]) { + double T_H = NAN; + if (i < lo[0]) { + T_H = T_L; + } else { + T_H = T_R; + } + + auto Erad_g = RadSystem::ComputeThermalRadiation(T_H, radBoundaries_g); + const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T_initial); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + consVar(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + consVar(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + consVar(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + } + + // gas boundary conditions are the same on both sides + consVar(i, j, k, RadSystem::gasEnergy_index) = Egas; + consVar(i, j, k, RadSystem::gasDensity_index) = rho0; + consVar(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + consVar(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + consVar(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + } +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const auto radBoundaries_g = RadSystem::radBoundaries_; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T_initial); + // const double Erad = a_rad * std::pow(T_initial, 4); + auto Erad_g = RadSystem::ComputeThermalRadiation(T_initial, radBoundaries_g); + + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +auto problem_main() -> int +{ + // This problem tests whether the M1 closure can handle a Marshak wave problem with variable opacity + // as accurately as ray-tracing method. + + // Problem parameters + const int max_timesteps = max_step_; + const double CFL_number = 0.8; + // const double initial_dt = 5.0e-12; // s + const double max_dt = 1.0; // s + const double max_time = 1.36e-7; + + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + BCs_cc[n].setLo(0, + amrex::BCType::ext_dir); // custom (Marshak) x1 + BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1 + for (int i = 1; i < AMREX_SPACEDIM; ++i) { + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + // sim.initDt_ = initial_dt; + sim.maxDt_ = max_dt; + sim.radiationCflNumber_ = CFL_number; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); + const int nx = static_cast(position.size()); + + // compare against diffusion solution + const int status = 0; + if (amrex::ParallelDescriptor::IOProcessor()) { + std::vector xs(nx); + std::vector Trad(nx); + std::vector Tgas(nx); + std::vector vgas(nx); + // define a vector of n_groups_ vectors + std::vector> Trad_g(n_groups_); + + int const n_item = n_groups_ / n_coll; + std::vector> Trad_coll(n_coll); + + for (int i = 0; i < nx; ++i) { + double Erad_t = 0.; + // const double Erad_t = values.at(RadSystem::radEnergy_index)[i]; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const double Egas_t = values.at(RadSystem::gasInternalEnergy_index)[i]; + const double rho = values.at(RadSystem::gasDensity_index)[i]; + amrex::Real const x = position[i]; + xs.at(i) = x; + Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho, Egas_t); + Trad.at(i) = std::pow(Erad_t / a_rad, 1. / 4.); + + // vgas + const auto x1GasMomentum = values.at(RadSystem::x1GasMomentum_index)[i]; + vgas.at(i) = x1GasMomentum / rho; + + int counter = 0; + int group_counter = 0; + double Erad_sum = 0.; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + auto Erad_g = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + Trad_g[g].push_back(std::pow(Erad_g / a_rad, 1. / 4.)); + + if (counter == 0) { + Erad_sum = 0.0; + } + Erad_sum += Erad_g; + if (counter == n_item - 1) { + Trad_coll[group_counter].push_back(std::pow(Erad_sum / a_rad, 1. / 4.)); + group_counter++; + counter = 0; + } else { + counter++; + } + } + } + + // // read in exact solution + + // std::vector xs_exact; + // std::vector Tmat_exact; + + // std::string filename = "../extern/marshak_similarity.csv"; + // std::ifstream fstream(filename, std::ios::in); + // AMREX_ALWAYS_ASSERT(fstream.is_open()); + + // std::string header; + // std::getline(fstream, header); + + // for (std::string line; std::getline(fstream, line);) { + // std::istringstream iss(line); + // std::vector values; + + // for (double value = NAN; iss >> value;) { + // values.push_back(value); + // } + // auto x_val = values.at(0); + // auto Tmat_val = values.at(1); + + // xs_exact.push_back(x_val); + // Tmat_exact.push_back(Tmat_val); + // } + + // // compute error norm + + // // interpolate numerical solution onto exact tabulated solution + // std::vector Tmat_interp(xs_exact.size()); + // interpolate_arrays(xs_exact.data(), Tmat_interp.data(), static_cast(xs_exact.size()), xs.data(), Tgas.data(), + // static_cast(xs.size())); + + // double err_norm = 0.; + // double sol_norm = 0.; + // for (size_t i = 0; i < xs_exact.size(); ++i) { + // err_norm += std::abs(Tmat_interp[i] - Tmat_exact[i]); + // sol_norm += std::abs(Tmat_exact[i]); + // } + + // const double error_tol = 0.09; + // const double rel_error = err_norm / sol_norm; + // amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + + // save data to file + std::ofstream fstream; + fstream.open("marshak_wave_Vaytet.csv"); + fstream << "# x, Tgas, Trad"; + for (int i = 0; i < n_groups_; ++i) { + fstream << ", " << "Trad_" << i; + } + for (int i = 0; i < nx; ++i) { + fstream << std::endl; + fstream << std::scientific << std::setprecision(14) << xs[i] << ", " << Tgas[i] << ", " << Trad[i]; + for (int j = 0; j < n_groups_; ++j) { + fstream << ", " << Trad_g[j][i]; + } + } + fstream.close(); + + // save Trad_coll to file + std::ofstream fstream_coll; + fstream_coll.open("marshak_wave_Vaytet_coll.csv"); + fstream_coll << "# x, Tgas, Trad"; + for (int i = 0; i < n_coll; ++i) { + fstream_coll << ", " << "Trad_" << i; + } + for (int i = 0; i < nx; ++i) { + fstream_coll << std::endl; + fstream_coll << std::scientific << std::setprecision(14) << xs[i] << ", " << Tgas[i] << ", " << Trad[i]; + for (int j = 0; j < n_coll; ++j) { + fstream_coll << ", " << Trad_coll[j][i]; + } + } + fstream_coll.close(); + + // // check if velocity is strictly zero + // const double error_v_tol = 1.0e-10; + // double error_v = 0.0; + // const double cs = std::sqrt(5. / 3. * C::k_B / mu * T_initial); // sound speed + // for (size_t i = 0; i < xs.size(); ++i) { + // error_v += std::abs(vgas[i]) / cs; + // } + // amrex::Print() << "Sum of abs(v) / cs = " << error_v << std::endl; + // if ((error_v > error_v_tol) || std::isnan(error_v)) { + // status = 1; + // } + +#ifdef HAVE_PYTHON + // plot results + matplotlibcpp::clf(); + std::map args; + args["label"] = "gas"; + args["linestyle"] = "-"; + args["color"] = "k"; + matplotlibcpp::plot(xs, Tgas, args); + args["label"] = "radiation"; + args["linestyle"] = "--"; + args["color"] = "k"; + // args["marker"] = "x"; + matplotlibcpp::plot(xs, Trad, args); + + for (int g = 0; g < n_coll; ++g) { + std::map Trad_coll_args; + Trad_coll_args["label"] = fmt::format("group {}", g); + Trad_coll_args["linestyle"] = "-"; + Trad_coll_args["color"] = "C" + std::to_string(g); + matplotlibcpp::plot(xs, Trad_coll[g], Trad_coll_args); + } + + // Tgas_exact_args["label"] = "gas temperature (exact)"; + // Tgas_exact_args["color"] = "C0"; + // // Tgas_exact_args["marker"] = "x"; + // matplotlibcpp::plot(xs_exact, Tmat_exact, Tgas_exact_args); + + matplotlibcpp::xlim(0.0, 12.0); // cm + matplotlibcpp::ylim(0.0, 1000.0); // K + matplotlibcpp::xlabel("length x (cm)"); + matplotlibcpp::ylabel("temperature (K)"); + matplotlibcpp::legend(); + // matplotlibcpp::title(fmt::format("time t = {:.4g}", sim.tNew_[0])); + if (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + matplotlibcpp::title("PC"); + } else if (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + matplotlibcpp::title("PPL-fixed"); + } else if (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + matplotlibcpp::title("PPL-free"); + } + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./marshak_wave_Vaytet.pdf"); +#endif // HAVE_PYTHON + } + + // Cleanup and exit + std::cout << "Finished." << std::endl; + + // if ((rel_error > error_tol) || std::isnan(rel_error)) { + // status = 1; + // } + return status; +} diff --git a/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.hpp b/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.hpp new file mode 100644 index 000000000..575081b1d --- /dev/null +++ b/src/RadMarshakVaytet/test_radiation_marshak_Vaytet.hpp @@ -0,0 +1,26 @@ +#ifndef TEST_RADIATION_MARSHAK_VAYTET_HPP_ // NOLINT +#define TEST_RADIATION_MARSHAK_VAYTET_HPP_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radiation_marshak_vaytet.hpp +/// \brief Defines a test problem for radiation in the static diffusion regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// function definitions + +#endif // TEST_RADIATION_MARSHAK_VAYTET_HPP_ diff --git a/src/RadPulse/test_radiation_pulse.cpp b/src/RadPulse/test_radiation_pulse.cpp index b98db27f4..db75ee5c4 100644 --- a/src/RadPulse/test_radiation_pulse.cpp +++ b/src/RadPulse/test_radiation_pulse.cpp @@ -81,19 +81,6 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const return ComputePlanckOpacity(rho, Tgas); } -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, - const double Tgas) -> quokka::valarray -{ - quokka::valarray opacity_deriv{}; - double opacity_deriv_scalar = 0.; - if (Tgas > 1.0) { - opacity_deriv_scalar = (kappa0 / rho) * (3.0 / T0) * std::pow(Tgas / T0, 2); - } - opacity_deriv.fillin(opacity_deriv_scalar); - return opacity_deriv; -} - template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) { // extract variables required from the geom object diff --git a/src/RadTube/test_radiation_tube.cpp b/src/RadTube/test_radiation_tube.cpp index 229438ad5..9ab40230d 100644 --- a/src/RadTube/test_radiation_tube.cpp +++ b/src/RadTube/test_radiation_tube.cpp @@ -4,7 +4,7 @@ // Released under the MIT license. See LICENSE file included in the GitHub repo. //============================================================================== /// \file test_radiation_tube.cpp -/// \brief Defines a test problem for radiation pressure terms. +/// \brief Defines a test problem for radiation pressure terms. This is also a trivial test for the PPL_fixed_slope opacity model. /// #include @@ -36,6 +36,7 @@ constexpr double rho0 = 1.0; // g cm^-3 constexpr double T0 = 2.75e7; // K constexpr double rho1 = 2.1940476649492044; // g cm^-3 constexpr double T1 = 2.2609633884436745e7; // K +constexpr double a_rad = C::a_rad; constexpr double a0 = 4.0295519855200705e7; // cm s^-1 @@ -60,28 +61,45 @@ template <> struct Physics_Traits { template <> struct RadSystem_Traits { static constexpr double c_light = c_light_cgs_; static constexpr double c_hat = 10.0 * a0; - static constexpr double radiation_constant = radiation_constant_cgs_; + static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::k_B; - static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0., 3.3 * T0, inf}; // Kelvin + static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0.01 * T0, 3.3 * T0, 1000. * T0}; // Kelvin + // static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0.01 * T0, 1000. * T0}; // Kelvin static constexpr int beta_order = 1; - static constexpr OpacityModel opacity_model = OpacityModel::user; + // static constexpr OpacityModel opacity_model = OpacityModel::single_group; + static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; + // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; }; template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double /*rho*/, + const double /*Tgas*/) -> amrex::GpuArray, 2> { - quokka::valarray kappaPVec{}; - for (int g = 0; g < nGroups_; ++g) { - kappaPVec[g] = kappa0; + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_ + 1; ++i) { + exponents_and_values[0][i] = 0.0; + exponents_and_values[1][i] = kappa0; } - return kappaPVec; + return exponents_and_values; } -template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray -{ - return ComputePlanckOpacity(rho, Tgas); -} +// template <> +// AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray +// { +// quokka::valarray kappaPVec{}; +// for (int g = 0; g < nGroups_; ++g) { +// kappaPVec[g] = kappa0; +// } +// return kappaPVec; +// } + +// template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +// { +// return ComputePlanckOpacity(rho, Tgas); +// } // declare global variables // initial conditions read from file @@ -276,8 +294,8 @@ auto problem_main() -> int // Problem initialization RadhydroSimulation sim(BCs_cc); - sim.radiationReconstructionOrder_ = 2; // PLM - sim.reconstructionOrder_ = 2; // PLM + sim.radiationReconstructionOrder_ = 3; // PPM + sim.reconstructionOrder_ = 3; // PPM sim.stopTime_ = tmax; sim.cflNumber_ = CFL_number; sim.radiationCflNumber_ = CFL_number; @@ -323,27 +341,27 @@ auto problem_main() -> int } Erad_exact_arr[i] = Erad_0; Erad_arr[i] = Erad_t; - const double Trad_exact = std::pow(Erad_0 / radiation_constant_cgs_, 1. / 4.); - const double Trad = std::pow(Erad_t / radiation_constant_cgs_, 1. / 4.); + const double Trad_exact = std::pow(Erad_0 / a_rad, 1. / 4.); + const double Trad = std::pow(Erad_t / a_rad, 1. / 4.); Trad_arr[i] = Trad; Trad_exact_arr[i] = Trad_exact; Trad_err[i] = (Trad - Trad_exact) / Trad_exact; - double Egas_exact = values0.at(RadSystem::gasEnergy_index)[i]; - double x1GasMom_exact = values0.at(RadSystem::x1GasMomentum_index)[i]; - double x2GasMom_exact = values0.at(RadSystem::x2GasMomentum_index)[i]; - double x3GasMom_exact = values0.at(RadSystem::x3GasMomentum_index)[i]; + double const Egas_exact = values0.at(RadSystem::gasEnergy_index)[i]; + double const x1GasMom_exact = values0.at(RadSystem::x1GasMomentum_index)[i]; + double const x2GasMom_exact = values0.at(RadSystem::x2GasMomentum_index)[i]; + double const x3GasMom_exact = values0.at(RadSystem::x3GasMomentum_index)[i]; - double Egas = values.at(RadSystem::gasEnergy_index)[i]; - double x1GasMom = values.at(RadSystem::x1GasMomentum_index)[i]; - double x2GasMom = values.at(RadSystem::x2GasMomentum_index)[i]; - double x3GasMom = values.at(RadSystem::x3GasMomentum_index)[i]; + double const Egas = values.at(RadSystem::gasEnergy_index)[i]; + double const x1GasMom = values.at(RadSystem::x1GasMomentum_index)[i]; + double const x2GasMom = values.at(RadSystem::x2GasMomentum_index)[i]; + double const x3GasMom = values.at(RadSystem::x3GasMomentum_index)[i]; - double Eint_exact = RadSystem::ComputeEintFromEgas(rho_exact, x1GasMom_exact, x2GasMom_exact, x3GasMom_exact, Egas_exact); - double Tgas_exact = quokka::EOS::ComputeTgasFromEint(rho_exact, Eint_exact); + double const Eint_exact = RadSystem::ComputeEintFromEgas(rho_exact, x1GasMom_exact, x2GasMom_exact, x3GasMom_exact, Egas_exact); + double const Tgas_exact = quokka::EOS::ComputeTgasFromEint(rho_exact, Eint_exact); - double Eint = RadSystem::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Egas); - double Tgas = quokka::EOS::ComputeTgasFromEint(rho, Eint); + double const Eint = RadSystem::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Egas); + double const Tgas = quokka::EOS::ComputeTgasFromEint(rho, Eint); Tgas_arr[i] = Tgas; Tgas_err[i] = (Tgas - Tgas_exact) / Tgas_exact; diff --git a/src/RadhydroBB/CMakeLists.txt b/src/RadhydroBB/CMakeLists.txt new file mode 100644 index 000000000..55b16b1fa --- /dev/null +++ b/src/RadhydroBB/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_radhydro_bb test_radhydro_bb.cpp ../fextract.cpp ../interpolate.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_bb) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadhydroBB COMMAND test_radhydro_bb RadhydroBB.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/RadhydroBB/test_radhydro_bb.cpp b/src/RadhydroBB/test_radhydro_bb.cpp new file mode 100644 index 000000000..26ab83cc3 --- /dev/null +++ b/src/RadhydroBB/test_radhydro_bb.cpp @@ -0,0 +1,492 @@ +/// \file test_radhydro_bb.cpp +/// \brief Defines a test problem for blackbody spectrum in a uniform advecting medium. +/// + +#include +#include + +#include "AMReX_Array.H" +#include "AMReX_BC_TYPES.H" + +#include "AMReX_BLassert.H" +#include "ArrayUtil.hpp" +#include "fextract.hpp" +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// #include "AMReX_BC_TYPES.H" +#include "AMReX_IntVect.H" +#include "AMReX_Print.H" +// #include "RadhydroSimulation.hpp" +// #include "fextract.hpp" +#include "physics_info.hpp" +// #include "radiation_system.hpp" +#include "test_radhydro_bb.hpp" + +static constexpr bool export_csv = true; + +struct PulseProblem { +}; // dummy type to allow compile-type polymorphism via template specialization + +constexpr int n_groups_ = 4; +// constexpr int n_groups_ = 8; +// constexpr int n_groups_ = 16; +// constexpr int n_groups_ = 64; + +constexpr amrex::GpuArray rad_boundaries_ = []() constexpr { + if constexpr (n_groups_ == 1) { + return amrex::GpuArray{0.0, inf}; + } else if constexpr (n_groups_ == 4) { + return amrex::GpuArray{1.00000000e-03, 1.77827941e-02, 3.16227766e-01, 5.62341325e+00, 1.00000000e+02}; + } else if constexpr (n_groups_ == 8) { + return amrex::GpuArray{1.00000000e-03, 4.21696503e-03, 1.77827941e-02, 7.49894209e-02, 3.16227766e-01, + 1.33352143e+00, 5.62341325e+00, 2.37137371e+01, 1.00000000e+02}; + } else if constexpr (n_groups_ == 16) { + return amrex::GpuArray{1.00000000e-03, 2.05352503e-03, 4.21696503e-03, 8.65964323e-03, 1.77827941e-02, 3.65174127e-02, + 7.49894209e-02, 1.53992653e-01, 3.16227766e-01, 6.49381632e-01, 1.33352143e+00, 2.73841963e+00, + 5.62341325e+00, 1.15478198e+01, 2.37137371e+01, 4.86967525e+01, 1.00000000e+02}; + } else if constexpr (n_groups_ == 64) { + return amrex::GpuArray{ + 1.00000000e-03, 1.19708503e-03, 1.43301257e-03, 1.71543790e-03, 2.05352503e-03, 2.45824407e-03, 2.94272718e-03, 3.52269465e-03, + 4.21696503e-03, 5.04806572e-03, 6.04296390e-03, 7.23394163e-03, 8.65964323e-03, 1.03663293e-02, 1.24093776e-02, 1.48550802e-02, + 1.77827941e-02, 2.12875166e-02, 2.54829675e-02, 3.05052789e-02, 3.65174127e-02, 4.37144481e-02, 5.23299115e-02, 6.26433537e-02, + 7.49894209e-02, 8.97687132e-02, 1.07460783e-01, 1.28639694e-01, 1.53992653e-01, 1.84342299e-01, 2.20673407e-01, 2.64164832e-01, + 3.16227766e-01, 3.78551525e-01, 4.53158364e-01, 5.42469094e-01, 6.49381632e-01, 7.77365030e-01, 9.30572041e-01, 1.11397386e+00, + 1.33352143e+00, 1.59633854e+00, 1.91095297e+00, 2.28757320e+00, 2.73841963e+00, 3.27812115e+00, 3.92418976e+00, 4.69758882e+00, + 5.62341325e+00, 6.73170382e+00, 8.05842188e+00, 9.64661620e+00, 1.15478198e+01, 1.38237223e+01, 1.65481710e+01, 1.98095678e+01, + 2.37137371e+01, 2.83873596e+01, 3.39820833e+01, 4.06794432e+01, 4.86967525e+01, 5.82941535e+01, 6.97830585e+01, 8.35362547e+01, + 1.00000000e+02}; + } +}(); + +constexpr double c = 1.0e8; +// model 0 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-4 * c; +// constexpr double kappa0 = 1.0e4; // dx = 1, tau = kappa0 * dx = 1e4 +// constexpr double chat = 1.0e7; +// model 1 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-4 * c; +// constexpr double kappa0 = 1.0e4; // dx = 1, tau = kappa0 * dx = 1e4 +// constexpr double chat = 1.0e8; +// model 2 +// constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 1e-2 * c; +// constexpr double kappa0 = 1.0e5; +// constexpr double chat = 1.0e8; +// model 3 +constexpr int beta_order_ = 1; // order of beta in the radiation four-force +// constexpr double v0 = 0.0; +// constexpr double v0 = 1e-2 * c; +// constexpr double v0 = 0.3 * c; +constexpr double v0 = 0.001 * c; +constexpr double kappa0 = 1.0e5; +constexpr double chat = c; + +constexpr double T0 = 1.0; // temperature +constexpr double rho0 = 1.0; // matter density +constexpr double a_rad = 1.0; +constexpr double mu = 1.0; +constexpr double k_B = 1.0; + +constexpr double nu_unit = 1.0; +constexpr double T_equilibrium = 0.768032502191; + +// static diffusion, beta = 1e-4, tau_cell = kappa0 * dx = 100, beta tau_cell = 1e-2 +// constexpr double kappa0 = 100.; // cm^2 g^-1 +// constexpr double v0 = 1.0e-4 * c; // advecting pulse + +// dynamic diffusion, beta = 1e-3, tau = kappa0 * dx = 1e5, beta tau = 100 +// constexpr double max_time = 10.0 / v0; +// constexpr double max_time = 1000.0 / (c * rho0 * kappa0); // dt >> 1 / (c * chi) +constexpr double max_time = 10.0 / (1e-2 * c); + +constexpr double Erad0 = a_rad * T0 * T0 * T0 * T0; +constexpr double Erad_beta2 = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad0; +constexpr double erad_floor = a_rad * 1e-30; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_groups_; +}; + +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr int beta_order = beta_order_; + static constexpr double energy_unit = nu_unit; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; + static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; +}; + +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double /*rho*/, + const double /*Tgas*/) -> amrex::GpuArray, 2> +{ + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_ + 1; ++i) { + exponents_and_values[0][i] = 0.0; + } + for (int i = 0; i < nGroups_ + 1; ++i) { + exponents_and_values[1][i] = kappa0; + } + return exponents_and_values; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray +{ + quokka::valarray kappaPVec{}; + for (int i = 0; i < nGroups_; ++i) { + kappaPVec[i] = kappa0; + } + return kappaPVec; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +{ + return ComputePlanckOpacity(rho, Tgas); +} + +AMREX_GPU_HOST_DEVICE +auto compute_exact_bb(const double nu, const double T) -> double +{ + double const x = nu_unit * nu / (k_B * T); + double const coeff = nu_unit / (k_B * T); + double const planck_integral = std::pow(x, 3) / (std::exp(x) - 1.0); + return coeff * planck_integral / (std::pow(PI, 4) / 15.0) * (a_rad * std::pow(T, 4)); +} + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +{ + // extract variables required from the geom object + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T0); + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int g = 0; g < n_groups_; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = erad_floor; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0.; + } + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho0 * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho0; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +auto problem_main() -> int +{ + // This problem is a *linear* radiation diffusion problem, i.e. + // parameters are chosen such that the radiation and gas temperatures + // should be near equilibrium, and the opacity is chosen to go as + // T^3, such that the radiation diffusion equation becomes linear in T. + + // This makes this problem a stringent test of the asymptotic- + // preserving property of the computational method, since the + // optical depth per cell at the peak of the temperature profile is + // of order 10^5. + + // Problem parameters + const int max_timesteps = 1e6; + const double CFL_number_gas = 0.8; + const double CFL_number_rad = 8.0; + + const double max_dt = 1.0; + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); + } + } + + // Problem initialization + RadhydroSimulation sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + sim.radiationCflNumber_ = CFL_number_rad; + sim.cflNumber_ = CFL_number_gas; + sim.maxDt_ = max_dt; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read output variables + auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); + const int nx = static_cast(position.size()); + + std::vector xs(nx); + std::vector Trad(nx); + std::vector Tgas(nx); + std::vector Erad(nx); + std::vector Egas(nx); + std::vector Vgas(nx); + std::vector Vgas_exact(nx); + std::vector rhogas(nx); + std::vector Trad_exact{}; + std::vector Tgas_exact{}; + std::vector Erad_exact{}; + + for (int i = 0; i < nx; ++i) { + amrex::Real const x = position[i]; + xs.at(i) = x; + // const auto Erad_t = values.at(RadSystem::radEnergy_index)[i]; + double Erad_t = 0.0; + for (int g = 0; g < n_groups_; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + } + const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); + const auto rho_t = values.at(RadSystem::gasDensity_index)[i]; + const auto v_t = values.at(RadSystem::x1GasMomentum_index)[i] / rho_t; + rhogas.at(i) = rho_t; + Erad.at(i) = Erad_t; + Trad.at(i) = Trad_t / T0; + Egas.at(i) = values.at(RadSystem::gasInternalEnergy_index)[i]; + Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas.at(i)) / T0; + Tgas_exact.push_back(1.0); + Vgas.at(i) = v_t; + Vgas_exact.at(i) = v0; + + auto const Erad_val = a_rad * std::pow(T0, 4); + Trad_exact.push_back(T_equilibrium); + Erad_exact.push_back(Erad_val); + } + + // compute spectrum + std::vector E_r{}; // radiation energy density + std::vector F_r{}; // flux density + std::vector E_r_spec{}; // specific radiation energy density + std::vector F_r_spec{}; // specific flux density + std::vector bin_center{}; + int const ii = 3; // a random grid + for (int g = 0; g < n_groups_; ++g) { + // bin_center.push_back(std::sqrt(rad_boundaries_[g] * rad_boundaries_[g + 1])); + bin_center.push_back(rad_boundaries_[g]); + const auto Erad_t = values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[ii]; + const auto Frad_t = values.at(RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g)[ii]; + const double bin_width = rad_boundaries_[g + 1] - rad_boundaries_[g]; + E_r.push_back(Erad_t); + F_r.push_back(Frad_t); + F_r_spec.push_back(Frad_t / bin_width); + E_r_spec.push_back(Erad_t / bin_width); + } + + // Read in exact solution + std::vector nu_exact; + std::vector F_read; + std::vector E_read; + + std::string const filename = "../extern/Doppler-spectrum/exact_flux_density.csv"; + int const nbins = 1024; + int const bins_per_group = nbins / n_groups_; + std::ifstream fstream(filename, std::ios::in); + AMREX_ASSERT(fstream.is_open()); + + double err_norm = 0.; + double sol_norm = 0.; + + for (std::string line; std::getline(fstream, line);) { + std::istringstream iss(line); + std::vector values_line; + std::string value; + + while (std::getline(iss, value, ',')) { + values_line.push_back(std::stod(value)); + } + double const nu_val = values_line.at(0); // dimensionless + double const Fnu_val = values_line.at(1); // dimensionless + nu_exact.push_back(nu_val); + F_read.push_back(Fnu_val); + double const enu = compute_exact_bb(nu_val, T_equilibrium); + E_read.push_back(enu); + } + // insert a dummy breakpoint + int aa = 0; + aa += 1; + std::cout << aa << std::endl; + + // assert nu_exact[0] = 0.001 and nu_exact[-1] = 100 + AMREX_ASSERT(nu_exact[0] == 0.001); + AMREX_ASSERT(nu_exact.back() == 100.0); + AMREX_ASSERT(nu_exact.size() == nbins + 1); + // compute group-integrated exact solution + std::vector Fnu_exact(n_groups_); + std::vector Enu_exact(n_groups_); + int count = 0; + double sum_F = 0.0; + double sum_E = 0.0; + for (int i = 0; i < nbins; ++i) { + sum_F += F_read[i] * (nu_exact[i + 1] - nu_exact[i]); + sum_E += E_read[i] * (nu_exact[i + 1] - nu_exact[i]); + if ((i + 1) % bins_per_group == 0) { + Fnu_exact[count] = sum_F / (rad_boundaries_[count + 1] - rad_boundaries_[count]); + Enu_exact[count] = sum_E / (rad_boundaries_[count + 1] - rad_boundaries_[count]); + sum_F = 0.0; + sum_E = 0.0; + count++; + } + } + + // compute error norm + + // compare temperature with T_equilibrium + double err_norm_T = 0.; + double sol_norm_T = 0.; + for (int i = 0; i < nx; ++i) { + err_norm_T += std::abs(Trad[i] - Trad_exact[i]); + sol_norm_T += std::abs(Trad_exact[i]); + } + const double rel_error_T = err_norm_T / sol_norm_T; + amrex::Print() << "Relative L1 error norm for T_gas = " << rel_error_T << std::endl; + + for (int g = 0; g < n_groups_; ++g) { + err_norm += std::abs(Fnu_exact[g] - F_r_spec[g]); + sol_norm += std::abs(Fnu_exact[g]); + } + + const double error_tol = 0.1; + const double rel_error = err_norm / sol_norm; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + +#ifdef HAVE_PYTHON + // plot temperature + matplotlibcpp::clf(); + std::map Trad_args; + std::map Tgas_args; + std::map Tradexact_args; + Trad_args["label"] = "T_rad (numerical)"; + Trad_args["linestyle"] = "-"; + Tradexact_args["label"] = "T (exact)"; + Tradexact_args["linestyle"] = "--"; + Tgas_args["label"] = "T_gas (numerical)"; + Tgas_args["linestyle"] = "-"; + matplotlibcpp::plot(xs, Trad, Trad_args); + matplotlibcpp::plot(xs, Trad_exact, Tradexact_args); + matplotlibcpp::plot(xs, Tgas, Tgas_args); + matplotlibcpp::xlabel("x (dimensionless)"); + matplotlibcpp::ylabel("temperature (dimensionless)"); + matplotlibcpp::legend(); + matplotlibcpp::ylim(0.0, 1.0); + matplotlibcpp::title(fmt::format("time ct = {:.4g}", sim.tNew_[0] * c)); + // if constexpr (beta_order_ == 1) { + // matplotlibcpp::ylim(1.0 - 1.0e-7, 1.0 + 1.0e-7); + // } + matplotlibcpp::tight_layout(); + // matplotlibcpp::save("./adv_temp.pdf"); + // save to adv_temp_{n_groups_}bins.pdf + matplotlibcpp::save(fmt::format("./adv_temp_{}_bins.pdf", n_groups_)); + + // plot spectrum + matplotlibcpp::clf(); + std::unordered_map spec_args; + spec_args["label"] = "spectrum"; + spec_args["color"] = "C0"; + matplotlibcpp::scatter(bin_center, E_r_spec, 10.0, spec_args); + std::map spec_exact_args; + spec_exact_args["label"] = "spectrum (exact)"; + spec_exact_args["linestyle"] = "-"; + spec_exact_args["color"] = "C1"; + matplotlibcpp::plot(bin_center, Enu_exact, spec_exact_args); + // log-log + matplotlibcpp::xscale("log"); + matplotlibcpp::yscale("log"); + matplotlibcpp::ylim(1.0e-8, 1.0e0); + matplotlibcpp::xlabel("frequency (dimensionless)"); + matplotlibcpp::ylabel("spectrum density (dimensionless)"); + // matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time ct = {:.4g}", sim.tNew_[0] * c)); + matplotlibcpp::tight_layout(); + // matplotlibcpp::save("./adv_spectrum.pdf"); + // save to adv_spectrum_{n_groups_}bins.pdf + matplotlibcpp::save(fmt::format("./adv_spectrum_{}_bins.pdf", n_groups_)); + + // plot flux spectrum + matplotlibcpp::clf(); + std::unordered_map Frad_args; + Frad_args["label"] = "flux (exact)"; + Frad_args["color"] = "C0"; + matplotlibcpp::scatter(bin_center, Fnu_exact, 10.0, Frad_args); + std::map Frad_exact_args; + Frad_exact_args["label"] = "flux (numerical)"; + Frad_exact_args["linestyle"] = "-"; + Frad_exact_args["color"] = "C1"; + matplotlibcpp::plot(bin_center, F_r_spec, Frad_exact_args); + // log-log + matplotlibcpp::xscale("log"); + matplotlibcpp::yscale("log"); + matplotlibcpp::ylim(1.0e-3, 1.0e5); + matplotlibcpp::xlabel("frequency (dimensionless)"); + matplotlibcpp::ylabel("flux density (dimensionless)"); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("time ct = {:.4g}", sim.tNew_[0] * c)); + matplotlibcpp::tight_layout(); + // matplotlibcpp::save("./adv_flux_spectrum.pdf"); + // save to adv_flux_spectrum_{n_groups_}bins.pdf + matplotlibcpp::save(fmt::format("./adv_flux_spectrum_{}_bins.pdf", n_groups_)); +#endif + + if (export_csv) { + std::ofstream file; + file.open("adv_spectrum.csv"); + file << "nu_Left, E_r (erg/cm^3/Hz)\n"; + for (int g = 0; g < n_groups_; ++g) { + file << std::scientific << std::setprecision(12) << rad_boundaries_[g] << "," << E_r_spec[g] << "\n"; + } + file.close(); + + file.open("adv_flux_spectrum.csv"); + file << "nu_Left, F_r (erg/cm^2/s/Hz)\n"; + for (int g = 0; g < n_groups_; ++g) { + file << std::scientific << std::setprecision(12) << rad_boundaries_[g] << "," << F_r_spec[g] << "\n"; + } + file.close(); + + file.open("adv_temp.csv"); + file << "xs,Tgas,Trad\n"; + for (size_t i = 0; i < xs.size(); ++i) { + file << std::scientific << std::setprecision(12) << xs[i] << "," << Tgas[i] << "," << Trad[i] << "\n"; + } + file.close(); + } + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error) || (rel_error_T > error_tol) || std::isnan(rel_error_T)) { + status = 1; + } + return status; +} diff --git a/src/RadhydroBB/test_radhydro_bb.hpp b/src/RadhydroBB/test_radhydro_bb.hpp new file mode 100644 index 000000000..3b859dc7e --- /dev/null +++ b/src/RadhydroBB/test_radhydro_bb.hpp @@ -0,0 +1,29 @@ +#ifndef TEST_RADHYDRO_BB_HPP_ // NOLINT +#define TEST_RADHYDRO_BB_HPP_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_radhydro_bb.hpp +/// \brief Defines a test problem for blackbody spectrum in a uniform advecting medium. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "RadhydroSimulation.hpp" +#include "hydro_system.hpp" +#include "interpolate.hpp" +#include "radiation_system.hpp" + +// function definitions +auto problem_main() -> int; + +#endif // TEST_RADHYDRO_BB_HPP_ diff --git a/src/RadhydroPulseMG/CMakeLists.txt b/src/RadhydroPulseMG/CMakeLists.txt deleted file mode 100644 index 3b9994821..000000000 --- a/src/RadhydroPulseMG/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -add_executable(test_radhydro_pulse_MG test_radhydro_pulse_MG.cpp ../fextract.cpp ${QuokkaObjSources}) - -if(AMReX_GPU_BACKEND MATCHES "CUDA") - setup_target_for_cuda_compilation(test_radhydro_pulse_MG) -endif(AMReX_GPU_BACKEND MATCHES "CUDA") - -add_test(NAME RadhydroPulseMG COMMAND test_radhydro_pulse_MG RadhydroPulse.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/RadhydroPulseMGconst/CMakeLists.txt b/src/RadhydroPulseMGconst/CMakeLists.txt new file mode 100644 index 000000000..2f716509a --- /dev/null +++ b/src/RadhydroPulseMGconst/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_radhydro_pulse_MG_const_kappa test_radhydro_pulse_MG_const_kappa.cpp ../fextract.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_radhydro_pulse_MG_const_kappa) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadhydroPulseMGconst COMMAND test_radhydro_pulse_MG_const_kappa RadhydroPulse.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp b/src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp similarity index 55% rename from src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp rename to src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp index 415d2ee3e..d64fecc0b 100644 --- a/src/RadhydroPulseMG/test_radhydro_pulse_MG.cpp +++ b/src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp @@ -1,17 +1,20 @@ -/// \file test_radhydro_pulse_MG.cpp -/// \brief Defines a test problem for multigroup radiation in the diffusion regime with advection by gas. +/// \file test_radhydro_pulse_MG_const_kappa.cpp +/// \brief Defines a test problem for multigroup radiation in the diffusion regime with advection by gas, running +/// with PPL_opacity_fixed_slope_spectrum opacity model. /// -#include "test_radhydro_pulse_MG.hpp" +#include "test_radhydro_pulse_MG_const_kappa.hpp" #include "AMReX_BC_TYPES.H" #include "AMReX_Print.H" #include "RadhydroSimulation.hpp" #include "fextract.hpp" #include "physics_info.hpp" -struct PulseProblem { +// Single-group problem +struct SGProblem { }; // dummy type to allow compile-type polymorphism via template specialization -struct AdvPulseProblem { +// Multi-group problem +struct MGproblem { }; constexpr int isconst = 0; @@ -20,9 +23,7 @@ constexpr int isconst = 0; // constexpr int n_groups_ = 2; // constexpr amrex::GpuArray rad_boundaries_{1e15, 1e17, 1e19}; constexpr int n_groups_ = 4; -constexpr amrex::GpuArray rad_boundaries_{1e16, 1e17, 1e18, 1e19, 1e20}; -// constexpr int n_groups_ = 8; -// constexpr amrex::GpuArray rad_boundaries_{1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19, 3.16e19, 1e20}; +constexpr amrex::GpuArray rad_boundaries_{1e15, 1e16, 1e17, 1e18, 1e19}; // constexpr int n_groups_ = 64; // constexpr amrex::GpuArray rad_boundaries_{1.00000000e+15, 1.15478198e+15, 1.33352143e+15, 1.53992653e+15, // 1.77827941e+15, 2.05352503e+15, 2.37137371e+15, 2.73841963e+15, @@ -42,57 +43,58 @@ constexpr amrex::GpuArray rad_boundaries_{1e16, 1e17, 1e1 // 5.62341325e+18, 6.49381632e+18, 7.49894209e+18, 8.65964323e+18, // 1.00000000e+19}; -constexpr double kappa0 = 180.; // cm^2 g^-1 -constexpr double kappa_const = 100.; // cm^2 g^-1 -constexpr double T0 = 1.0e7; // K (temperature) -constexpr double T1 = 2.0e7; // K (temperature) -constexpr double rho0 = 1.2; // g cm^-3 (matter density) +constexpr double kappa0 = 100.; // cm^2 g^-1 +constexpr double T0 = 1.0e7; // K (temperature) +constexpr double T1 = 2.0e7; // K (temperature) +constexpr double rho0 = 1.2; // g cm^-3 (matter density) constexpr double a_rad = C::a_rad; constexpr double c = C::c_light; // speed of light (cgs) constexpr double chat = c; constexpr double width = 24.0; // cm, width of the pulse -constexpr double erad_floor = a_rad * T0 * T0 * T0 * T0 * 1.0e-14; +constexpr double Erad0 = a_rad * T0 * T0 * T0 * T0; +constexpr double erad_floor = Erad0 * 1.0e-14; constexpr double mu = 2.33 * C::m_u; constexpr double h_planck = C::hplanck; constexpr double k_B = C::k_B; -constexpr double v0_nonadv = 0.; // non-advecting pulse -// static diffusion: (for single group) tau = 2e3, beta = 3e-5, beta tau = 6e-2 -constexpr double v0_adv = 1.0e6; // advecting pulse -constexpr double max_time = 4.8e-5; // max_time = 0.02 * width / v1; -constexpr int64_t max_timesteps = 30; +// testing parameters: +// dynamic diffusion limit: tau = 2e3, beta = 7e-3, beta tau = 14 +constexpr double v0 = 2.0e8; // advecting speed +constexpr double max_time = 4.8e-5; // diffusion distance is about a few pulse width +constexpr int64_t max_timesteps = 1e2; // for fast testing // dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 // constexpr double kappa0 = 1000.; // cm^2 g^-1 -// constexpr double v0_adv = 1.0e8; // advecting pulse +// constexpr double v0 = 1.0e8; // advecting pulse // constexpr double max_time = 1.2e-4; // max_time = 2.0 * width / v1; constexpr double T_ref = T0; constexpr double nu_ref = 1.0e18; // Hz constexpr double coeff_ = h_planck * nu_ref / (k_B * T_ref); -template <> struct quokka::EOS_Traits { - static constexpr double mean_molecular_weight = mu; - static constexpr double boltzmann_constant = k_B; - static constexpr double gamma = 5. / 3.; -}; -template <> struct quokka::EOS_Traits { +AMREX_GPU_HOST_DEVICE +auto compute_initial_Tgas(const double x) -> double +{ + // compute temperature profile for Gaussian radiation pulse + const double sigma = width; + return T0 + (T1 - T0) * std::exp(-x * x / (2.0 * sigma * sigma)); +} + +AMREX_GPU_HOST_DEVICE +auto compute_exact_rho(const double x) -> double +{ + // compute density profile for Gaussian radiation pulse + auto T = compute_initial_Tgas(x); + return rho0 * T0 / T + (a_rad * mu / 3. / k_B) * (std::pow(T0, 4) / T - std::pow(T, 3)); +} + +template <> struct quokka::EOS_Traits { static constexpr double mean_molecular_weight = mu; static constexpr double boltzmann_constant = k_B; static constexpr double gamma = 5. / 3.; }; -template <> struct Physics_Traits { - // cell-centred - static constexpr bool is_hydro_enabled = true; - static constexpr int numMassScalars = 0; // number of mass scalars - static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars - static constexpr bool is_radiation_enabled = true; - // face-centred - static constexpr bool is_mhd_enabled = false; - static constexpr int nGroups = n_groups_; -}; -template <> struct Physics_Traits { +template <> struct Physics_Traits { // cell-centred static constexpr bool is_hydro_enabled = true; static constexpr int numMassScalars = 0; // number of mass scalars @@ -100,132 +102,33 @@ template <> struct Physics_Traits { static constexpr bool is_radiation_enabled = true; // face-centred static constexpr bool is_mhd_enabled = false; - static constexpr int nGroups = n_groups_; + static constexpr int nGroups = 1; }; -template <> struct RadSystem_Traits { - static constexpr double c_light = c; - static constexpr double c_hat = chat; - static constexpr double radiation_constant = a_rad; - static constexpr double Erad_floor = erad_floor; - static constexpr double energy_unit = h_planck; - static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; - static constexpr int beta_order = 1; -}; -template <> struct RadSystem_Traits { +template <> struct RadSystem_Traits { static constexpr double c_light = c; static constexpr double c_hat = chat; static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; - static constexpr double energy_unit = h_planck; - static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr int beta_order = 1; }; -AMREX_GPU_HOST_DEVICE -auto compute_initial_Tgas(const double x) -> double -{ - // compute temperature profile for Gaussian radiation pulse - const double sigma = width; - return T0 + (T1 - T0) * std::exp(-x * x / (2.0 * sigma * sigma)); -} - -AMREX_GPU_HOST_DEVICE -auto compute_exact_rho(const double x) -> double -{ - // compute density profile for Gaussian radiation pulse - auto T = compute_initial_Tgas(x); - return rho0 * T0 / T + (a_rad * mu / 3. / k_B) * (std::pow(T0, 4) / T - std::pow(T, 3)); -} - -AMREX_GPU_HOST_DEVICE -auto compute_kappa(const double nu, const double Tgas) -> double -{ - // cm^-1 - auto T_ = Tgas / T_ref; - auto nu_ = nu / nu_ref; - return kappa0 * std::pow(T_, -0.5) * std::pow(nu_, -3) * (1.0 - std::exp(-coeff_ * nu_ / T_)); -} - -AMREX_GPU_HOST_DEVICE -auto compute_repres_nu(amrex::GpuArray rad_boundaries) -> quokka::valarray -{ - // return the geometrical mean as the representative frequency for each group - quokka::valarray nu_rep{}; - if (n_groups_ == 1) { - nu_rep[0] = nu_ref; - } else { - for (int g = 0; g < n_groups_; ++g) { - nu_rep[g] = std::sqrt(rad_boundaries[g] * rad_boundaries[g + 1]); - } - } - return nu_rep; -} - -template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray { quokka::valarray kappaPVec{}; - auto nu_rep = compute_repres_nu(rad_boundaries_); for (int g = 0; g < nGroups_; ++g) { - kappaPVec[g] = compute_kappa(nu_rep[g], Tgas) / rho; - } - if constexpr (isconst == 1) { - kappaPVec.fillin(kappa_const); + kappaPVec[g] = kappa0; } return kappaPVec; } -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray -{ - return RadSystem::ComputePlanckOpacity(rho, Tgas); -} -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray -{ - return ComputePlanckOpacity(rho, Tgas); -} -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray { return ComputePlanckOpacity(rho, Tgas); } -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, - const double Tgas) -> quokka::valarray -{ - quokka::valarray opacity_deriv{}; - const auto nu_rep = compute_repres_nu(rad_boundaries_); - const auto T = Tgas / T0; - for (int g = 0; g < nGroups_; ++g) { - const auto nu = nu_rep[g] / nu_ref; - opacity_deriv[g] = - 1. / T0 * kappa0 * (-0.5 * std::pow(T, -1.5) * (1. - std::exp(-coeff_ * nu / T)) - coeff_ * std::pow(T, -2.5) * std::exp(-coeff_ * nu / T)); - opacity_deriv[g] /= rho; - } - if constexpr (isconst == 1) { - opacity_deriv.fillin(0.0); - } - return opacity_deriv; -} -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double rho, - const double Tgas) -> quokka::valarray -{ - return RadSystem::ComputePlanckOpacityTempDerivative(rho, Tgas); -} - -template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double -{ - return (1. / 3.); // Eddington approximation -} -template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputeEddingtonFactor(double /*f*/) -> double -{ - return (1. / 3.); // Eddington approximation -} - -template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) { // extract variables required from the geom object amrex::GpuArray const dx = grid_elem.dx_; @@ -236,35 +139,74 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(qu amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); - const auto radBoundaries_g = RadSystem_Traits::radBoundaries; - // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; const double Trad = compute_initial_Tgas(x - x0); const double rho = compute_exact_rho(x - x0); - const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); - const double v0 = v0_nonadv; + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); + const double Erad = a_rad * Trad * Trad * Trad * Trad; + + state_cc(i, j, k, RadSystem::radEnergy_index) = Erad; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::x2RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index) = 0; + + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} - auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; - state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; - state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; - state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; - state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; - } +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars + static constexpr bool is_radiation_enabled = true; + // face-centred + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = n_groups_; +}; - state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; - state_cc(i, j, k, RadSystem::gasDensity_index) = rho; - state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; - state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; - state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; - state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; - }); +template <> struct RadSystem_Traits { + static constexpr double c_light = c; + static constexpr double c_hat = chat; + static constexpr double radiation_constant = a_rad; + static constexpr double Erad_floor = erad_floor; + static constexpr double energy_unit = h_planck; + static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; + static constexpr int beta_order = 1; + // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; + static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; + // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; +}; + +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray const /* rad_boundaries */, const double /* rho */, + const double /* Tgas */) -> amrex::GpuArray, 2> +{ + amrex::GpuArray exponents{}; + amrex::GpuArray kappa_lower{}; + for (int g = 0; g < nGroups_ + 1; ++g) { + exponents[g] = 0.; + kappa_lower[g] = kappa0; + } + amrex::GpuArray, 2> const exponents_and_values{exponents, kappa_lower}; + return exponents_and_values; } -template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) + +template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) { // extract variables required from the geom object amrex::GpuArray const dx = grid_elem.dx_; @@ -275,32 +217,30 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid amrex::Real const x0 = prob_lo[0] + 0.5 * (prob_hi[0] - prob_lo[0]); - const auto radBoundaries_g = RadSystem_Traits::radBoundaries; + const auto radBoundaries_g = RadSystem_Traits::radBoundaries; // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { amrex::Real const x = prob_lo[0] + (i + static_cast(0.5)) * dx[0]; const double Trad = compute_initial_Tgas(x - x0); const double rho = compute_exact_rho(x - x0); - const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); - const double v0 = v0_adv; + const double Egas = quokka::EOS::ComputeEintFromTgas(rho, Trad); - auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - // state_cc(i, j, k, RadSystem::radEnergy_index) = (1. + 4. / 3. * (v0 * v0) / (c * c)) * Erad; - state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; - state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; - state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; - state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; } - state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; - state_cc(i, j, k, RadSystem::gasDensity_index) = rho; - state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; - state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; - state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; - state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; }); } @@ -316,8 +256,10 @@ auto problem_main() -> int const double max_dt = 1e-3; // t_cr = 2 cm / cs = 7e-8 s + // Problem 1: pulse with grey radiation + // Boundary conditions - constexpr int nvars = RadSystem::nvar_; + constexpr int nvars = RadSystem::nvar_; amrex::Vector BCs_cc(nvars); for (int n = 0; n < nvars; ++n) { // periodic boundary condition in the x-direction will not work @@ -329,10 +271,8 @@ auto problem_main() -> int } } - // Problem 1: non-advecting pulse - // Problem initialization - RadhydroSimulation sim(BCs_cc); + RadhydroSimulation sim(BCs_cc); sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; @@ -351,8 +291,6 @@ auto problem_main() -> int // read output variables auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); const int nx = static_cast(position.size()); - amrex::GpuArray prob_lo = sim.geom[0].ProbLoArray(); - amrex::GpuArray prob_hi = sim.geom[0].ProbHiArray(); std::vector xs(nx); std::vector Trad(nx); @@ -363,30 +301,44 @@ auto problem_main() -> int for (int i = 0; i < nx; ++i) { amrex::Real const x = position[i]; xs.at(i) = x; - // const auto Erad_t = values.at(RadSystem::radEnergy_index)[i]; + // const auto Erad_t = values.at(RadSystem::radEnergy_index)[i]; double Erad_t = 0.0; - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; } const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); - const auto rho_t = values.at(RadSystem::gasDensity_index)[i]; - const auto v_t = values.at(RadSystem::x1GasMomentum_index)[i] / rho_t; - const auto Egas = values.at(RadSystem::gasInternalEnergy_index)[i]; + const auto rho_t = values.at(RadSystem::gasDensity_index)[i]; + const auto v_t = values.at(RadSystem::x1GasMomentum_index)[i] / rho_t; + const auto Egas = values.at(RadSystem::gasInternalEnergy_index)[i]; rhogas.at(i) = rho_t; Trad.at(i) = Trad_t; - Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); + Tgas.at(i) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); Vgas.at(i) = 1e-5 * v_t; } // END OF PROBLEM 1 // Problem 2: advecting pulse + // Boundary conditions + constexpr int nvars2 = RadSystem::nvar_; + amrex::Vector BCs_cc2(nvars2); + for (int n = 0; n < nvars2; ++n) { + // periodic boundary condition in the x-direction will not work + BCs_cc2[n].setLo(0, amrex::BCType::foextrap); // extrapolate + BCs_cc2[n].setHi(0, amrex::BCType::foextrap); + for (int i = 1; i < AMREX_SPACEDIM; ++i) { + BCs_cc2[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc2[n].setHi(i, amrex::BCType::int_dir); + } + } + // Problem initialization - RadhydroSimulation sim2(BCs_cc); + RadhydroSimulation sim2(BCs_cc2); sim2.radiationReconstructionOrder_ = 3; // PPM sim2.stopTime_ = max_time; sim2.radiationCflNumber_ = CFL_number; + sim2.cflNumber_ = CFL_number; sim2.maxDt_ = max_dt; sim2.maxTimesteps_ = max_timesteps; sim2.plotfileInterval_ = -1; @@ -399,52 +351,54 @@ auto problem_main() -> int // read output variables auto [position2, values2] = fextract(sim2.state_new_cc_[0], sim2.Geom(0), 0, 0.0); - prob_lo = sim2.geom[0].ProbLoArray(); - prob_hi = sim2.geom[0].ProbHiArray(); + const int nx2 = static_cast(position2.size()); + + amrex::GpuArray prob_lo = sim2.geom[0].ProbLoArray(); + amrex::GpuArray prob_hi = sim2.geom[0].ProbHiArray(); // compute the pixel size - const double dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx); - const double move = v0_adv * sim2.tNew_[0]; + const double dx = (prob_hi[0] - prob_lo[0]) / static_cast(nx2); + const double move = v0 * sim2.tNew_[0]; const int n_p = static_cast(move / dx); - const int half = static_cast(nx / 2.0); + const int half = static_cast(nx2 / 2.0); const double drift = move - static_cast(n_p) * dx; - const int shift = n_p - static_cast((n_p + half) / nx) * nx; + const int shift = n_p - static_cast((n_p + half) / nx2) * nx2; - std::vector xs2(nx); - std::vector Trad2(nx); - std::vector Tgas2(nx); - std::vector Vgas2(nx); - std::vector rhogas2(nx); + std::vector xs2(nx2); + std::vector Trad2(nx2); + std::vector Tgas2(nx2); + std::vector Vgas2(nx2); + std::vector rhogas2(nx2); - for (int i = 0; i < nx; ++i) { + for (int i = 0; i < nx2; ++i) { int index_ = 0; if (shift >= 0) { if (i < shift) { - index_ = nx - shift + i; + index_ = nx2 - shift + i; } else { index_ = i - shift; } } else { - if (i <= nx - 1 + shift) { + if (i <= nx2 - 1 + shift) { index_ = i - shift; } else { - index_ = i - (nx + shift); + index_ = i - (nx2 + shift); } } const amrex::Real x = position2[i]; - // const auto Erad_t = values2.at(RadSystem::radEnergy_index)[i]; + // const auto Erad_t = values2.at(RadSystem::radEnergy_index)[i]; double Erad_t = 0.0; - for (int g = 0; g < Physics_Traits::nGroups; ++g) { - Erad_t += values2.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_t += values2.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; } const auto Trad_t = std::pow(Erad_t / a_rad, 1. / 4.); - const auto rho_t = values2.at(RadSystem::gasDensity_index)[i]; - const auto v_t = values2.at(RadSystem::x1GasMomentum_index)[i] / rho_t; - const auto Egas = values2.at(RadSystem::gasInternalEnergy_index)[i]; + const auto rho_t = values2.at(RadSystem::gasDensity_index)[i]; + const auto v_t = values2.at(RadSystem::x1GasMomentum_index)[i] / rho_t; + const auto Egas = values2.at(RadSystem::gasInternalEnergy_index)[i]; xs2.at(i) = x - drift; rhogas2.at(index_) = rho_t; Trad2.at(index_) = Trad_t; - Tgas2.at(index_) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); - Vgas2.at(index_) = 1e-5 * (v_t - v0_adv); + Tgas2.at(index_) = quokka::EOS::ComputeTgasFromEint(rho_t, Egas); + Vgas2.at(index_) = 1e-5 * (v_t - v0); } // END OF PROBLEM 2 @@ -457,7 +411,7 @@ auto problem_main() -> int err_norm += std::abs(Tgas2[i] - Trad[i]); sol_norm += std::abs(Trad[i]) * 3.0; } - const double error_tol = 0.008; + const double error_tol = 0.006; const double rel_error = err_norm / sol_norm; amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; diff --git a/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp b/src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.hpp similarity index 62% rename from src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp rename to src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.hpp index eed311866..ecc3caac8 100644 --- a/src/RadhydroPulseMG/test_radhydro_pulse_MG.hpp +++ b/src/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.hpp @@ -1,6 +1,6 @@ -#ifndef TEST_RADHYDRO_PULSE_MG_HPP_ // NOLINT -#define TEST_RADHYDRO_PULSE_MG_HPP_ -/// \file test_radiation_marshak.hpp +#ifndef TEST_RADHYDRO_PULSE_MG_CONST_HPP_ // NOLINT +#define TEST_RADHYDRO_PULSE_MG_CONST_HPP_ +/// \file test_radhydro_pulse_MG_const_kappa.hpp /// \brief Defines a test problem for radiation in the static diffusion regime. /// @@ -18,4 +18,4 @@ // function definitions -#endif // TEST_RADHYDRO_PULSE_MG_HPP_ +#endif // TEST_RADHYDRO_PULSE_MG_CONST_HPP_ diff --git a/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp index fb63bbc6f..18f630180 100644 --- a/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp +++ b/src/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp @@ -3,6 +3,8 @@ /// #include "test_radhydro_pulse_MG_int.hpp" +#include "AMReX.H" +#include "AMReX_Array.H" #include "AMReX_BC_TYPES.H" #include "AMReX_Print.H" #include "RadhydroSimulation.hpp" @@ -16,32 +18,47 @@ struct MGProblem { struct ExactProblem { }; -// A fixed power law for radiation quantities; for testing purpose only -AMREX_GPU_MANAGED double spec_power = -1.0; // NOLINT -static constexpr bool export_csv = true; - // constexpr int n_groups_ = 2; -// constexpr amrex::GpuArray rad_boundaries_{1e16, 1e18, 1e20}; constexpr int n_groups_ = 4; -constexpr amrex::GpuArray rad_boundaries_{1e16, 1e17, 1e18, 1e19, 1e20}; // constexpr int n_groups_ = 8; -// constexpr amrex::GpuArray rad_boundaries_{1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19, 3.16e19, 1e20}; // constexpr int n_groups_ = 16; -// constexpr amrex::GpuArray -// rad_boundaries_{1.00000000e+16, 1.77827941e+16, 3.16227766e+16, 5.62341325e+16, 1.00000000e+17, 1.77827941e+17, 3.16227766e+17, 5.62341325e+17, 1.00000000e+18, -// 1.77827941e+18, 3.16227766e+18, 5.62341325e+18, 1.00000000e+19, 1.77827941e+19, 3.16227766e+19, 5.62341325e+19, 1.00000000e+20}; constexpr int n_groups_ = -// 32; constexpr amrex::GpuArray -// rad_boundaries_{1.00000000e+16, 1.33352143e+16, 1.77827941e+16, 2.37137371e+16, 3.16227766e+16, 4.21696503e+16, 5.62341325e+16, 7.49894209e+16, 1.00000000e+17, -// 1.33352143e+17, 1.77827941e+17, 2.37137371e+17, 3.16227766e+17, 4.21696503e+17, 5.62341325e+17, 7.49894209e+17, 1.00000000e+18, 1.33352143e+18, 1.77827941e+18, -// 2.37137371e+18, 3.16227766e+18, 4.21696503e+18, 5.62341325e+18, 7.49894209e+18, 1.00000000e+19, 1.33352143e+19, 1.77827941e+19, 2.37137371e+19, 3.16227766e+19, -// 4.21696503e+19, 5.62341325e+19, 7.49894209e+19, 1.00000000e+20}; constexpr int n_groups_ = 64; constexpr amrex::GpuArray -// rad_boundaries_{1.00000000e+16, 1.15478198e+16, 1.33352143e+16, 1.53992653e+16, 1.77827941e+16, 2.05352503e+16, 2.37137371e+16, 2.73841963e+16, 3.16227766e+16, -// 3.65174127e+16, 4.21696503e+16, 4.86967525e+16, 5.62341325e+16, 6.49381632e+16, 7.49894209e+16, 8.65964323e+16, 1.00000000e+17, 1.15478198e+17, 1.33352143e+17, -// 1.53992653e+17, 1.77827941e+17, 2.05352503e+17, 2.37137371e+17, 2.73841963e+17, 3.16227766e+17, 3.65174127e+17, 4.21696503e+17, 4.86967525e+17, 5.62341325e+17, -// 6.49381632e+17, 7.49894209e+17, 8.65964323e+17, 1.00000000e+18, 1.15478198e+18, 1.33352143e+18, 1.53992653e+18, 1.77827941e+18, 2.05352503e+18, 2.37137371e+18, -// 2.73841963e+18, 3.16227766e+18, 3.65174127e+18, 4.21696503e+18, 4.86967525e+18, 5.62341325e+18, 6.49381632e+18, 7.49894209e+18, 8.65964323e+18, 1.00000000e+19, -// 1.15478198e+19, 1.33352143e+19, 1.53992653e+19, 1.77827941e+19, 2.05352503e+19, 2.37137371e+19, 2.73841963e+19, 3.16227766e+19, 3.65174127e+19, 4.21696503e+19, -// 4.86967525e+19, 5.62341325e+19, 6.49381632e+19, 7.49894209e+19, 8.65964323e+19, 1.00000000e+20}; +// constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity; +constexpr OpacityModel opacity_model_ = OpacityModel::PPL_opacity_fixed_slope_spectrum; +// constexpr OpacityModel opacity_model_ = OpacityModel::PPL_opacity_full_spectrum; + +constexpr amrex::GpuArray rad_boundaries_ = []() constexpr { + if constexpr (n_groups_ == 2) { + return amrex::GpuArray{1e15, 1e17, 1e19}; + } else if constexpr (n_groups_ == 4) { + return amrex::GpuArray{1e15, 1e16, 1e17, 1e18, 1e19}; + } else if constexpr (n_groups_ == 8) { + return amrex::GpuArray{1e15, 3.16e15, 1e16, 3.16e16, 1e17, 3.16e17, 1e18, 3.16e18, 1e19}; + } else if constexpr (n_groups_ == 16) { + return amrex::GpuArray{1.00000000e+15, 1.77827941e+15, 3.16227766e+15, 5.62341325e+15, 1.00000000e+16, 1.77827941e+16, + 3.16227766e+16, 5.62341325e+16, 1.00000000e+17, 1.77827941e+17, 3.16227766e+17, 5.62341325e+17, + 1.00000000e+18, 1.77827941e+18, 3.16227766e+18, 5.62341325e+18, 1.00000000e+19}; + } else if constexpr (n_groups_ == 32) { + return amrex::GpuArray{1.00000000e+15, 1.33352143e+15, 1.77827941e+15, 2.37137371e+15, 3.16227766e+15, 4.21696503e+15, + 5.62341325e+15, 7.49894209e+15, 1.00000000e+16, 1.33352143e+16, 1.77827941e+16, 2.37137371e+16, + 3.16227766e+16, 4.21696503e+16, 5.62341325e+16, 7.49894209e+16, 1.00000000e+17, 1.33352143e+17, + 1.77827941e+17, 2.37137371e+17, 3.16227766e+17, 4.21696503e+17, 5.62341325e+17, 7.49894209e+17, + 1.00000000e+18, 1.33352143e+18, 1.77827941e+18, 2.37137371e+18, 3.16227766e+18, 4.21696503e+18, + 5.62341325e+18, 7.49894209e+18, 1.00000000e+19}; + } else if constexpr (n_groups_ == 64) { + return amrex::GpuArray{ + 1.00000000e+15, 1.15478198e+15, 1.33352143e+15, 1.53992653e+15, 1.77827941e+15, 2.05352503e+15, 2.37137371e+15, 2.73841963e+15, + 3.16227766e+15, 3.65174127e+15, 4.21696503e+15, 4.86967525e+15, 5.62341325e+15, 6.49381632e+15, 7.49894209e+15, 8.65964323e+15, + 1.00000000e+16, 1.15478198e+16, 1.33352143e+16, 1.53992653e+16, 1.77827941e+16, 2.05352503e+16, 2.37137371e+16, 2.73841963e+16, + 3.16227766e+16, 3.65174127e+16, 4.21696503e+16, 4.86967525e+16, 5.62341325e+16, 6.49381632e+16, 7.49894209e+16, 8.65964323e+16, + 1.00000000e+17, 1.15478198e+17, 1.33352143e+17, 1.53992653e+17, 1.77827941e+17, 2.05352503e+17, 2.37137371e+17, 2.73841963e+17, + 3.16227766e+17, 3.65174127e+17, 4.21696503e+17, 4.86967525e+17, 5.62341325e+17, 6.49381632e+17, 7.49894209e+17, 8.65964323e+17, + 1.00000000e+18, 1.15478198e+18, 1.33352143e+18, 1.53992653e+18, 1.77827941e+18, 2.05352503e+18, 2.37137371e+18, 2.73841963e+18, + 3.16227766e+18, 3.65174127e+18, 4.21696503e+18, 4.86967525e+18, 5.62341325e+18, 6.49381632e+18, 7.49894209e+18, 8.65964323e+18, + 1.00000000e+19}; + } +}(); + +static constexpr bool export_csv = true; constexpr double T0 = 1.0e7; // K (temperature) constexpr double T1 = 2.0e7; // K (temperature) @@ -56,10 +73,13 @@ constexpr double h_planck = C::hplanck; constexpr double k_B = C::k_B; // static diffusion: (for single group) tau = 2e3, beta = 3e-5, beta tau = 6e-2 -constexpr double kappa0 = 180.; // cm^2 g^-1 -constexpr double v0_adv = 1.0e6; // advecting pulse -constexpr double max_time = 4.8e-5; // max_time = 2 * width / v1; -constexpr int64_t max_timesteps = 30; // to make 3D test run fast on GPUs +constexpr double kappa0 = 180.; // cm^2 g^-1 +constexpr double scaleup = 1.; +constexpr double v0_adv = 1.0e6; // advecting pulse +constexpr double max_time = 4.8e-5; // max_time = 2 * width / v1; +// constexpr double max_time = 2e-5; // max_time = 2 * width / v1; +constexpr int64_t max_timesteps = 1e2; // to make 3D test run fast on GPUs +// constexpr int64_t max_timesteps = 1e8; // full run // dynamic diffusion: tau = 2e4, beta = 3e-3, beta tau = 60 // constexpr double kappa0 = 1000.; // cm^2 g^-1 @@ -111,7 +131,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = h_planck; static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr int beta_order = 1; - static constexpr OpacityModel opacity_model = OpacityModel::piecewisePowerLaw; + static constexpr OpacityModel opacity_model = opacity_model_; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -120,22 +140,9 @@ template <> struct RadSystem_Traits { static constexpr double Erad_floor = erad_floor; static constexpr bool compute_v_over_c_terms = true; static constexpr int beta_order = 1; - static constexpr OpacityModel opacity_model = OpacityModel::user; + static constexpr OpacityModel opacity_model = OpacityModel::single_group; }; -template <> -template -AMREX_GPU_HOST_DEVICE auto -RadSystem::ComputeRadQuantityExponents(ArrayType const & /*quant*/, - amrex::GpuArray const & /*boundaries*/) -> amrex::GpuArray -{ - amrex::GpuArray exponents{}; - for (int g = 0; g < nGroups_; ++g) { - exponents[g] = spec_power; - } - return exponents; -} - AMREX_GPU_HOST_DEVICE auto compute_initial_Tgas(const double x) -> double { @@ -158,31 +165,47 @@ auto compute_kappa(const double nu, const double Tgas) -> double // cm^-1 auto T_ = Tgas / T_ref; auto nu_ = nu / nu_ref; - return kappa0 * std::pow(T_, -0.5) * std::pow(nu_, -3) * (1.0 - std::exp(-coeff_ * nu_ / T_)); + return scaleup * kappa0 * std::pow(T_, -0.5) * std::pow(nu_, -3) * (1.0 - std::exp(-coeff_ * nu_ / T_)); } template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray const rad_boundaries, const double rho, - const double Tgas) -> amrex::GpuArray, 2> + const double Tgas) -> amrex::GpuArray, 2> { - amrex::GpuArray exponents{}; - amrex::GpuArray kappa_lower{}; - for (int g = 0; g < nGroups_; ++g) { - auto kappa_up = compute_kappa(rad_boundaries[g + 1], Tgas); - auto kappa_down = compute_kappa(rad_boundaries[g], Tgas); - exponents[g] = std::log(kappa_up / kappa_down) / std::log(rad_boundaries[g + 1] / rad_boundaries[g]); - kappa_lower[g] = kappa_down / rho; + amrex::GpuArray exponents{}; + amrex::GpuArray kappa_lower{}; + + for (int g = 0; g < nGroups_ + 1; ++g) { + if constexpr (RadSystem_Traits::opacity_model == OpacityModel::piecewise_constant_opacity) { + exponents[g] = 0.0; + if (g < n_groups_) { + auto nu_center = std::sqrt(rad_boundaries[g] * rad_boundaries[g + 1]); + kappa_lower[g] = compute_kappa(nu_center, Tgas) / rho; + } + // note kappa_lower[nGroups_] is not used + } else { + if (g == n_groups_) { + exponents[g] = 0.; + kappa_lower[g] = compute_kappa(rad_boundaries[g], Tgas) / rho; + } else { + auto kappa_up = compute_kappa(rad_boundaries[g + 1], Tgas); + auto kappa_down = compute_kappa(rad_boundaries[g], Tgas); + exponents[g] = std::log(kappa_up / kappa_down) / std::log(rad_boundaries[g + 1] / rad_boundaries[g]); + kappa_lower[g] = kappa_down / rho; + } + } AMREX_ASSERT(!std::isnan(exponents[g])); AMREX_ASSERT(kappa_lower[g] >= 0.); } - amrex::GpuArray, 2> const exponents_and_values{exponents, kappa_lower}; + + amrex::GpuArray, 2> const exponents_and_values{exponents, kappa_lower}; return exponents_and_values; } template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray { - const double sigma = 3063.96 * std::pow(Tgas / T0, -3.5); + const double sigma = scaleup * 3063.96 * std::pow(Tgas / T0, -3.5); quokka::valarray kappaPVec{}; kappaPVec.fillin(sigma / rho); return kappaPVec; @@ -191,7 +214,7 @@ template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpa template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray { - const double sigma = 101.248 * std::pow(Tgas / T0, -3.5); + const double sigma = scaleup * 101.248 * std::pow(Tgas / T0, -3.5); quokka::valarray kappaPVec{}; kappaPVec.fillin(sigma / rho); return kappaPVec; @@ -219,10 +242,13 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokk const double v0 = v0_adv; auto Erad_g = RadSystem::ComputeThermalRadiation(Trad, radBoundaries_g); + auto Frad_g = RadSystem::ComputeFluxInDiffusionLimit(radBoundaries_g, Trad, v0); for (int g = 0; g < Physics_Traits::nGroups; ++g) { state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g]; - state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + // OLD, correct if you ignore the (delta nu B) term + // state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 4. / 3. * v0 * Erad_g[g]; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frad_g[g]; state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; } @@ -280,7 +306,6 @@ auto problem_main() -> int const double max_dt = 1e-3; // t_cr = 2 cm / cs = 7e-8 s amrex::ParmParse const pp("rad"); - pp.query("spec_power", spec_power); // Boundary conditions constexpr int nvars = RadSystem::nvar_; @@ -444,10 +469,30 @@ auto problem_main() -> int err_norm += std::abs(Tgas[i] - Tgas2[i]); sol_norm += std::abs(Tgas2[i]); } - const double error_tol = 0.008; + const double error_tol = 0.02; const double rel_error = err_norm / sol_norm; amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + // symmetry check + double symm_err = 0.; + double symm_norm = 0.; + const double symm_err_tol = 0.02; + for (size_t i = 0; i < xs.size(); ++i) { + symm_err += std::abs(Tgas[i] - Tgas[xs.size() - 1 - i]); + symm_norm += std::abs(Tgas[i]); + } + const double symm_rel_error_1 = symm_err / symm_norm; + amrex::Print() << "Symmetry L1 error norm of the MG pulse = " << symm_rel_error_1 << std::endl; + + symm_err = 0.; + symm_norm = 0.; + for (size_t i = 0; i < xs2.size(); ++i) { + symm_err += std::abs(Tgas2[i] - Tgas2[xs2.size() - 1 - i]); + symm_norm += std::abs(Tgas2[i]); + } + const double symm_rel_error_2 = symm_err / symm_norm; + amrex::Print() << "Symmetry L1 error norm of the exact (grey) pulse = " << symm_rel_error_2 << std::endl; + #ifdef HAVE_PYTHON // plot temperature matplotlibcpp::clf(); @@ -511,6 +556,7 @@ auto problem_main() -> int matplotlibcpp::plot(xs2, Vgas2, vgas_args); matplotlibcpp::xlabel("length x (cm)"); matplotlibcpp::ylabel("gas velocity (km s^-1)"); + matplotlibcpp::ylim(-5., 5.); matplotlibcpp::legend(); matplotlibcpp::title(fmt::format("nGroups = {}, time t = {:.4g}", n_groups_, sim.tNew_[0])); matplotlibcpp::tight_layout(); @@ -547,7 +593,8 @@ auto problem_main() -> int // Cleanup and exit int status = 0; - if ((rel_error > error_tol) || std::isnan(rel_error)) { + if ((rel_error > error_tol) || std::isnan(rel_error) || (symm_rel_error_1 > symm_err_tol) || (symm_rel_error_2 > symm_err_tol) || + std::isnan(symm_rel_error_1) || std::isnan(symm_rel_error_2)) { status = 1; } return status; diff --git a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index e37bb1408..0035e284f 100644 --- a/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -63,7 +63,9 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.00000000e+15, 1.00000000e+16, 1.00000000e+17, 1.00000000e+18, 1.00000000e+19, 1.00000000e+20}; static constexpr int beta_order = 1; - static constexpr OpacityModel opacity_model = OpacityModel::piecewisePowerLaw; + // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; + static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; + // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; }; template <> struct quokka::EOS_Traits { @@ -75,13 +77,11 @@ template <> struct quokka::EOS_Traits { template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, - const double /*Tgas*/) -> amrex::GpuArray, 2> + const double /*Tgas*/) -> amrex::GpuArray, 2> { - amrex::GpuArray, 2> exponents_and_values{}; - for (int i = 0; i < nGroups_; ++i) { + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_ + 1; ++i) { exponents_and_values[0][i] = 0.0; - } - for (int i = 0; i < nGroups_; ++i) { exponents_and_values[1][i] = kappa / rho; } return exponents_and_values; @@ -260,6 +260,8 @@ auto problem_main() -> int // Problem initialization RadhydroSimulation sim(BCs_cc); + sim.radiationReconstructionOrder_ = 3; // PPM + sim.reconstructionOrder_ = 3; // PPM sim.cflNumber_ = CFL_number; sim.radiationCflNumber_ = CFL_number; sim.maxTimesteps_ = max_timesteps; diff --git a/src/radiation_system.hpp b/src/radiation_system.hpp index 102ad115e..0a18d4d65 100644 --- a/src/radiation_system.hpp +++ b/src/radiation_system.hpp @@ -20,6 +20,7 @@ #include "AMReX_BLassert.H" #include "AMReX_GpuQualifiers.H" #include "AMReX_IParser_Y.H" +#include "AMReX_IntVect.H" #include "AMReX_REAL.H" // internal headers @@ -32,8 +33,20 @@ // Hyper parameters of the radiation solver +static constexpr bool include_delta_B = true; +static constexpr bool use_diffuse_flux_mean_opacity = true; +static constexpr bool special_edge_bin_slopes = false; // Use 2 and -4 as the slopes for the first and last bins, respectively +static constexpr bool force_rad_floor_in_iteration = false; // force radiation energy density to be positive (and above the floor value) in the Newton iteration + +static constexpr bool use_diffuse_flux_mean_opacity_incl_kappaE = true; +static const int max_ite_to_update_alpha_E = 5; // Apply to the PPL_opacity_full_spectrum only. Only update alpha_E for the first max_ite_to_update_alpha_E +// iterations of the Newton iteration + static constexpr bool include_work_term_in_source = true; static constexpr bool use_D_as_base = true; +static const bool PPL_free_slope_st_total = false; // PPL with free slopes for all, but subject to the constraint sum_g alpha_g B_g = - sum_g B_g. Not working + // well -- Newton iteration convergence issue. +static const bool use_wavespeed_correction = false; // Optional: include a wavespeed correction term in the radiation flux to suppress instability // Time integration scheme // IMEX PD-ARS @@ -48,14 +61,13 @@ static constexpr double c_light_cgs_ = C::c_light; // cgs static constexpr double radiation_constant_cgs_ = C::a_rad; // cgs static constexpr double inf = std::numeric_limits::max(); -// Optional: include a wavespeed correction term in the radiation flux to suppress instability -static const bool use_wavespeed_correction = false; - // enum for opacity_model enum class OpacityModel { - user = 0, // user-defined opacity for each group, given as a function of density and temperature. - piecewisePowerLaw // piecewise power-law opacity model with piecewise power-law fitting to a user-defined opacity function and on-the-fly piecewise - // power-law fitting to radiation energy density and flux. + single_group = 0, // user-defined opacity for each group, given as a function of density and temperature. + piecewise_constant_opacity, + PPL_opacity_fixed_slope_spectrum, + PPL_opacity_full_spectrum // piecewise power-law opacity model with piecewise power-law fitting to a user-defined opacity function and on-the-fly + // piecewise power-law fitting to radiation energy density and flux. }; // this struct is specialized by the user application code @@ -68,7 +80,7 @@ template struct RadSystem_Traits { static constexpr double energy_unit = C::ev2erg; static constexpr amrex::GpuArray::nGroups + 1> radBoundaries = {0., inf}; static constexpr double beta_order = 1; - static constexpr OpacityModel opacity_model = OpacityModel::user; + static constexpr OpacityModel opacity_model = OpacityModel::single_group; }; // A struct to hold the results of the ComputeRadPressure function. @@ -148,10 +160,13 @@ template class RadSystem : public HyperbolicSystem::value) { return RadSystem_Traits::opacity_model; } else { - return OpacityModel::user; + return OpacityModel::single_group; } }(); + static_assert(!(nGroups_ < 3 && opacity_model_ == OpacityModel::PPL_opacity_full_spectrum), + "PPL_opacity_full_spectrum requires at least 3 photon groups."); // NOLINT + static constexpr double mean_molecular_mass_ = quokka::EOS_Traits::mean_molecular_mass; static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; static constexpr double gamma_ = quokka::EOS_Traits::gamma; @@ -195,13 +210,26 @@ template class RadSystem : public HyperbolicSystem quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeEnergyMeanOpacity(double rho, double Tgas) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto DefineOpacityExponentsAndLowerValues(amrex::GpuArray rad_boundaries, double rho, - double Tgas) -> amrex::GpuArray, 2>; - AMREX_GPU_HOST_DEVICE static auto ComputeGroupMeanOpacity(amrex::GpuArray, 2> kappa_expo_and_lower_value, - amrex::GpuArray radBoundaryRatios, - amrex::GpuArray alpha_quant) -> quokka::valarray; - AMREX_GPU_HOST_DEVICE static auto ComputePlanckOpacityTempDerivative(double rho, double Tgas) -> quokka::valarray; + double Tgas) -> amrex::GpuArray, 2>; + AMREX_GPU_HOST_DEVICE static auto ComputeGroupMeanOpacity(amrex::GpuArray, 2> const &kappa_expo_and_lower_value, + amrex::GpuArray const &radBoundaryRatios, + amrex::GpuArray const &alpha_quant) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto + ComputeBinCenterOpacity(amrex::GpuArray rad_boundaries, + amrex::GpuArray, 2> kappa_expo_and_lower_value) -> quokka::valarray; + // AMREX_GPU_HOST_DEVICE static auto + // ComputeGroupMeanOpacityWithMinusOneSlope(amrex::GpuArray, 2> kappa_expo_and_lower_value, + // amrex::GpuArray radBoundaryRatios) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeEintFromEgas(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Etot) -> double; AMREX_GPU_HOST_DEVICE static auto ComputeEgasFromEint(double density, double X1GasMom, double X2GasMom, double X3GasMom, double Eint) -> double; + AMREX_GPU_HOST_DEVICE static auto PlanckFunction(double nu, double T) -> double; + AMREX_GPU_HOST_DEVICE static auto + ComputeDiffusionFluxMeanOpacity(quokka::valarray kappaPVec, quokka::valarray kappaEVec, + quokka::valarray fourPiBoverC, amrex::GpuArray delta_nu_kappa_B_at_edge, + amrex::GpuArray delta_nu_B_at_edge, + amrex::GpuArray kappa_slope) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto ComputeFluxInDiffusionLimit(amrex::GpuArray rad_boundaries, double T, + double vel) -> amrex::GpuArray; template AMREX_GPU_HOST_DEVICE static auto @@ -248,7 +276,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(am { quokka::valarray radEnergyFractions{}; if constexpr (nGroups_ == 1) { - // TODO(CCH): allow the total radEnergyFraction to be smaller than 1. One usage case is to allow, say, a single group from 13.6 eV to 24.6 eV. + // TODO(CCH): allow the total radEnergyFraction to be smaller than 1. One usage case is to allow, say, a single group representing IR radiation. radEnergyFractions[0] = 1.0; return radEnergyFractions; } else { @@ -886,11 +914,11 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const doub template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray { - quokka::valarray kappaFVec{}; + quokka::valarray kappa{}; for (int g = 0; g < nGroups_; ++g) { - kappaFVec[g] = NAN; + kappa[g] = NAN; } - return kappaFVec; + return kappa; } template @@ -902,9 +930,9 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEnergyMeanOpacity(const template AMREX_GPU_HOST_DEVICE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double /*rho*/, - const double /*Tgas*/) -> amrex::GpuArray, 2> + const double /*Tgas*/) -> amrex::GpuArray, 2> { - amrex::GpuArray, 2> exponents_and_values{}; + amrex::GpuArray, 2> exponents_and_values{}; return exponents_and_values; } @@ -915,6 +943,9 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeRadQuantityExponents(Arr { // Compute the exponents for the radiation energy density, radiation flux, radiation pressure, or Planck function. + // static assert nGroups_ is at least 3 + static_assert(nGroups_ >= 3); + // Note: Could save some memory by using bin_center_previous and bin_center_current amrex::GpuArray bin_center{}; amrex::GpuArray quant_mean{}; @@ -939,37 +970,83 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeRadQuantityExponents(Arr AMREX_ASSERT(!std::isnan(logslopes[g - 1])); } } + for (int g = 0; g < nGroups_; ++g) { - if (g == 0 || g == nGroups_ - 1) { - exponents[g] = 0.0; + if (g == 0) { + if constexpr (!special_edge_bin_slopes) { + exponents[g] = -1.0; + } else { + exponents[g] = 2.0; + } + } else if (g == nGroups_ - 1) { + if constexpr (!special_edge_bin_slopes) { + exponents[g] = -1.0; + } else { + exponents[g] = -4.0; + } } else { exponents[g] = minmod_func(logslopes[g - 1], logslopes[g]); } AMREX_ASSERT(!std::isnan(exponents[g])); - AMREX_ASSERT(std::abs(exponents[g]) < 100); } + if constexpr (PPL_free_slope_st_total) { + int peak_idx = 0; // index of the peak of logslopes + for (; peak_idx < nGroups_; ++peak_idx) { + if (peak_idx == nGroups_ - 1) { + peak_idx += 0; + break; + } + if (exponents[peak_idx] >= 0.0 && exponents[peak_idx + 1] < 0.0) { + break; + } + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(peak_idx < nGroups_ - 1, + "Peak index not found. Here peak_index is the index at which the exponent changes its sign."); + double quant_sum = 0.0; + double part_sum = 0.0; + for (int g = 0; g < nGroups_; ++g) { + quant_sum += quant[g]; + if (g == peak_idx) { + continue; + } + part_sum += exponents[g] * quant[g]; + } + if (quant[peak_idx] > 0.0 && quant_sum > 0.0) { + exponents[peak_idx] = (-quant_sum - part_sum) / quant[peak_idx]; + AMREX_ASSERT(!std::isnan(exponents[peak_idx])); + } + } return exponents; } template -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeGroupMeanOpacity(amrex::GpuArray, 2> kappa_expo_and_lower_value, - amrex::GpuArray radBoundaryRatios, - amrex::GpuArray alpha_quant) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto +RadSystem::ComputeGroupMeanOpacity(amrex::GpuArray, 2> const &kappa_expo_and_lower_value, + amrex::GpuArray const &radBoundaryRatios, + amrex::GpuArray const &alpha_quant) -> quokka::valarray { - amrex::GpuArray const &alpha_kappa = kappa_expo_and_lower_value[0]; - amrex::GpuArray const &kappa_lower = kappa_expo_and_lower_value[1]; + amrex::GpuArray const &alpha_kappa = kappa_expo_and_lower_value[0]; + amrex::GpuArray const &kappa_lower = kappa_expo_and_lower_value[1]; quokka::valarray kappa{}; for (int g = 0; g < nGroups_; ++g) { double alpha = alpha_quant[g] + 1.0; + if (alpha > 100.) { + kappa[g] = kappa_lower[g] * std::pow(radBoundaryRatios[g], kappa_expo_and_lower_value[0][g]); + continue; + } + if (alpha < -100.) { + kappa[g] = kappa_lower[g]; + continue; + } double part1 = 0.0; if (std::abs(alpha) < 1e-8) { part1 = std::log(radBoundaryRatios[g]); } else { part1 = (std::pow(radBoundaryRatios[g], alpha) - 1.0) / alpha; } - alpha = alpha_quant[g] + alpha_kappa[g] + 1.0; + alpha += alpha_kappa[g]; double part2 = 0.0; if (std::abs(alpha) < 1e-8) { part2 = std::log(radBoundaryRatios[g]); @@ -982,15 +1059,6 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeGroupMeanOpacity(amrex:: return kappa; } -template -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacityTempDerivative(const double /* rho */, - const double /* Tgas */) -> quokka::valarray -{ - quokka::valarray kappa{}; - kappa.fillin(0.0); - return kappa; -} - template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEintFromEgas(const double density, const double X1GasMom, const double X2GasMom, const double X3GasMom, const double Etot) -> double @@ -1012,6 +1080,76 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeEgasFromEint(const doubl return Etot; } +template AMREX_GPU_HOST_DEVICE auto RadSystem::PlanckFunction(const double nu, const double T) -> double +{ + // returns 4 pi B(nu) / c + double const coeff = RadSystem_Traits::energy_unit / (boltzmann_constant_ * T); + double const x = coeff * nu; + if (x > 100.) { + return 0.0; + } + double planck_integral = NAN; + if (x <= 1.0e-10) { + // Taylor series + planck_integral = x * x - x * x * x / 2.; + } else { + planck_integral = std::pow(x, 3) / (std::exp(x) - 1.0); + } + return coeff / (std::pow(PI, 4) / 15.0) * (radiation_constant_ * std::pow(T, 4)) * planck_integral; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDiffusionFluxMeanOpacity( + const quokka::valarray kappaPVec, const quokka::valarray kappaEVec, + const quokka::valarray fourPiBoverC, const amrex::GpuArray delta_nu_kappa_B_at_edge, + const amrex::GpuArray delta_nu_B_at_edge, const amrex::GpuArray kappa_slope) -> quokka::valarray +{ + quokka::valarray kappaF{}; + for (int g = 0; g < nGroups_; ++g) { + // kappaF[g] = 4. / 3. * kappaPVec[g] * fourPiBoverC[g] + 1. / 3. * kappa_slope[g] * kappaPVec[g] * fourPiBoverC[g] - 1. / 3. * + // delta_nu_kappa_B_at_edge[g]; + kappaF[g] = (kappaPVec[g] + 1. / 3. * kappaEVec[g]) * fourPiBoverC[g] + + 1. / 3. * (kappa_slope[g] * kappaEVec[g] * fourPiBoverC[g] - delta_nu_kappa_B_at_edge[g]); + AMREX_ASSERT(kappaF[g] > 0.0); + auto const denom = 4. / 3. * fourPiBoverC[g] - 1. / 3. * delta_nu_B_at_edge[g]; + AMREX_ASSERT(denom > 0.0); + kappaF[g] /= denom; + } + return kappaF; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeBinCenterOpacity(amrex::GpuArray rad_boundaries, + amrex::GpuArray, 2> kappa_expo_and_lower_value) + -> quokka::valarray +{ + quokka::valarray kappa_center{}; + for (int g = 0; g < nGroups_; ++g) { + kappa_center[g] = + kappa_expo_and_lower_value[1][g] * std::pow(rad_boundaries[g + 1] / rad_boundaries[g], 0.5 * kappa_expo_and_lower_value[0][g]); + } + return kappa_center; +} + +template +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxInDiffusionLimit(const amrex::GpuArray rad_boundaries, const double T, + const double vel) -> amrex::GpuArray +{ + double const coeff = RadSystem_Traits::energy_unit / (boltzmann_constant_ * T); + amrex::GpuArray edge_values{}; + amrex::GpuArray flux{}; + for (int g = 0; g < nGroups_ + 1; ++g) { + auto x = coeff * rad_boundaries[g]; + edge_values[g] = 4. / 3. * integrate_planck_from_0_to_x(x) - 1. / 3. * x * (std::pow(x, 3) / (std::exp(x) - 1.0)) / gInf; + // test: reproduce the Planck function + // edge_values[g] = 4. / 3. * integrate_planck_from_0_to_x(x); + } + for (int g = 0; g < nGroups_; ++g) { + flux[g] = vel * radiation_constant_ * std::pow(T, 4) * (edge_values[g + 1] - edge_values[g]); + } + return flux; +} + template void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt_radiation, const int stage) @@ -1026,7 +1164,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::GpuArray radBoundaries_g = radBoundaries_; amrex::GpuArray radBoundaryRatios{}; if constexpr (nGroups_ > 1) { - if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + if constexpr (static_cast(opacity_model_) > 0) { for (int g = 0; g < nGroups_; ++g) { radBoundaryRatios[g] = radBoundaries_g[g + 1] / radBoundaries_g[g]; } @@ -1080,25 +1218,23 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray kappaPVec{}; quokka::valarray kappaEVec{}; quokka::valarray kappaFVec{}; - amrex::GpuArray, 2> kappa_expo_and_lower_value{}; + amrex::GpuArray, 2> kappa_expo_and_lower_value{}; amrex::GpuArray alpha_B{}; amrex::GpuArray alpha_E{}; - amrex::GpuArray alpha_F{}; quokka::valarray kappaPoverE{}; quokka::valarray tau0{}; // optical depth across c * dt at old state quokka::valarray tau{}; // optical depth across c * dt at new state quokka::valarray D{}; // D = S / tau0 quokka::valarray work{}; quokka::valarray work_prev{}; - amrex::GpuArray, 3> frad{}; amrex::GpuArray dMomentum{}; amrex::GpuArray, 3> Frad_t1{}; + amrex::GpuArray delta_nu_kappa_B_at_edge{}; + amrex::GpuArray delta_nu_B_at_edge{}; work.fillin(0.0); work_prev.fillin(0.0); - EradVec_guess = Erad0Vec; - if constexpr (gamma_ != 1.0) { Egas0 = ComputeEintFromEgas(rho, x1GasMom0, x2GasMom0, x3GasMom0, Egastot0); Etot0 = Egas0 + (c / chat) * (Erad0 + sum(Src)); @@ -1114,6 +1250,23 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } + // define a list of alpha_quant for the model PPL_opacity_fixed_slope_spectrum + amrex::GpuArray alpha_quant_minus_one{}; + if constexpr ((opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) || + (gamma_ == 1.0 && opacity_model_ == OpacityModel::PPL_opacity_full_spectrum)) { + if constexpr (!special_edge_bin_slopes) { + for (int g = 0; g < nGroups_; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } else { + alpha_quant_minus_one[0] = 2.0; + alpha_quant_minus_one[nGroups_ - 1] = -4.0; + for (int g = 1; g < nGroups_ - 1; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + } + } + amrex::Real gas_update_factor = 1.0; if (stage == 1) { gas_update_factor = IMEX_a32; @@ -1124,7 +1277,10 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (; ite < max_ite; ++ite) { quokka::valarray Rvec{}; + EradVec_guess = Erad0Vec; + if constexpr (gamma_ != 1.0) { + Egas_guess = Egas0; Ekin0 = Egastot0 - Egas0; AMREX_ASSERT(min(Src) >= 0.0); @@ -1169,21 +1325,57 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // dF_{D,i} / dE_g = 1 / (chat * C_v) * (kappa_{P,i} / kappa_{E,i}) * d/dT (4 \pi B_i) // dF_{D,i} / dD_i = - (1 / (chat * dt * rho * kappa_{E,i}) + 1) * tau0_i = - ((1 / tau_i)(kappa_Pi / kappa_Ei) + 1) * tau0_i - Egas_guess = Egas0; - T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas0, massScalars); AMREX_ASSERT(T_gas >= 0.); fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - if constexpr (opacity_model_ == OpacityModel::user) { + if constexpr (opacity_model_ == OpacityModel::single_group) { kappaPVec = ComputePlanckOpacity(rho, T_gas); kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = kappa_expo_and_lower_value[1][g]; + kappaEVec[g] = kappa_expo_and_lower_value[1][g]; + kappaFVec[g] = kappa_expo_and_lower_value[1][g]; + } + } else { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + for (int g = 0; g < nGroups_; ++g) { + auto const nu_L = radBoundaries_g_copy[g]; + auto const nu_R = radBoundaries_g_copy[g + 1]; + auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const kappa_L = kappa_expo_and_lower_value[1][g]; + auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); + delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; + delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; + } + if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + kappaEVec = kappaPVec; + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = + ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, + delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + } else { + kappaFVec = kappaPVec; + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = + ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, + delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + } else { // fall back to bin-center opacity. For simplicity we don't do PPL on the flux mean opacity because it + // requires fitting to each components of the flux. + kappaFVec = ComputeBinCenterOpacity(radBoundaries_g_copy, kappa_expo_and_lower_value); + } + } } AMREX_ASSERT(!kappaPVec.hasnan()); AMREX_ASSERT(!kappaEVec.hasnan()); @@ -1201,35 +1393,33 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // compute the work term at the old state // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); if (ite == 0) { - if constexpr (opacity_model_ == OpacityModel::user) { + if constexpr (opacity_model_ == OpacityModel::single_group) { for (int g = 0; g < nGroups_; ++g) { - // work[g] = dt * chat * rho * kappaPVec[g] * (Erad0Vec[g] - fourPiBoverC[g]); const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); // work = v * F * chi work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * - (2.0 * kappaEVec[g] - kappaFVec[g]); - work[g] *= chat / (c * c) * lorentz_factor_v * dt; + (2.0 * kappaEVec[g] - kappaFVec[g]) * chat / (c * c) * lorentz_factor_v * dt; } - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { for (int g = 0; g < nGroups_; ++g) { - frad[0][g] = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); - frad[1][g] = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); - frad[2][g] = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); - work[g] = 0.0; - } - for (int n = 0; n < 3; ++n) { - alpha_F = ComputeRadQuantityExponents(frad[n], radBoundaries_g_copy); - kappaFVec = - ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); - for (int g = 0; g < nGroups_; ++g) { - work[g] += - (kappa_expo_and_lower_value[0][g] + 1.0) * gasMtm0[n] * kappaFVec[g] * frad[n][g]; - } + const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + // work = v * F * chi + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * kappaFVec[g] * chat / + (c * c) * dt; } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { for (int g = 0; g < nGroups_; ++g) { - work[g] *= chat / (c * c) * dt; + const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + // work = v * F * chi + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; } } } @@ -1256,47 +1446,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray deltaD{}; quokka::valarray F_D{}; + tau = tau0; const double resid_tol = 1.0e-11; // 1.0e-15; const int maxIter = 400; int n = 0; for (; n < maxIter; ++n) { - // compute material temperature - T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); - AMREX_ASSERT(T_gas >= 0.); - // compute opacity, emissivity - fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - - if constexpr (opacity_model_ == OpacityModel::user) { - kappaPVec = ComputePlanckOpacity(rho, T_gas); - kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); - } - AMREX_ASSERT(!kappaPVec.hasnan()); - AMREX_ASSERT(!kappaEVec.hasnan()); - - for (int g = 0; g < nGroups_; ++g) { - if (kappaEVec[g] > 0.0) { - kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; - } else { - kappaPoverE[g] = 1.0; - } - } - - tau = dt * rho * kappaEVec * chat * lorentz_factor; - if constexpr (use_D_as_base) { - Rvec = tau0 * D; - } - for (int g = 0; g < nGroups_; ++g) { - // If tau = 0.0, Erad_guess shouldn't change - if (tau[g] > 0.0) { - EradVec_guess[g] = kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work[g]) / tau[g]); - } - } // F_G = Egas_guess - Egas0 + (c / chat) * sum(Rvec); F_G = Egas_guess - Egas0; F_D = EradVec_guess - Erad0Vec - (Rvec + Src); @@ -1308,6 +1462,13 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } + if (n > 200) { +#ifndef NDEBUG + std::cout << "n = " << n << ", F_G = " << F_G << ", F_D_abs_sum = " << F_D_abs_sum + << ", F_D_abs_sum / Etot0 = " << F_D_abs_sum / Etot0 << std::endl; +#endif + } + // check relative convergence of the residuals if ((std::abs(F_G / Etot0) < resid_tol) && ((c / chat) * F_D_abs_sum / Etot0 < resid_tol)) { break; @@ -1349,6 +1510,63 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne Rvec += deltaD; } + // compute material temperature + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + // compute opacity, emissivity + fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaPVec = ComputePlanckOpacity(rho, T_gas); + kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = kappa_expo_and_lower_value[1][g]; + kappaEVec[g] = kappa_expo_and_lower_value[1][g]; + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + // kappaPVec = ComputeGroupMeanOpacityWithMinusOneSlope(kappa_expo_and_lower_value, radBoundaryRatios_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + kappaEVec = kappaPVec; + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + if (n < max_ite_to_update_alpha_E) { + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); + } + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + } + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); + + for (int g = 0; g < nGroups_; ++g) { + if (kappaEVec[g] > 0.0) { + kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; + } else { + kappaPoverE[g] = 1.0; + } + } + + tau = dt * rho * kappaEVec * chat * lorentz_factor; + if constexpr (use_D_as_base) { + Rvec = tau0 * D; + } + for (int g = 0; g < nGroups_; ++g) { + // If tau = 0.0, Erad_guess shouldn't change + if (tau[g] > 0.0) { + EradVec_guess[g] = kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work[g]) / tau[g]); + if constexpr (force_rad_floor_in_iteration) { + if (EradVec_guess[g] < 0.0) { + Egas_guess -= (c_light_ / c_hat_) * (Erad_floor_ - EradVec_guess[g]); + EradVec_guess[g] = Erad_floor_; + } + } + } + } + // check relative and absolute convergence of E_r // if (std::abs(deltaEgas / Egas_guess) < 1e-7) { // break; @@ -1373,28 +1591,81 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); } - if constexpr (opacity_model_ == OpacityModel::user) { + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note that kappaFVec is used no matter what the value of gamma is if constexpr (gamma_ != 1.0) { kappaPVec = ComputePlanckOpacity(rho, T_gas); kappaEVec = ComputeEnergyMeanOpacity(rho, T_gas); AMREX_ASSERT(!kappaPVec.hasnan()); AMREX_ASSERT(!kappaEVec.hasnan()); } - kappaFVec = ComputeFluxMeanOpacity(rho, T_gas); // note that kappaFVec is used no matter what the value of gamma is - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { + } else { kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); if constexpr (gamma_ != 1.0) { - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); - AMREX_ASSERT(!kappaPVec.hasnan()); - AMREX_ASSERT(!kappaEVec.hasnan()); + for (int g = 0; g < nGroups_; ++g) { + auto const nu_L = radBoundaries_g_copy[g]; + auto const nu_R = radBoundaries_g_copy[g + 1]; + auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + delta_nu_kappa_B_at_edge[g] = kappa_expo_and_lower_value[1][g] * (nu_R * B_R - nu_L * B_L); + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + auto const kappa_L = kappa_expo_and_lower_value[1][g]; + auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); + delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; + delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; + } + } + } + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + for (int g = 0; g < nGroups_; ++g) { + kappaFVec[g] = kappa_expo_and_lower_value[1][g]; + if constexpr (gamma_ != 1.0) { + kappaPVec[g] = kappaFVec[g]; + kappaEVec[g] = kappaFVec[g]; + } + } + } else { + if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + if constexpr (gamma_ != 1.0) { + kappaPVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + kappaEVec = kappaPVec; + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, + delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, + kappa_expo_and_lower_value[0]); + } else { + kappaFVec = kappaPVec; + } + } else { + kappaFVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + // Note that alpha_F has not been changed in the Newton iteration + if constexpr (gamma_ != 1.0) { + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, + delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, + kappa_expo_and_lower_value[0]); + } else { // fall back to bin-center opacity + kappaFVec = ComputeBinCenterOpacity(radBoundaries_g_copy, kappa_expo_and_lower_value); + } + } else { + kappaFVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + } + } } - // Note that alpha_F has not been changed in the Newton iteration - kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); } - AMREX_ASSERT(!kappaFVec.hasnan()); dMomentum = {0., 0., 0.}; @@ -1409,44 +1680,47 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne std::array gasVel{}; std::array v_terms{}; - auto Fx = Frad_t0[0]; - auto Fy = Frad_t0[1]; - auto Fz = Frad_t0[2]; - auto fx = Fx / (c_light_ * erad); - auto fy = Fy / (c_light_ * erad); - auto fz = Fz / (c_light_ * erad); + auto fx = Frad_t0[0] / (c_light_ * erad); + auto fy = Frad_t0[1] / (c_light_ * erad); + auto fz = Frad_t0[2] / (c_light_ * erad); double F_coeff = chat * rho * kappaFVec[g] * dt * lorentz_factor; auto Tedd = ComputeEddingtonTensor(fx, fy, fz); for (int n = 0; n < 3; ++n) { // compute thermal radiation term - double v_term = NAN; - - if constexpr (opacity_model_ == OpacityModel::user) { - v_term = kappaPVec[g] * fourPiBoverC[g] * lorentz_factor_v; + double Planck_term = NAN; + if constexpr (opacity_model_ == OpacityModel::single_group) { + Planck_term = kappaPVec[g] * fourPiBoverC[g] * lorentz_factor_v; // compute (kappa_F - kappa_E) term if (kappaFVec[g] != kappaEVec[g]) { - v_term += (kappaFVec[g] - kappaEVec[g]) * erad * std::pow(lorentz_factor_v, 3); + Planck_term += (kappaFVec[g] - kappaEVec[g]) * erad * std::pow(lorentz_factor_v, 3); + } + } else { + if constexpr (include_delta_B) { + Planck_term = kappaPVec[g] * fourPiBoverC[g] - 1.0 / 3.0 * delta_nu_kappa_B_at_edge[g]; + } else { + Planck_term = kappaPVec[g] * fourPiBoverC[g]; } - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { - v_term = kappaPVec[g] * fourPiBoverC[g] * (2.0 - kappa_expo_and_lower_value[0][g] - alpha_B[g]) / 3.0; } - - v_term *= chat * dt * gasMtm0[n]; + Planck_term *= chat * dt * gasMtm0[n]; // compute radiation pressure double pressure_term = 0.0; for (int z = 0; z < 3; ++z) { pressure_term += gasMtm0[z] * Tedd[n][z] * erad; } - if constexpr (opacity_model_ == OpacityModel::user) { + if constexpr (opacity_model_ == OpacityModel::single_group) { pressure_term *= chat * dt * kappaFVec[g] * lorentz_factor_v; - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { - pressure_term *= chat * dt * kappaEVec[g] * (kappa_expo_and_lower_value[0][g] + 1.0); + } else { + // Simplification: assuming Eddington tensors are the same for all groups, we have kappaP = kappaE + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + pressure_term *= chat * dt * kappaEVec[g]; + } else { + pressure_term *= chat * dt * (1.0 + kappa_expo_and_lower_value[0][g]) * kappaEVec[g]; + } } - v_term += pressure_term; - v_terms[n] = v_term; + v_terms[n] = Planck_term + pressure_term; } if constexpr (beta_order_ == 1) { @@ -1571,30 +1845,25 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // copy work to work_prev work_prev[g] = work[g]; // compute new work term from the updated radiation flux and velocity - if constexpr (opacity_model_ == OpacityModel::user) { + // work = v * F * chi + if constexpr (opacity_model_ == OpacityModel::single_group) { work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * chat / (c * c) * lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; - } else if constexpr (opacity_model_ == OpacityModel::piecewisePowerLaw) { - for (int n = 0; n < 3; ++n) { - work[n] = 0.0; - } - for (int n = 0; n < 3; ++n) { - alpha_F = ComputeRadQuantityExponents(Frad_t1[n], radBoundaries_g_copy); - kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_F); - for (int g = 0; g < nGroups_; ++g) { - work[g] += (kappa_expo_and_lower_value[0][g] + 1.0) * gasMtm0[n] * kappaFVec[g] * Frad_t1[n][g]; - } - } - for (int g = 0; g < nGroups_; ++g) { - work[g] *= chat * dt / (c * c); - } + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * kappaFVec[g] * chat / + (c * c) * dt; + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; } } // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop const double lag_tol = 1.0e-13; if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) / Etot0 < lag_tol) || - (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec))) { + (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || + (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { break; } } // end full-step iteration diff --git a/tests/MarshakVaytet.in b/tests/MarshakVaytet.in new file mode 100644 index 000000000..5e2659372 --- /dev/null +++ b/tests/MarshakVaytet.in @@ -0,0 +1,24 @@ +cfl = 0.4 + +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 20.0 1.0 1.0 +geometry.is_periodic = 0 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 64 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 1 \ No newline at end of file diff --git a/tests/RadhydroBB.in b/tests/RadhydroBB.in new file mode 100644 index 000000000..eafd1c9c8 --- /dev/null +++ b/tests/RadhydroBB.in @@ -0,0 +1,23 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 64.0 1.0 1.0 +geometry.is_periodic = 1 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 8 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this +amr.max_grid_size = 64 + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 1 \ No newline at end of file diff --git a/tests/radshockMG.in b/tests/radshockMG.in index 63ec0aeca..79aa2fde8 100644 --- a/tests/radshockMG.in +++ b/tests/radshockMG.in @@ -13,7 +13,7 @@ amr.v = 0 # verbosity in Amr # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 128 4 4 +amr.n_cell = 64 4 4 amr.max_level = 0 # number of levels = max_level + 1 amr.blocking_factor = 4 # grid size must be divisible by this amr.max_grid_size_x = 128