diff --git a/Util/exact_riemann/ci-benchmarks/test1-riemann.out b/Util/exact_riemann/ci-benchmarks/test1-riemann.out index e9595d449f..4c6ec142cf 100644 --- a/Util/exact_riemann/ci-benchmarks/test1-riemann.out +++ b/Util/exact_riemann/ci-benchmarks/test1-riemann.out @@ -27,92 +27,92 @@ 26 199218.75 10000000 0 8.39953973425e+23 100000000 1.57524902122e+17 1.45276631226 27 207031.25 10000000 0 8.39953973425e+23 100000000 1.57524902122e+17 1.45276631226 28 214843.75 10000000 0 8.39953973425e+23 100000000 1.57524902122e+17 1.45276631226 - 29 222656.25 9937642.65309 2183663.02514 8.32354752664e+23 99662325.1703 1.57000231185e+17 1.45306476475 - 30 230468.75 9709901.70513 10251910.7507 8.04771591954e+23 98420371.8931 1.55068583081e+17 1.45417444944 - 31 238281.25 9486385.89679 18317060.0138 7.77963804416e+23 97188380.487 1.53148584392e+17 1.45529422204 - 32 246093.75 9266966.51457 26381560.7449 7.51905422388e+23 95965424.1262 1.51239653975e+17 1.45642406962 - 33 253906.25 9051793.54434 34437540.8547 7.26603668573e+23 94752214.9833 1.49343647326e+17 1.45756305571 - 34 261718.75 8840609.37756 42492257.1495 7.02017795716e+23 93548005.7883 1.47458863753e+17 1.45871305459 - 35 269531.25 8633401.34303 50543942.3276 6.78136325302e+23 92353089.2744 1.45585728968e+17 1.45987397189 - 36 277343.75 8430076.379 58593934.6682 6.54938678046e+23 91167662.8803 1.43723950117e+17 1.46104536843 - 37 285156.25 8230715.49759 66636610.8903 6.32424252216e+23 89991950.432 1.4187481839e+17 1.46222608033 - 38 292968.75 8034832.96548 74689465.5594 6.10529026499e+23 88823050.309 1.40034330633e+17 1.4634189711 - 39 300781.25 7843276.70737 82714958.8575 5.89337968843e+23 87667204.6781 1.38211065865e+17 1.46461939352 - 40 308593.75 7655129.37331 90748789.6847 5.68739786019e+23 86518471.2351 1.36396884033e+17 1.46583253085 - 41 316406.25 7470624.82291 98778876.9484 5.48751535566e+23 85379441.8221 1.34594557178e+17 1.46705592313 - 42 324218.75 7289705.25339 106805289.986 5.29358018796e+23 84249099.2209 1.32804057278e+17 1.46828932618 - 43 332031.25 7112199.44582 114833307.874 5.10532342594e+23 83126697.3171 1.31024222179e+17 1.46953405957 - 44 339843.75 6938388.96123 122847561.466 4.92295668399e+23 82015768.6131 1.29258477347e+17 1.47078862532 - 45 347656.25 6767903.17226 130862507.974 4.74600318924e+23 80912514.0778 1.27503605304e+17 1.47205472605 - 46 355468.75 6600803.96836 138872804.144 4.57444660615e+23 79818761.6342 1.25760804705e+17 1.47333091693 - 47 363281.25 6436864.48488 146886990.391 4.40797484146e+23 78732196.5312 1.24028216561e+17 1.47461852775 - 48 371093.75 6276476.38528 154883169.503 4.24690379729e+23 77656178.4583 1.22310573218e+17 1.47591487426 - 49 378906.25 6118495.13624 162916603.313 4.09001480639e+23 76583886.5211 1.20596074142e+17 1.477229381 - 50 386718.75 5965013.66448 170877440.433 3.93930119858e+23 75531305.6649 1.18908133569e+17 1.47854322379 - 51 394531.25 5813726.09716 178882079.158 3.79241466765e+23 74479511.5774 1.17221977469e+17 1.47987508489 - 52 402343.75 5665975.12416 186857492.862 3.65059344482e+23 73439553.2492 1.15553059891e+17 1.48121323952 - 53 410156.25 5521067.30265 194837976.618 3.513094079e+23 72406739.4884 1.13894173352e+17 1.48256442274 - 54 417968.75 5378402.95374 202855367.503 3.37929158716e+23 71377997.8179 1.12238824303e+17 1.48393366593 - 55 425781.25 5240033.70371 210790514.204 3.25103057011e+23 70370282.3137 1.10611565975e+17 1.48529979966 - 56 433593.75 5103573.55859 218776898.035 3.12602167336e+23 69361237.9415 1.08984897717e+17 1.48668564763 - 57 441406.25 4970427.23891 226730339.082 3.00549389365e+23 68364621.0996 1.07376071864e+17 1.48807764382 - 58 449218.75 4839918.84486 234687829.979 2.88876341941e+23 67374864.5081 1.05777552776e+17 1.48948243693 - 59 457031.25 4711261.74129 242696475.442 2.77507944248e+23 66387937.5771 1.04180042678e+17 1.49090755103 - 60 464843.75 4586905.30072 250599537.404 2.6665322397e+23 65424617.9561 1.02614733863e+17 1.49232462878 - 61 472656.25 4464458.411 258544516.094 2.56095593839e+23 64461080.2071 1.01052210978e+17 1.49376061964 - 62 480468.75 4344587.32012 266486808.786 2.45887720238e+23 63504668.1356 9.9501378356e+16 1.49520837014 - 63 488281.25 4227166.6049 274432114.632 2.36013209786e+23 62555780.9675 9.79611624117e+16 1.4966680197 - 64 496093.75 4112462.57556 282359199.256 2.26488729767e+23 61622622.8145 9.64357526277e+16 1.49813506882 - 65 503906.25 4000144.68079 290287932.496 2.17280446758e+23 60691094.8085 9.49211318404e+16 1.49961347261 - 66 511718.75 3890255.59174 298212431.422 2.08386680121e+23 59768976.3974 9.34185476389e+16 1.50110328489 - 67 519531.25 3782705.88447 306136345.798 1.99794677181e+23 58853482.3226 9.19272837241e+16 1.50260437684 - 68 527343.75 3677089.80954 314087675.696 1.91467447198e+23 57942974.054 9.04421710374e+16 1.50412131456 - 69 535156.25 3574803.70947 321956681.817 1.83509085676e+23 57055441.9326 8.8983691651e+16 1.50563336672 - 70 542968.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 71 550781.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 72 558593.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 73 566406.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 74 574218.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 75 582031.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 76 589843.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 77 597656.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 78 605468.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 79 613281.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 80 621093.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 81 628906.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 82 636718.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 83 644531.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 84 652343.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 85 660156.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 86 667968.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 87 675781.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 88 683593.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 89 691406.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 90 699218.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 91 707031.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 92 714843.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 93 722656.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 94 730468.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 95 738281.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 96 746093.75 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 97 753906.25 3566780.04365 322581308.247 1.82889078257e+23 56979872.1834 8.88683143667e+16 1.50575384017 - 98 761718.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 99 769531.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 100 777343.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 101 785156.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 102 792968.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 103 800781.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 104 808593.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 105 816406.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 106 824218.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 107 832031.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 108 839843.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 109 847656.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 110 855468.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 111 863281.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 112 871093.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 113 878906.25 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 - 114 886718.75 2952883.40476 322581308.247 1.82889078257e+23 758657904.323 1.0975004821e+17 1.49554982554 + 29 222656.25 9928583.54116 2501799.6711 8.31252420485e+23 99613180.0658 1.56923858619e+17 1.45310831422 + 30 230468.75 9704593.18704 10441751.4975 8.04131860434e+23 98391261.8227 1.55023263939e+17 1.45420068832 + 31 238281.25 9486357.88507 18318080.0803 7.7796046161e+23 97188228.2969 1.53148342292e+17 1.45529436435 + 32 246093.75 9230622.7828 27731832.085 7.47614190455e+23 95761494.9436 1.50921101437e+17 1.45661424075 + 33 253906.25 9042199.64069 34800232.9967 7.25481418396e+23 94697746.0256 1.49258542366e+17 1.45761458884 + 34 261718.75 8837352.22682 42617669.1639 7.01640532398e+23 93529275.3566 1.47429603369e+17 1.4587310511 + 35 269531.25 8633402.37919 50543902.8796 6.78136437374e+23 92353017.6409 1.45585736953e+17 1.45987396536 + 36 277343.75 8424166.80801 58830179.3826 6.54267980773e+23 91132958.3428 1.43669477048e+17 1.46107990247 + 37 285156.25 8230629.92822 66640095.7209 6.32414638323e+23 89991443.7522 1.41874019579e+17 1.46222659419 + 38 292968.75 8034949.87676 74684613.1174 6.10542029385e+23 88823783.0186 1.40035436837e+17 1.4634182492 + 39 300781.25 7843096.36042 82722591.5607 5.8931809781e+23 87665807.5108 1.38209332502e+17 1.46462053761 + 40 308593.75 7649307.12199 90999860.6892 5.68105804809e+23 86482508.2736 1.3634036201e+17 1.46587061839 + 41 316406.25 7470587.5927 98780517.0137 5.4874750423e+23 85378953.7541 1.34594186351e+17 1.46705617138 + 42 324218.75 7288269.80728 106869595.863 5.29204964271e+23 84239948.9418 1.32789754767e+17 1.46829925003 + 43 332031.25 7112320.19366 114827793.24 5.10545083595e+23 83127515.6114 1.31025441806e+17 1.46953320094 + 44 339843.75 6934566.20002 123025594.781 4.91896733388e+23 81990387.1944 1.29219365222e+17 1.47081661934 + 45 347656.25 6750486.78572 131690155.024 4.72803428499e+23 80798483.3299 1.27323012723e+17 1.47218608441 + 46 355468.75 6600412.29664 138891775.957 4.57404624767e+23 79815503.9137 1.25756679944e+17 1.4733339451 + 47 363281.25 6437003.88511 146880115.475 4.40811534033e+23 78732708.7965 1.24029691873e+17 1.474617414 + 48 371093.75 6276474.80815 154883249.082 4.24690223072e+23 77656181.1487 1.22310556462e+17 1.4759148873 + 49 378906.25 6112565.87908 163221240.401 4.08416108169e+23 76543467.7092 1.20531279715e+17 1.47727945729 + 50 386718.75 5948385.03461 171749502.523 3.92307380155e+23 75415166.7442 1.18723878247e+17 1.47868778629 + 51 394531.25 5813983.33262 178868346.808 3.79266247747e+23 74480433.3208 1.17224847464e+17 1.47987278158 + 52 402343.75 5632712.05483 188675252.44 3.61889014912e+23 73203050.8303 1.15174221081e+17 1.48151991526 + 53 410156.25 5521066.68926 194838013.088 3.51309344074e+23 72406628.3682 1.13894164306e+17 1.48256442775 + 54 417968.75 5378917.49304 202826156.038 3.37977138172e+23 71381811.9689 1.1224483654e+17 1.48392865739 + 55 425781.25 5236495.59633 210995568.911 3.24776964761e+23 70342394.3125 1.10569632197e+17 1.48533522384 + 56 433593.75 5103854.20568 218760327.567 3.12627654366e+23 69361980.2004 1.08988241021e+17 1.48668274662 + 57 441406.25 4956749.23672 227556710.807 2.9931935141e+23 68260486.4015 1.0720953507e+17 1.48822297337 + 58 449218.75 4839745.89329 234698487.502 2.88860961111e+23 67373438.3099 1.0577541791e+17 1.48948432483 + 59 457031.25 4690792.87241 243986151.57 2.75712171442e+23 66229884.3651 1.03923846886e+17 1.49113806553 + 60 464843.75 4584236.44142 250770998.196 2.66421598655e+23 65401136.6284 1.02580855826e+17 1.49235546707 + 61 472656.25 4464443.70532 258545525.124 2.56094219284e+23 64458438.244 1.01051975264e+17 1.49376077765 + 62 480468.75 4344308.5769 266505498.133 2.45864071824e+23 63501059.8528 9.94977217217e+16 1.495211777 + 63 488281.25 4198311.4832 276410556.084 2.33605802435e+23 62321207.5161 9.75793895051e+16 1.49703316141 + 64 496093.75 4111606.10228 282419102.001 2.26417890286e+23 61611339.532 9.642420376e+16 1.49814615416 + 65 503906.25 4000041.78464 290295327.218 2.17271946838e+23 60687326.4757 9.49196806953e+16 1.49961482884 + 66 511718.75 3890253.04293 298212667.589 2.08386366995e+23 59766222.6724 9.34184619422e+16 1.50110330311 + 67 519531.25 3782769.06389 306131662.364 1.99799647911e+23 58852893.6708 9.19281449819e+16 1.50260347548 + 68 527343.75 3677627.42215 314046760.692 1.9150955856e+23 57947741.148 9.04497857546e+16 1.5041134824 + 69 535156.25 3573422.36381 322064216.687 1.83402150589e+23 57038358.4286 8.89637619826e+16 1.50565406349 + 70 542968.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 71 550781.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 72 558593.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 73 566406.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 74 574218.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 75 582031.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 76 589843.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 77 597656.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 78 605468.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 79 613281.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 80 621093.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 81 628906.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 82 636718.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 83 644531.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 84 652343.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 85 660156.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 86 667968.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 87 675781.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 88 683593.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 89 691406.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 90 699218.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 91 707031.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 92 714843.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 93 722656.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 94 730468.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 95 738281.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 96 746093.75 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 97 753906.25 3566771.3761 322581533.648 1.82889264617e+23 57003324.3672 8.8868625651e+16 1.50575410561 + 98 761718.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 99 769531.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 100 777343.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 101 785156.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 102 792968.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 103 800781.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 104 808593.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 105 816406.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 106 824218.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 107 832031.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 108 839843.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 109 847656.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 110 855468.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 111 863281.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 112 871093.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 113 878906.25 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 + 114 886718.75 2952884.63461 322581533.648 1.82889264617e+23 758658820.738 1.09750124478e+17 1.49554976549 115 894531.25 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 116 902343.75 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 117 910156.25 1000000 0 2.55457320691e+22 999999.999997 4.08260803009e+16 1.57201404928 diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H new file mode 100644 index 0000000000..9592225fa0 --- /dev/null +++ b/Util/exact_riemann/riemann_rarefaction.H @@ -0,0 +1,432 @@ +#ifndef RIEMANN_RAREFACTION_H +#define RIEMANN_RAREFACTION_H + +#include + +#include + +#include +#include + + +AMREX_INLINE +void +riemann_invariant_rhs(const amrex::Real p, const amrex::Real tau, const amrex::Real u, + const amrex::Real* xn, const int iwave, + amrex::Real& T, + amrex::Real& dtaudp, amrex::Real& dudp) { + + amrex::ignore_unused(u); + + // here, p is out independent variable, and tau, u are the + // dependent variables. We return the derivatives of these + // wrt p for integration. + + // T should be an initial guess coming in + + // get the thermodynamics + + eos_rep_t eos_state; + eos_state.rho = 1.0_rt / tau; + eos_state.p = p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = T; + + eos(eos_input_rp, eos_state); + + T = eos_state.T; + + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); + + dtaudp = -1.0_rt / (C * C); + + if (iwave == 1) { + dudp = -1.0_rt / C; + } else if (iwave == 3) { + dudp = 1.0_rt / C; + } + +} + + +AMREX_INLINE +void +riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex::Real p, + const amrex::Real* xn, const int iwave, + amrex::Real& T, + amrex::Real& dtaudu, amrex::Real& dpdu) { + + // here, u is out independent variable, and tau, p are the + // dependent variables. We return the derivatives of these + // wrt u for integration. + + // here T is an initial guess + + // get the thermodynamics + + amrex::ignore_unused(u); + + std::cout << "in riemann_invariant_rhs2: " << u << " " << tau << " " << p << std::endl; + eos_rep_t eos_state; + eos_state.rho = 1.0_rt / tau; + eos_state.p = p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = T; + + eos(eos_input_rp, eos_state); + + T = eos_state.T; + + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); + + if (iwave == 3) { + dpdu = C; + dtaudu = -1.0_rt / C; + + } else if (iwave == 1) { + dpdu = -C; + dtaudu = 1.0_rt / C; + } +} + + +AMREX_INLINE +std::pair +single_step_p(const amrex::Real pstar0, const amrex::Real dp, + const amrex::Real u0, const amrex::Real tau0, + const amrex::Real* xn, const int iwave, amrex::Real& T) { + + // this takes a single step of the Riemann invariant system where + // p is the independent variable + + amrex::Real dtaudp1, dudp1; + riemann_invariant_rhs(pstar0, tau0, u0, xn, iwave, T, dtaudp1, dudp1); + + amrex::Real dtaudp2, dudp2; + riemann_invariant_rhs(pstar0+0.5*dp, tau0+0.5*dp*dtaudp1, u0+0.5*dp*dudp1, + xn, iwave, T, dtaudp2, dudp2); + + amrex::Real dtaudp3, dudp3; + riemann_invariant_rhs(pstar0+0.5*dp, tau0+0.5*dp*dtaudp2, u0+0.5*dp*dudp2, + xn, iwave, T, dtaudp3, dudp3); + + amrex::Real dtaudp4, dudp4; + riemann_invariant_rhs(pstar0+dp, tau0+dp*dtaudp3, u0+dp*dudp3, + xn, iwave, T, dtaudp4, dudp4); + + amrex::Real u = u0 + (1.0_rt/6.0_rt) * dp * (dudp1 + 2.0_rt * dudp2 + 2.0_rt * dudp3 + dudp4); + amrex::Real tau = tau0 + (1.0_rt/6.0_rt) * dp * (dtaudp1 + 2.0_rt * dtaudp2 + 2.0_rt * dtaudp3 + dtaudp4); + + return {u, tau}; +} + + +AMREX_INLINE +std::pair +single_step_u(const amrex::Real u0, const amrex::Real du, + const amrex::Real tau0, const amrex::Real p0, + const amrex::Real* xn, const int iwave, amrex::Real& T) { + + // this takes a single step of the Riemann invariant system where + // u is the independent variable + + amrex::Real dtaudu1, dpdu1; + riemann_invariant_rhs2(u0, tau0, p0, xn, iwave, T, dtaudu1, dpdu1); + + amrex::Real dtaudu2, dpdu2; + riemann_invariant_rhs2(u0+0.5*du, tau0+0.5*du*dtaudu1, p0+0.5*du*dpdu1, + xn, iwave, T, dtaudu2, dpdu2); + + amrex::Real dtaudu3, dpdu3; + riemann_invariant_rhs2(u0+0.5*du, tau0+0.5*du*dtaudu2, p0+0.5*du*dpdu2, + xn, iwave, T, dtaudu3, dpdu3); + + amrex::Real dtaudu4, dpdu4; + riemann_invariant_rhs2(u0+du, tau0+du*dtaudu3, p0+du*dpdu3, + xn, iwave, T, dtaudu4, dpdu4); + + amrex::Real p = p0 + (1.0_rt/6.0_rt) * du * (dpdu1 + 2.0_rt * dpdu2 + 2.0_rt * dpdu3 + dpdu4); + amrex::Real tau = tau0 + (1.0_rt/6.0_rt) * du * (dtaudu1 + 2.0_rt * dtaudu2 + 2.0_rt * dtaudu3 + dtaudu4); + + return {tau, p}; + +} + +AMREX_INLINE +void +rarefaction(const amrex::Real pstar, + const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, + const amrex::Real* xn, const int iwave, + amrex::Real& Z_s, amrex::Real& W_s, amrex::Real& rhostar) { + + // Compute Z_s = C(p*, rho*) for a rarefaction connecting the + // state to the star region by integrating the Riemann invariant + // from p_s to pstar. + // + // The system is: + // + // dtau/dp = -1/C**2 + // du/dp = +/- 1/C (+ for 1-wave, - for 3-wave) + + const double S1{0.9}; + const double S2{4.0}; + + constexpr amrex::Real tol = 1.e-5; + + // initial conditions + + amrex::Real tau = 1.0_rt / rho_s; + amrex::Real u = u_s; + amrex::Real p = p_s; + + // initial guess at integration step + amrex::Real dp = (pstar - p_s) / static_cast(100); + + amrex::Real T = problem::initial_temp_guess; + + // adaptive RK4 loop + + // note: we should only be in here if p > pstar, since otherwise we'd have a shock + + amrex::Real dp_new{dp}; + int nstep{0}; + + while (p > pstar) { + + amrex::Real u0{u}; + amrex::Real tau0{tau}; + + amrex::Real u_new; + amrex::Real tau_new; + + amrex::Real rel_error = std::numeric_limits::max(); + + // take a step + + while (rel_error > tol) { + dp = dp_new; + if (p + dp < pstar) { + dp = pstar - p; + } + + // take 2 half steps + amrex::Real u_tmp; + amrex::Real tau_tmp; + std::tie(u_tmp, tau_tmp) = single_step_p(p, 0.5*dp, u0, tau0, xn, iwave, T); + std::tie(u_new, tau_new) = single_step_p(p + 0.5*dp, 0.5*dp, u_tmp, tau_tmp, xn, iwave, T); + + // now take a single step to cover dp + auto [u_single, tau_single] = single_step_p(p, dp, u0, tau0, xn, iwave, T); + + // estimate the relative error + amrex::Real u_err = std::abs(u_new - u_single); + if (u0 != 0.0) { + u_err /= u0; + } + amrex::Real tau_err = std::abs((tau_new - tau_single) / tau0); + + rel_error = std::max(u_err, tau_err); + + // adaptive step estimate + + double dp_est = S1 * dp * std::pow(std::abs(tol/rel_error), 0.2); + dp_new = dp_est; //std::clamp(S1*dp_est, dp/S2, S2*dp); + + } + + // success + p += dp; + u = u_new; + tau = tau_new; + + nstep++; + } + + // Z_s is just the Lagrangian sound speed + + eos_rep_t eos_state; + eos_state.rho = 1.0_rt / tau; + eos_state.p = p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = problem::initial_temp_guess; + + eos(eos_input_rp, eos_state); + + Z_s = std::sqrt(eos_state.gam1 * p / tau); + + // also need W_s -- this is C&G Eq. 16. u above is ustar_s. + + if (u == u_s) { + W_s = Z_s; + } else { + W_s = std::abs(pstar - p_s) / std::abs(u - u_s); + } + + rhostar = 1.0_rt / tau; +} + + +AMREX_INLINE +void +rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, + const amrex::Real* xn, const int iwave, const amrex::Real xi, + amrex::Real& rho, amrex::Real& p, amrex::Real& u) { + + const int npts = 1000; + + // here we integrate the Riemann invariants for a rarefaction up to + // some intermediate u (between u_s and ustar). This accounts for + // the fact that we are inside the rarefaction. + // + // We reformulate the system of ODEs from C&G Eq. 13 to make u the + // dependent variable. Now we solve: + // + // dp/du = C; dtau/du = -1/C for the 1-wave + // dp/du = -C; dtau/du = 1/C for the 3-wave + // + // we actually don't know the stopping point. For the 1-wave, we + // stop at u = xi + c, for the 3-wave, we stop at u = xi - c, where + // c is computed as we step. + + const double S1{0.9}; + + constexpr amrex::Real tol = 1.e-8; + + // initial conditions + + amrex::Real tau = 1.0_rt / rho_s; + u = u_s; + p = p_s; + + // compute c and estimate the velocity for the endpoint of integration + + eos_rep_t eos_state; + eos_state.rho = 1.0_rt / tau; + eos_state.p = p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = problem::initial_temp_guess; + + eos(eos_input_rp, eos_state); + + amrex::Real c = std::sqrt(eos_state.gam1 * p * tau); + + amrex::Real ustop = (iwave == 1) ? xi + c : xi - c; + + // initial guess at the step size -- we'll adapt below + + amrex::Real du = (ustop - u_s) / static_cast(10); + + bool finished = false; + + std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl; + + // this will be used as an initial guess to accelerate the EOS inversions + + amrex::Real T = eos_state.T; + + amrex::Real du_new{du}; + + while (! finished) { + + amrex::Real p0{p}; + amrex::Real tau0{tau}; + + amrex::Real p_new; + amrex::Real tau_new; + + amrex::Real rel_error = std::numeric_limits::max(); + + // take a step + + while (rel_error > tol) { + du = du_new; + + // take 2 half steps + amrex::Real p_tmp; + amrex::Real tau_tmp; + std::tie(tau_tmp, p_tmp) = single_step_u(u, 0.5*du, tau0, p0, xn, iwave, T); + std::tie(tau_new, p_new) = single_step_u(u + 0.5*du, 0.5*du, tau_tmp, p_tmp, xn, iwave, T); + + // now take a single step to cover du + auto [tau_single, p_single] = single_step_u(u, du, tau0, p0, xn, iwave, T); + + // estimate the relative error + amrex::Real p_err = std::abs((p_new - p_single) / p0); + amrex::Real tau_err = std::abs((tau_new - tau_single) / tau0); + + rel_error = std::max(p_err, tau_err); + + // adaptive step estimate + + // if our original step was too small, the error might be zero, in which case, + // just bump up the step for next time + + amrex::Real du_est{}; + if (rel_error > 0) { + du_est = S1 * du * std::pow(std::abs(tol / rel_error), 0.2); + } else { + du_est = 10.0 * du; + } + du_new = du_est; + + } + + // success + u += du; + p = p_new; + tau = tau_new; + + // compute c and re-estimate the endpoint of integration + + eos_rep_t eos_state; + eos_state.rho = 1.0_rt/tau; + eos_state.p = p; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = problem::initial_temp_guess; + + eos(eos_input_rp, eos_state); + + c = std::sqrt(eos_state.gam1 * p * tau); + + ustop = (iwave == 1) ? xi + c : xi - c; + + // check the step size for the next step to ensure we don't + // go past the stopping velocity, ustop + + if (du_new * u > 0.0_rt) { + while (std::abs(u + du_new) > std::abs(ustop) && du_new != 0.0_rt) { + du_new *= 0.5_rt; + } + } else { + if (u > 0.0_rt) { + while (u + du_new < ustop && du_new != 0.0_rt) { + du_new *= 0.5_rt; + } + + } else { + while (u + du_new > ustop && du_new != 0.0_rt) { + du_new *= 0.5_rt; + } + } + } + + if (std::abs(du_new) < riemann_solver::tol * std::abs(u) + riemann_solver::tol) { + finished = true; + } + + } + + rho = 1.0_rt / tau; +} + +#endif diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index c2580dd258..9db07b34df 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -4,6 +4,7 @@ #include #include +#include using namespace amrex::literals; diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H new file mode 100644 index 0000000000..93fc4c17ac --- /dev/null +++ b/Util/exact_riemann/riemann_shock.H @@ -0,0 +1,204 @@ +#ifndef RIEMANN_SHOCK_H +#define RIEMANN_SHOCK_H + +#include + +#include +#include + + +AMREX_INLINE +void +W_s_shock(const amrex::Real W_s, const amrex::Real pstar, + const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, + const amrex::Real* xn, + amrex::Real& rhostar_s, eos_rep_t& eos_state, amrex::Real& f, amrex::Real& fprime) { + + // we need rhostar -- get it from the R-H conditions + + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + rhostar_s = 1.0_rt / taustar_s; + + // get the thermodynamics + + eos_state.rho = rhostar_s; + eos_state.p = pstar; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = problem::initial_temp_guess; + + eos(eos_input_rp, eos_state); + + // compute the correction + + f = W_s * W_s * (eos_state.e - e_s) - 0.5_rt * (pstar * pstar - p_s * p_s); + + // we need de/drho at constant p -- this is not returned by the EOS + + amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; + + fprime = 2.0_rt * W_s * (eos_state.e - e_s) - + 2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s; + +} + + +AMREX_INLINE +void +newton_shock(amrex::Real& W_s, const amrex::Real pstar, + const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, + const amrex::Real* xn, + amrex::Real* rhostar_hist, amrex::Real* Ws_hist, + eos_rep_t& eos_state, bool& converged) { + + // Newton iterations -- we are zeroing the energy R-H jump condition + // W^2 [e] = 1/2 [p^2] + // + // we write f(W) = W^2 (e(pstar, rhostar_s) - e_s) - (1/2)(pstar^2 - p_s) + // + // and then compute f' + + converged = false; + + int iter = 1; + while (! converged && iter < riemann_solver::max_iters) { + + amrex::Real rhostar_s, f, fprime; + + W_s_shock(W_s, pstar, rho_s, p_s, e_s, xn, + rhostar_s, eos_state, f, fprime); + + amrex::Real dW = -f / fprime; + + if (std::abs(dW) < riemann_solver::tol * W_s) { + converged = true; + } + + W_s = amrex::min(2.0_rt * W_s, amrex::max(0.5_rt * W_s, W_s + dW)); + + // store some history + + rhostar_hist[iter] = rhostar_s; + Ws_hist[iter] = W_s; + + iter++; + } +} + +AMREX_INLINE +void +shock(const amrex::Real pstar, + const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, + const amrex::Real* xn, + const amrex::Real gammaE_bar, const amrex::Real gammaC_bar, + amrex::Real& Z_s, amrex::Real& W_s) { + + amrex::ignore_unused(u_s); + + amrex::Real rhostar_hist[riemann_solver::max_iters], Ws_hist[riemann_solver::max_iters]; + + // compute the Z_s function for a shock following C&G Eq. 20 and + // 23. Here the "_s" variables are the state (L or R) that we are + // connecting to the star region through a shock. + + // first we need to compute W_s -- this is iterative because of + // the nonlinear EOS. We use the R-H jump conditions + the EOS + + // get the s-state energy, e_s + + eos_rep_t eos_state; + eos_state.rho = rho_s; + eos_state.p = p_s; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = xn[n]; + } + eos_state.T = problem::initial_temp_guess; + + eos(eos_input_rp, eos_state); + + amrex::Real e_s = eos_state.e; + + // to kick things off, we need a guess for W_s. We'll use the + // approximation from Colella & Glaz (Eq. 34), which in turn + // makes an approximation for gammaE_star, using equation 31. + + amrex::Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; + + amrex::Real gammaE_star = gammaE_s + + 2.0_rt * (1.0_rt - gammaE_bar / gammaC_bar) * (gammaE_bar - 1.0_rt) * + (pstar - p_s) / (pstar + p_s); + + std::cout << "pstar, ps = " << pstar << " " << p_s << " " << gammaE_s << " " << gammaE_star << std::endl; + std::cout << (pstar/rho_s - (gammaE_star - 1.0_rt)/(gammaE_s - 1.0_rt) * p_s/rho_s); + std::cout << (pstar + 0.5_rt * (gammaE_star - 1.0_rt) * (pstar + p_s)); + + // there is a pathological case that if p_s - pstar ~ 0, the root finding + // just doesn't work. In this case, we use the ideas from CG, Eq. 35, and + // take W = sqrt(gamma p rho) + + if (pstar - p_s < riemann_solver::tol_p * p_s) { + W_s = std::sqrt(eos_state.gam1 * p_s * rho_s); + } else { + W_s = std::sqrt((pstar - p_s) * + (pstar + 0.5_rt * (gammaE_star - 1.0_rt) * (pstar + p_s)) / + (pstar / rho_s - (gammaE_star - 1.0_rt) / (gammaE_s - 1.0_rt) * p_s / rho_s)); + } + + // we need rhostar -- get it from the R-H conditions + + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + amrex::Real rhostar_s; + + if (taustar_s < 0.0_rt) { + rhostar_s = riemann_solver::smallrho; + W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s)); + } + + // newton + + bool converged; + newton_shock(W_s, pstar, rho_s, p_s, e_s, xn, + rhostar_hist, Ws_hist, + eos_state, converged); + + + // now did we converge? + + if (! converged) { + for (int i = 0; i < riemann_solver::max_iters; ++i) { + std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl; + } + + amrex::Error("shock did not converge"); + } + + + // now that we have W_s, we can get rhostar from the R-H conditions + // (C&G Eq. 12) + + taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + rhostar_s = 1.0_rt / taustar_s; + + // next we compute the derivative dW_s/dpstar -- the paper gives + // dW**2/dpstar (Eq. 23), so we take 1/2W of that + + amrex::Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); + + amrex::Real p_e = eos_state.dpdT / eos_state.dedT; + amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; + + amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho; + + amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / + ((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s)); + + amrex::Real dWdpstar = 0.5_rt * dW2dpstar / W_s; + + // finally compute Z_s + + Z_s = W_s * W_s / (W_s - dWdpstar * (pstar - p_s)); + +} + +#endif diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index b29c81dd50..47d1ee4781 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -4,6 +4,8 @@ #include #include +#include +#include using namespace amrex::literals; @@ -15,7 +17,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: const amrex::Real tol = 1.e-10_rt; const amrex::Real smallp = 1.e-8_rt; - const amrex::Real SMALL = 1.e-13_rt; + const amrex::Real SMALL = 1.e-8_rt; // get the initial sound speeds @@ -69,7 +71,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: if (W_l == W_r) { pstar = 0.5_rt * (p_l + p_r + W_l * (u_l - u_r)); } else { - pstar = ((W_r*p_l + W_l*p_r) + W_l*W_r*(u_l - u_r))/(W_l + W_r); + pstar = ((W_r * p_l + W_l * p_r) + W_l * W_r * (u_l - u_r)) / (W_l + W_r); } pstar = amrex::max(pstar, smallp); @@ -77,10 +79,33 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: // find the exact pstar and ustar - std::cout << "solving for star state: " << rho_l << " " << u_l << " " << p_l << " " << rho_r << " " << u_r << " " << p_r << std::endl; + std::cout << "solving for star state: " << std::endl; + std::cout << rho_l << " " << u_l << " " << p_l << std::endl; + std::cout << rho_r << " " << u_r << " " << p_r << std::endl; // this procedure follows directly from Colella & Glaz 1985, section 1 + // the basic idea is that we want to find the pstar that satisfies: + // + // ustar_l(pstar) - ustar_r(pstar) = 0 + // + // where ustar_l connects the left state to the star state across + // the left wave. this wave can be a shock or rarefaction, so we + // need to use either the Riemann invariants (rarefaction) or + // Rankine-Hugoniot jump conditions (shock). + // + // We use a Newton method. It takes the form: + // + // pstar <-- pstar - Z_l Z_r (ustar_r - ustar_l) / (Z_l + Z_r) + // + // where Z_s = | dpstar / dustar_s | is the derivative in the p-u plane + // + // and we get ustar via: + // + // ustar_s = u_s +/- (pstar - p_s) / W_s + // + // where W_s is the wave speed (this also depends on shock or + // rarefaction) amrex::Real ustar_l, ustar_r; bool converged = false; @@ -88,11 +113,12 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: int iter = 1; while (! converged && iter < problem::riemann_max_iter) { - // compute Z_l and Z_r -- the form of these depend on whether the - // wave is a shock or a rarefaction + // compute Z_l, W_l and Z_r, W_r -- the form of these depend + // on whether the wave is a shock or a rarefaction amrex::Real Z_l, Z_r; + std::cout << "iteration: " << iter << " " << pstar << " " << p_l << " " << p_r << std::endl; // left wave if (pstar - p_l > SMALL * p_l) { diff --git a/Util/exact_riemann/riemann_support.H b/Util/exact_riemann/riemann_support.H index 9a96c23dab..56fd6b8aaf 100644 --- a/Util/exact_riemann/riemann_support.H +++ b/Util/exact_riemann/riemann_support.H @@ -3,11 +3,6 @@ #include -#include -#include - -using namespace amrex::literals; - namespace riemann_solver { const int max_iters = 100; @@ -16,478 +11,4 @@ namespace riemann_solver { const amrex::Real tol_p = 1.e-6_rt; } - -AMREX_INLINE -void -W_s_shock(const amrex::Real W_s, const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, const amrex::Real* xn, - amrex::Real& rhostar_s, eos_rep_t& eos_state, amrex::Real& f, amrex::Real& fprime) { - - // we need rhostar -- get it from the R-H conditions - - amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); - rhostar_s = 1.0_rt / taustar_s; - - // get the thermodynamics - - eos_state.rho = rhostar_s; - eos_state.p = pstar; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = problem::initial_temp_guess; - - eos(eos_input_rp, eos_state); - - // compute the correction - - f = W_s * W_s * (eos_state.e - e_s) - 0.5_rt * (pstar * pstar - p_s * p_s); - - // we need de/drho at constant p -- this is not returned by the EOS - - amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; - - fprime = 2.0_rt * W_s * (eos_state.e - e_s) - 2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s; - -} - - -AMREX_INLINE -void -newton_shock(amrex::Real& W_s, const amrex::Real pstar, - const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, - const amrex::Real* xn, - amrex::Real* rhostar_hist, amrex::Real* Ws_hist, - eos_rep_t& eos_state, bool& converged) { - - // Newton iterations -- we are zeroing the energy R-H jump condition - // W^2 [e] = 1/2 [p^2] - // - // we write f(W) = W^2 (e(pstar, rhostar_s) - e_s) - (1/2)(pstar^2 - p_s) - // - // and then compute f' - - converged = false; - - int iter = 1; - while (! converged && iter < riemann_solver::max_iters) { - - amrex::Real rhostar_s, f, fprime; - - W_s_shock(W_s, pstar, rho_s, p_s, e_s, xn, - rhostar_s, eos_state, f, fprime); - - amrex::Real dW = -f / fprime; - - if (std::abs(dW) < riemann_solver::tol * W_s) { - converged = true; - } - - W_s = amrex::min(2.0_rt * W_s, amrex::max(0.5_rt * W_s, W_s + dW)); - - // store some history - - rhostar_hist[iter] = rhostar_s; - Ws_hist[iter] = W_s; - - iter++; - } -} - -AMREX_INLINE -void -shock(const amrex::Real pstar, - const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, - const amrex::Real* xn, - const amrex::Real gammaE_bar, const amrex::Real gammaC_bar, - amrex::Real& Z_s, amrex::Real& W_s) { - - amrex::ignore_unused(u_s); - - amrex::Real rhostar_hist[riemann_solver::max_iters], Ws_hist[riemann_solver::max_iters]; - - // compute the Z_s function for a shock following C&G Eq. 20 and - // 23. Here the "_s" variables are the state (L or R) that we are - // connecting to the star region through a shock. - - // first we need to compute W_s -- this is iterative because of - // the nonlinear EOS. We use the R-H jump conditions + the EOS - - // get the s-state energy, e_s - - eos_rep_t eos_state; - eos_state.rho = rho_s; - eos_state.p = p_s; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = problem::initial_temp_guess; - - eos(eos_input_rp, eos_state); - - amrex::Real e_s = eos_state.e; - - // to kick things off, we need a guess for W_s. We'll use the - // approximation from Colella & Glaz (Eq. 34), which in turn - // makes an approximation for gammaE_star, using equation 31. - - amrex::Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; - - amrex::Real gammaE_star = gammaE_s + - 2.0_rt * (1.0_rt - gammaE_bar / gammaC_bar) * (gammaE_bar - 1.0_rt) * - (pstar - p_s) / (pstar + p_s); - - std::cout << "pstar, ps = " << pstar << " " << p_s << " " << gammaE_s << " " << gammaE_star << std::endl; - std::cout << (pstar/rho_s - (gammaE_star - 1.0_rt)/(gammaE_s - 1.0_rt) * p_s/rho_s); - std::cout << (pstar + 0.5_rt * (gammaE_star - 1.0_rt) * (pstar + p_s)); - - // there is a pathological case that if p_s - pstar ~ 0, the root finding - // just doesn't work. In this case, we use the ideas from CG, Eq. 35, and - // take W = sqrt(gamma p rho) - - if (pstar - p_s < riemann_solver::tol_p * p_s) { - W_s = std::sqrt(eos_state.gam1 * p_s * rho_s); - } else { - W_s = std::sqrt((pstar - p_s) * - (pstar + 0.5_rt * (gammaE_star - 1.0_rt) * (pstar + p_s)) / - (pstar / rho_s - (gammaE_star - 1.0_rt) / (gammaE_s - 1.0_rt) * p_s / rho_s)); - } - - // we need rhostar -- get it from the R-H conditions - - amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); - amrex::Real rhostar_s; - - if (taustar_s < 0.0_rt) { - rhostar_s = riemann_solver::smallrho; - W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s)); - } - - // newton - - bool converged; - newton_shock(W_s, pstar, rho_s, p_s, e_s, xn, - rhostar_hist, Ws_hist, - eos_state, converged); - - - // now did we converge? - - if (! converged) { - for (int i = 0; i < riemann_solver::max_iters; ++i) { - std::cout << i << " " << rhostar_hist[i] << " " << Ws_hist[i] << std::endl; - } - - amrex::Error("shock did not converge"); - } - - - // now that we have W_s, we can get rhostar from the R-H conditions - // (C&G Eq. 12) - - taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); - rhostar_s = 1.0_rt / taustar_s; - - // next we compute the derivative dW_s/dpstar -- the paper gives - // dW**2/dpstar (Eq. 23), so we take 1/2W of that - - amrex::Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); - - amrex::Real p_e = eos_state.dpdT / eos_state.dedT; - amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; - - amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho; - - amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / - ((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s)); - - amrex::Real dWdpstar = 0.5_rt * dW2dpstar / W_s; - - // finally compute Z_s - - Z_s = W_s * W_s / (W_s - dWdpstar * (pstar - p_s)); - -} - - -AMREX_INLINE -void -riemann_invariant_rhs(const amrex::Real p, const amrex::Real tau, const amrex::Real u, - const amrex::Real* xn, const int iwave, - amrex::Real& T, - amrex::Real& dtaudp, amrex::Real& dudp) { - - amrex::ignore_unused(u); - - // here, p is out independent variable, and tau, u are the - // dependent variables. We return the derivatives of these - // wrt p for integration. - - // T should be an initial guess coming in - - // get the thermodynamics - - eos_rep_t eos_state; - eos_state.rho = 1.0_rt / tau; - eos_state.p = p; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = T; - - eos(eos_input_rp, eos_state); - - T = eos_state.T; - - amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); - - dtaudp = -1.0_rt / (C * C); - - if (iwave == 1) { - dudp = -1.0_rt / C; - } else if (iwave == 3) { - dudp = 1.0_rt / C; - } - -} - - -AMREX_INLINE -void -riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex::Real p, - const amrex::Real* xn, const int iwave, - amrex::Real& T, - amrex::Real& dtaudu, amrex::Real& dpdu) { - - // here, u is out independent variable, and tau, p are the - // dependent variables. We return the derivatives of these - // wrt u for integration. - - // here T is an initial guess - - // get the thermodynamics - - amrex::ignore_unused(u); - - eos_rep_t eos_state; - eos_state.rho = 1.0_rt / tau; - eos_state.p = p; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = T; - - eos(eos_input_rp, eos_state); - - T = eos_state.T; - - amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); - - if (iwave == 3) { - dpdu = C; - dtaudu = -1.0_rt / C; - - } else if (iwave == 1) { - dpdu = -C; - dtaudu = 1.0_rt / C; - } -} - - -AMREX_INLINE -void -rarefaction(const amrex::Real pstar, - const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, - const amrex::Real* xn, const int iwave, - amrex::Real& Z_s, amrex::Real& W_s, amrex::Real& rhostar) { - - const int npts = 1000; - - // Compute Z_s = C for a rarefaction connecting the state to the star - // region by integrating the Riemann invariant from p_s to pstar. - // This means solving a system of ODEs. We use RK45. - - amrex::Real tau = 1.0_rt / rho_s; - amrex::Real u = u_s; - amrex::Real p = p_s; - - amrex::Real dp = (pstar - p_s) / static_cast(npts); - amrex::Real dp2 = 0.5_rt * dp; - - amrex::Real T = problem::initial_temp_guess; - - for (int i = 1; i <= npts; ++i) { - - // do 4th-order RT - - amrex::Real dtaudp1, dudp1; - riemann_invariant_rhs(p, tau, u, xn, iwave, T, dtaudp1, dudp1); - - amrex::Real dtaudp2, dudp2; - riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp1, u+dp2*dudp1, xn, iwave, T, dtaudp2, dudp2); - - amrex::Real dtaudp3, dudp3; - riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp2, u+dp2*dudp2, xn, iwave, T, dtaudp3, dudp3); - - amrex::Real dtaudp4, dudp4; - riemann_invariant_rhs(p+dp, tau+dp*dtaudp3, u+dp*dudp3, xn, iwave, T, dtaudp4, dudp4); - - p += dp; - u += (1.0_rt/6.0_rt) * dp * (dudp1 + 2.0_rt * dudp2 + 2.0_rt * dudp3 + dudp4); - tau += (1.0_rt/6.0_rt) * dp * (dtaudp1 + 2.0_rt * dtaudp2 + 2.0_rt * dtaudp3 + dtaudp4); - - } - - // Z_s is just the Lagrangian sound speed - - eos_rep_t eos_state; - eos_state.rho = 1.0_rt / tau; - eos_state.p = p; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = problem::initial_temp_guess; - - eos(eos_input_rp, eos_state); - - Z_s = std::sqrt(eos_state.gam1 * p / tau); - - // also need W_s -- this is C&G Eq. 16. u above is ustar_s. - - if (u == u_s) { - W_s = Z_s; - } else { - W_s = std::abs(pstar - p_s) / std::abs(u - u_s); - } - - rhostar = 1.0_rt / tau; -} - - -AMREX_INLINE -void -rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, - const amrex::Real* xn, const int iwave, const amrex::Real xi, - amrex::Real& rho, amrex::Real& p, amrex::Real& u) { - - const int npts = 1000; - - // here we integrate the Riemann invariants for a rarefaction up to - // some intermediate u (between u_s and ustar). This accounts for - // the fact that we are inside the rarefaction. - // - // We reformulate the system of ODEs from C&G Eq. 13 to make u the - // dependent variable. Now we solve: - // - // dp/du = C; dtau/du = -1/C for the 1-wave - // dp/du = -C; dtau/du = 1/C for the 3-wave - // - // we actually don't know the stopping point. For the 1-wave, we - // stop at u = xi + c, for the 3-wave, we stop at u = xi - c, where - // c is computed as we step. - - - amrex::Real tau = 1.0_rt / rho_s; - u = u_s; - p = p_s; - - // estimate - // compute c - - eos_rep_t eos_state; - eos_state.rho = 1.0_rt / tau; - eos_state.p = p; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = problem::initial_temp_guess; - - eos(eos_input_rp, eos_state); - - amrex::Real c = std::sqrt(eos_state.gam1 * p * tau); - - amrex::Real ustop; - if (iwave == 1) { - ustop = xi + c; - } else if (iwave == 3) { - ustop = xi - c; - } - - amrex::Real du = (ustop - u_s) / static_cast(npts); - - bool finished = false; - - std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl; - - amrex::Real du2 = 0.5_rt * du; - - amrex::Real T = eos_state.T; - - while (! finished) { - - // do 4th-order RT - - amrex::Real dtaudu1, dpdu1; - riemann_invariant_rhs2(u, tau, p, xn, iwave, T, dtaudu1, dpdu1); - - amrex::Real dtaudu2, dpdu2; - riemann_invariant_rhs2(u+du2, tau+du2*dtaudu1, p+du2*dpdu1, xn, iwave, T, dtaudu2, dpdu2); - - amrex::Real dtaudu3, dpdu3; - riemann_invariant_rhs2(u+du2, tau+du2*dtaudu2, p+du2*dpdu2, xn, iwave, T, dtaudu3, dpdu3); - - amrex::Real dtaudu4, dpdu4; - riemann_invariant_rhs2(u+du, tau+du*dtaudu3, p+du*dpdu3, xn, iwave, T, dtaudu4, dpdu4); - - u += du; - p += (1.0_rt/6.0_rt) * du * (dpdu1 + 2.0_rt * dpdu2 + 2.0_rt * dpdu3 + dpdu4); - tau += (1.0_rt/6.0_rt) * du * (dtaudu1 + 2.0_rt * dtaudu2 + 2.0_rt * dtaudu3 + dtaudu4); - - // compute c - - eos_rep_t eos_state; - eos_state.rho = 1.0_rt/tau; - eos_state.p = p; - for (int n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = xn[n]; - } - eos_state.T = problem::initial_temp_guess; - - eos(eos_input_rp, eos_state); - - c = std::sqrt(eos_state.gam1 * p * tau); - - // check the step size - - if (iwave == 1) { - ustop = xi + c; - } else if (iwave == 3) { - ustop = xi - c; - } - - if (du * u > 0.0_rt) { - while (std::abs(u + du) > std::abs(ustop) && du != 0.0_rt) { - du = 0.5_rt * du; - } - } else { - if (u > 0.0_rt) { - while (u + du < ustop && du != 0.0_rt) { - du = 0.5_rt * du; - } - - } else { - while (u + du > ustop && du != 0.0_rt) { - du = 0.5_rt * du; - } - } - } - - if (std::abs(du) < riemann_solver::tol * std::abs(u)) { - finished = true; - } - - } - - rho = 1.0_rt / tau; -} - #endif