diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py index 5ef0d4df7..86c13a50e 100755 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -173,7 +173,7 @@ print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') print(f'RHILL_PRESENT yes') - +print(f'MTINY 1e-12') diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py b/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py index 82850837d..18ef4ce48 100755 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/init_cond.py @@ -38,6 +38,7 @@ sim.param['GR'] = 'NO' sim.param['CHK_CLOSE'] = 'YES' sim.param['RHILL_PRESENT'] = 'YES' +sim.param['MTINY'] = 1.0e-12 sim.param['MU2KG'] = swiftest.MSun sim.param['TU2S'] = swiftest.JD2S diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in index 06edc324b..e9ed6376c 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/param.swiftest.in @@ -33,4 +33,4 @@ ENERGY NO GR NO YARKOVSKY NO YORP NO -MTINY 0.0 +MTINY 1e-12 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in index cd3c71538..fea43297f 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.in @@ -1,33 +1,33 @@ 8 -1 4.9125474498983623693e-11 0.0014751237493860230134 +1 4.9125474498983623693e-11 0.0014751238438755500459 1.6306381826061645943e-05 --0.032433320146471017464 0.30732647407569840814 0.0280888997405028297 --0.033622812072399158034 -0.0019305604712619159943 0.0029264451427202888868 -2 7.243452483873646905e-10 0.006759082196678506012 +-0.065841771551149230746 0.30388831943526661838 0.030872485461978960153 +-0.033141166233939436947 -0.0049297226604189817514 0.0026371811668407158825 +2 7.243452483873646905e-10 0.006759080797928606587 4.0453784346544178454e-05 --0.6608991468450423623 -0.28805695486041710263 0.034183953683804932377 -0.007943018642097033136 -0.018635382188272479886 -0.00071410720992500279457 -3 8.9970113821660187435e-10 0.010044863223462002622 +-0.65269716062695148917 -0.3065765656441301057 0.033456491497379246824 +0.008459831335658639026 -0.0184014319837384685 -0.0007407193515014080928 +3 8.9970113821660187435e-10 0.010044868190633438806 4.25875607065040958e-05 -0.5665449483756358484 -0.84285201543201082597 3.8152874628327130158e-05 -0.0139986033055793102076 0.009533392738922031109 -5.008237574040859916e-07 -4 9.549535102761465607e-11 0.0072467110395904559343 +0.58046286084934750615 -0.8332000042504307258 3.7646553415201541957e-05 +0.013836557832279990782 0.009770187318278569788 -5.1179589633921335467e-07 +4 9.549535102761465607e-11 0.0072467082986392815006 2.265740805092889601e-05 --1.5854600237231359916 0.50600057977052448344 0.049495356229978339224 --0.0037325822023031099417 -0.0121364162752466003825 -0.00016278089870573419053 -5 2.825345908631354893e-07 0.35527078496549785303 +-1.5891417403740180081 0.4938480736359250889 0.049330990309104823244 +-0.0036308073545784510204 -0.012168467501132099878 -0.00016594932370266260858 +5 2.825345908631354893e-07 0.3552707649709459117 0.00046732617030490929307 -4.1105798235203270252 -2.9003636368897538489 -0.07992066204197022239 -0.0042645403767569648595 0.006527961423420942065 -0.00012252307659855749943 -6 8.459715183006415395e-08 0.43765573308845887078 +4.1148395833578952363 -2.8938323061728068453 -0.080043092204059404504 +0.0042549773877191511204 0.006534697671907701254 -0.00012233719535540690457 +6 8.459715183006415395e-08 0.43765596788571493287 0.00038925687730393611812 -6.3549393159832749944 -7.6568459312514027815 -0.11978932080537739446 -0.0039872926987931916337 0.0035567518157804990653 -0.00022047226166396519348 -7 1.2920249163736673626e-08 0.46957395507687206725 +6.3589256477393849565 -7.653288021415167286 -0.12000977499446359442 +0.003985370599203661747 0.0035590677039893160206 -0.00022043610541731448703 +7 1.2920249163736673626e-08 0.46957663585116591335 0.00016953449859497231466 -14.81940372833062014 13.046490834898889943 -0.14356031024960910769 --0.002623943559850705834 0.002775224845039696818 4.4157032104701469965e-05 -8 1.5243589003230834323e-08 0.7813323455417420909 +14.816779495279050138 13.049265812461410263 -0.14351615042000470668 +-0.0026245225263081049631 0.002774730265364384104 4.416262654344997005e-05 +8 1.5243589003230834323e-08 0.7813355837717117843 0.000164587904124493665 -29.563994989459040141 -4.5855881090096284325 -0.5869609072731380994 -0.00046517035659338968386 0.0031282283842968541462 -7.504927375628088796e-05 +29.564459991843019537 -4.5824598513731222837 -0.5870359532621901577 +0.0004648344125208179762 0.0031282868879460171488 -7.5042704502708602616e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in index e42b53b68..79614abb4 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swifter.in @@ -2,35 +2,35 @@ 0 0.00029591220819207775568 0.0 0.0 0.0 0.0 0.0 0.0 -1 4.9125474498983623693e-11 0.0014751237493860230134 +1 4.9125474498983623693e-11 0.0014751238438755500459 1.6306381826061645943e-05 --0.032433320146471017464 0.30732647407569840814 0.0280888997405028297 --0.033622812072399158034 -0.0019305604712619159943 0.0029264451427202888868 -2 7.243452483873646905e-10 0.006759082196678506012 +-0.065841771551149230746 0.30388831943526661838 0.030872485461978960153 +-0.033141166233939436947 -0.0049297226604189817514 0.0026371811668407158825 +2 7.243452483873646905e-10 0.006759080797928606587 4.0453784346544178454e-05 --0.6608991468450423623 -0.28805695486041710263 0.034183953683804932377 -0.007943018642097033136 -0.018635382188272479886 -0.00071410720992500279457 -3 8.9970113821660187435e-10 0.010044863223462002622 +-0.65269716062695148917 -0.3065765656441301057 0.033456491497379246824 +0.008459831335658639026 -0.0184014319837384685 -0.0007407193515014080928 +3 8.9970113821660187435e-10 0.010044868190633438806 4.25875607065040958e-05 -0.5665449483756358484 -0.84285201543201082597 3.8152874628327130158e-05 -0.0139986033055793102076 0.009533392738922031109 -5.008237574040859916e-07 -4 9.549535102761465607e-11 0.0072467110395904559343 +0.58046286084934750615 -0.8332000042504307258 3.7646553415201541957e-05 +0.013836557832279990782 0.009770187318278569788 -5.1179589633921335467e-07 +4 9.549535102761465607e-11 0.0072467082986392815006 2.265740805092889601e-05 --1.5854600237231359916 0.50600057977052448344 0.049495356229978339224 --0.0037325822023031099417 -0.0121364162752466003825 -0.00016278089870573419053 -5 2.825345908631354893e-07 0.35527078496549785303 +-1.5891417403740180081 0.4938480736359250889 0.049330990309104823244 +-0.0036308073545784510204 -0.012168467501132099878 -0.00016594932370266260858 +5 2.825345908631354893e-07 0.3552707649709459117 0.00046732617030490929307 -4.1105798235203270252 -2.9003636368897538489 -0.07992066204197022239 -0.0042645403767569648595 0.006527961423420942065 -0.00012252307659855749943 -6 8.459715183006415395e-08 0.43765573308845887078 +4.1148395833578952363 -2.8938323061728068453 -0.080043092204059404504 +0.0042549773877191511204 0.006534697671907701254 -0.00012233719535540690457 +6 8.459715183006415395e-08 0.43765596788571493287 0.00038925687730393611812 -6.3549393159832749944 -7.6568459312514027815 -0.11978932080537739446 -0.0039872926987931916337 0.0035567518157804990653 -0.00022047226166396519348 -7 1.2920249163736673626e-08 0.46957395507687206725 +6.3589256477393849565 -7.653288021415167286 -0.12000977499446359442 +0.003985370599203661747 0.0035590677039893160206 -0.00022043610541731448703 +7 1.2920249163736673626e-08 0.46957663585116591335 0.00016953449859497231466 -14.81940372833062014 13.046490834898889943 -0.14356031024960910769 --0.002623943559850705834 0.002775224845039696818 4.4157032104701469965e-05 -8 1.5243589003230834323e-08 0.7813323455417420909 +14.816779495279050138 13.049265812461410263 -0.14351615042000470668 +-0.0026245225263081049631 0.002774730265364384104 4.416262654344997005e-05 +8 1.5243589003230834323e-08 0.7813355837717117843 0.000164587904124493665 -29.563994989459040141 -4.5855881090096284325 -0.5869609072731380994 -0.00046517035659338968386 0.0031282283842968541462 -7.504927375628088796e-05 +29.564459991843019537 -4.5824598513731222837 -0.5870359532621901577 +0.0004648344125208179762 0.0031282868879460171488 -7.5042704502708602616e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in index 8bc636061..fea43297f 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in @@ -1,33 +1,33 @@ 8 -201 4.9125474498983623693e-11 0.0014751237493860230134 +1 4.9125474498983623693e-11 0.0014751238438755500459 1.6306381826061645943e-05 --0.032433320146471017464 0.30732647407569840814 0.0280888997405028297 --0.033622812072399158034 -0.0019305604712619159943 0.0029264451427202888868 -2 7.243452483873646905e-10 0.006759082196678506012 +-0.065841771551149230746 0.30388831943526661838 0.030872485461978960153 +-0.033141166233939436947 -0.0049297226604189817514 0.0026371811668407158825 +2 7.243452483873646905e-10 0.006759080797928606587 4.0453784346544178454e-05 --0.6608991468450423623 -0.28805695486041710263 0.034183953683804932377 -0.007943018642097033136 -0.018635382188272479886 -0.00071410720992500279457 -1003 8.9970113821660187435e-10 0.010044863223462002622 +-0.65269716062695148917 -0.3065765656441301057 0.033456491497379246824 +0.008459831335658639026 -0.0184014319837384685 -0.0007407193515014080928 +3 8.9970113821660187435e-10 0.010044868190633438806 4.25875607065040958e-05 -0.5665449483756358484 -0.84285201543201082597 3.8152874628327130158e-05 -0.0139986033055793102076 0.009533392738922031109 -5.008237574040859916e-07 -4 9.549535102761465607e-11 0.0072467110395904559343 +0.58046286084934750615 -0.8332000042504307258 3.7646553415201541957e-05 +0.013836557832279990782 0.009770187318278569788 -5.1179589633921335467e-07 +4 9.549535102761465607e-11 0.0072467082986392815006 2.265740805092889601e-05 --1.5854600237231359916 0.50600057977052448344 0.049495356229978339224 --0.0037325822023031099417 -0.0121364162752466003825 -0.00016278089870573419053 -5 2.825345908631354893e-07 0.35527078496549785303 +-1.5891417403740180081 0.4938480736359250889 0.049330990309104823244 +-0.0036308073545784510204 -0.012168467501132099878 -0.00016594932370266260858 +5 2.825345908631354893e-07 0.3552707649709459117 0.00046732617030490929307 -4.1105798235203270252 -2.9003636368897538489 -0.07992066204197022239 -0.0042645403767569648595 0.006527961423420942065 -0.00012252307659855749943 -6 8.459715183006415395e-08 0.43765573308845887078 +4.1148395833578952363 -2.8938323061728068453 -0.080043092204059404504 +0.0042549773877191511204 0.006534697671907701254 -0.00012233719535540690457 +6 8.459715183006415395e-08 0.43765596788571493287 0.00038925687730393611812 -6.3549393159832749944 -7.6568459312514027815 -0.11978932080537739446 -0.0039872926987931916337 0.0035567518157804990653 -0.00022047226166396519348 -7 1.2920249163736673626e-08 0.46957395507687206725 +6.3589256477393849565 -7.653288021415167286 -0.12000977499446359442 +0.003985370599203661747 0.0035590677039893160206 -0.00022043610541731448703 +7 1.2920249163736673626e-08 0.46957663585116591335 0.00016953449859497231466 -14.81940372833062014 13.046490834898889943 -0.14356031024960910769 --0.002623943559850705834 0.002775224845039696818 4.4157032104701469965e-05 -8 1.5243589003230834323e-08 0.7813323455417420909 +14.816779495279050138 13.049265812461410263 -0.14351615042000470668 +-0.0026245225263081049631 0.002774730265364384104 4.416262654344997005e-05 +8 1.5243589003230834323e-08 0.7813355837717117843 0.000164587904124493665 -29.563994989459040141 -4.5855881090096284325 -0.5869609072731380994 -0.00046517035659338968386 0.0031282283842968541462 -7.504927375628088796e-05 +29.564459991843019537 -4.5824598513731222837 -0.5870359532621901577 +0.0004648344125208179762 0.0031282868879460171488 -7.5042704502708602616e-05 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in index e7424ed6f..7c1b15bd6 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/tp.in @@ -1,49 +1,49 @@ 16 101 --0.032410259440139366216 0.30734953478203003163 0.0280888997405028297 --0.031019262600390007378 -0.0019305604712619159943 0.0029264451427202888868 +-0.0658187108448175795 0.30391138014159824188 0.030872485461978960153 +-0.030537616761930286291 -0.0049297226604189817514 0.0026371811668407158825 102 --0.032456380852802668713 0.30730341336936678465 0.0280888997405028297 --0.03622636154440830869 -0.0019305604712619159943 0.0029264451427202888868 +-0.065864832257480881994 0.3038652587289349949 0.030872485461978960153 +-0.035744715705948587603 -0.0049297226604189817514 0.0026371811668407158825 103 --0.6608419365545701307 -0.28799974456994492655 0.034183953683804932377 -0.014290252347442427075 -0.018635382188272479886 -0.00071410720992500279457 +-0.6526399503364792576 -0.30651935535365792962 0.033456491497379246824 +0.014807065041004032965 -0.0184014319837384685 -0.0007407193515014080928 104 --0.6609563571355145939 -0.2881141651508892787 0.034183953683804932377 -0.0015957849367516391964 -0.018635382188272479886 -0.00071410720992500279457 +-0.65275437091742372075 -0.30663377593460228177 0.033456491497379246824 +0.0021125976303132450868 -0.0184014319837384685 -0.0007407193515014080928 105 -0.5666051762815753623 -0.8427917875260713121 3.8152874628327130158e-05 -0.020893043539853186491 0.009533392738922031109 -5.008237574040859916e-07 +0.58052308875528702004 -0.8331397763444912119 3.7646553415201541957e-05 +0.020730998066553867065 0.009770187318278569788 -5.1179589633921335467e-07 106 -0.56648472046969633453 -0.84291224333795033985 3.8152874628327130158e-05 -0.0071041630713054347915 0.009533392738922031109 -5.008237574040859916e-07 +0.58040263294340799227 -0.83326023215637023966 3.7646553415201541957e-05 +0.0069421175980061153657 0.009770187318278569788 -5.1179589633921335467e-07 107 --1.5854279813093821172 0.50603262218427835784 0.049495356229978339224 --0.00065310310407921696313 -0.0121364162752466003825 -0.00016278089870573419053 +-1.5891096979602641337 0.49388011604967890777 0.049330990309104823244 +-0.00055132825635455804184 -0.012168467501132099878 -0.00016594932370266260858 108 --1.585492066136889866 0.50596853735677060904 0.049495356229978339224 --0.0068120613005270029203 -0.0121364162752466003825 -0.00016278089870573419053 +-1.5891737827877718825 0.49381603122217127 0.049330990309104823244 +-0.0067102864528023435653 -0.012168467501132099878 -0.00016594932370266260858 109 -4.1112407225284242074 -2.8997027378816566667 -0.07992066204197022239 -0.041146757578404392908 0.006527961423420942065 -0.00012252307659855749943 +4.1155004823659924185 -2.893171407164709663 -0.080043092204059404504 +0.0411371945893665783 0.006534697671907701254 -0.00012233719535540690457 110 -4.109918924512229843 -2.901024535897851031 -0.07992066204197022239 --0.032617676824890466658 0.006527961423420942065 -0.00012252307659855749943 +4.114178684349798054 -2.8944932051809040274 -0.080043092204059404504 +-0.032627239813928281265 0.006534697671907701254 -0.00012233719535540690457 111 -6.3554898083384054885 -7.6562954388962722874 -0.11978932080537739446 -0.026100461582326782428 0.0035567518157804990653 -0.00022047226166396519348 +6.3594761400945154506 -7.652737529060036792 -0.12000977499446359442 +0.02609853948273724994 0.0035590677039893160206 -0.00022043610541731448703 112 -6.3543888236281445003 -7.6573964236065332756 -0.11978932080537739446 --0.018125876184740395691 0.0035567518157804990653 -0.00022047226166396519348 +6.3583751553842544624 -7.65383851377029778 -0.12000977499446359442 +-0.01812779828432992818 0.0035590677039893160206 -0.00022043610541731448703 113 -14.819643486317822578 13.046730592886092381 -0.14356031024960910769 -0.010470819978811187617 0.002775224845039696818 4.4157032104701469965e-05 +14.817019253266252576 13.049505570448612701 -0.14351615042000470668 +0.010470241012353788054 0.002774730265364384104 4.416262654344997005e-05 114 -14.819163970343417702 13.046251076911687505 -0.14356031024960910769 --0.015718707098512599285 0.002775224845039696818 4.4157032104701469965e-05 +14.8165397372918477 13.049026054474207825 -0.14351615042000470668 +-0.015719286064969997113 0.002774730265364384104 4.416262654344997005e-05 115 -29.56422775190525698 -4.58535534656341337 -0.5869609072731380994 -0.014900806169871520443 0.0031282283842968541462 -7.504927375628088796e-05 +29.564692754289236376 -4.5822270889269072214 -0.5870359532621901577 +0.014900470225798949711 0.0031282868879460171488 -7.5042704502708602616e-05 116 -29.563762227012823303 -4.585820871455843495 -0.5869609072731380994 --0.013970465456684741726 0.0031282283842968541462 -7.504927375628088796e-05 +29.564227229396802699 -4.582692613819337346 -0.5870359532621901577 +-0.013970801400757312458 0.0031282868879460171488 -7.5042704502708602616e-05 diff --git a/src/io/io.f90 b/src/io/io.f90 index 82c87dfc7..7c50242cb 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -259,10 +259,10 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Determine if the GR flag is set correctly for this integrator select case(integrator) - case(WHM, RMVS) + case(WHM, RMVS, SYMBA) write(*,*) "GR = ", self%lgr case default - write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' + if (self%lgr) write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' end select end associate diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 4ed7cf3fe..4c6bccc72 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -10,9 +10,9 @@ program swiftest_driver implicit none class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated - type(swiftest_parameters) :: param + class(swiftest_parameters), allocatable :: param !! Run configuration parameters integer(I4B) :: integrator !! Integrator type code (see swiftest_globals for symbolic names) - character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters + character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters integer(I4B) :: ierr !! I/O error code integer(I8B) :: iloop !! Loop counter integer(I8B) :: idump !! Dump cadence counter @@ -31,6 +31,12 @@ program swiftest_driver end if !$ start_wall_time = omp_get_wtime() !> Read in the user-defined parameters file and the initial conditions of the system + select case(integrator) + case(symba) + allocate(symba_parameters :: param) + case default + allocate(swiftest_parameters :: param) + end select param%integrator = integrator call setup_construct_system(nbody_system, param) call param%read_from_file(param_file_name) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 0a6dcb290..98108e9df 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -89,10 +89,12 @@ module symba_classes type(symba_particle_info), dimension(:), allocatable :: info contains private - procedure, public :: discard => symba_discard_pl !! Process massive body discards - procedure, public :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other - procedure, public :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure, public :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle + procedure, public :: discard => symba_discard_pl !! Process massive body discards + procedure, public :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure, public :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure, public :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle + procedure, public :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen + procedure, public :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods end type symba_pl !******************************************************************************************************************************** @@ -105,9 +107,11 @@ module symba_classes integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step contains private - procedure, public :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body - procedure, public :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles - procedure, public :: setup => symba_setup_tp !! Constructor method - Allocates space for number of particle + procedure, public :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body + procedure, public :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles + procedure, public :: setup => symba_setup_tp !! Constructor method - Allocates space for number of particle + procedure, public :: sort => symba_util_sort_tp !! Sorts body arrays by a sortable componen + procedure, public :: rearrange => symba_util_sort_rearrange_tp !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods end type symba_tp !******************************************************************************************************************************** @@ -390,5 +394,32 @@ module subroutine symba_util_resize_pltpenc(self, nrequested) integer(I4B), intent(in) :: nrequested !! New size of list needed end subroutine symba_util_resize_pltpenc + module subroutine symba_util_sort_pl(self, sortby, ascending) + implicit none + class(symba_pl), intent(inout) :: self !! Symba massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + end subroutine symba_util_sort_pl + + module subroutine symba_util_sort_tp(self, sortby, ascending) + implicit none + class(symba_tp), intent(inout) :: self !! Swiftest test particle object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + end subroutine symba_util_sort_tp + + module subroutine symba_util_sort_rearrange_pl(self, ind) + implicit none + class(symba_pl), intent(inout) :: self !! Symba massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + end subroutine symba_util_sort_rearrange_pl + + module subroutine symba_util_sort_rearrange_tp(self, ind) + implicit none + class(symba_tp), intent(inout) :: self !! Symba massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + end subroutine symba_util_sort_rearrange_tp + + end interface end module symba_classes \ No newline at end of file diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index bebb225b5..acc3aabf9 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -51,6 +51,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms call random_seed(size = nseeds) if (allocated(param%seed)) deallocate(param%seed) allocate(param%seed(nseeds)) + rewind(unit) do read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line line_trim = trim(adjustl(line)) @@ -121,8 +122,10 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms return end if end associate - return + iostat = 0 + + return end subroutine symba_io_param_reader module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 5ac26c220..9efb37e9a 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -104,6 +104,7 @@ module subroutine symba_setup_system(self, param) !! author: David A. Minton !! !! Initialize an SyMBA nbody system from files and sets up the planetocentric structures. + !! This subroutine will also sort the massive bodies in descending order by mass !! implicit none ! Arguments @@ -121,6 +122,7 @@ module subroutine symba_setup_system(self, param) call system%plplenc_list%setup(1) select type(pl => system%pl) class is (symba_pl) + call pl%sort("mass", ascending=.false.) select type(param) class is (symba_parameters) pl%lmtiny(:) = pl%Gmass(:) > param%MTINY diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 031ae4ae5..24a59a9cc 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -65,5 +65,166 @@ module subroutine symba_util_resize_pltpenc(self, nrequested) return end subroutine symba_util_resize_pltpenc + module subroutine symba_util_sort_pl(self, sortby, ascending) + !! author: David A. Minton + !! + !! Sort a Swiftest test particle object in-place. + !! sortby is a string indicating which array component to sort. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! Symba massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + ! Internals + integer(I4B), dimension(self%nbody) :: ind + + associate(pl => self, npl => self%nbody) + select case(sortby) + case("nplenc") + if (ascending) then + call util_sort(pl%nplenc(1:npl), ind(1:npl)) + else + call util_sort(-pl%nplenc(1:npl), ind(1:npl)) + end if + case("ntpenc") + if (ascending) then + call util_sort(pl%ntpenc(1:npl), ind(1:npl)) + else + call util_sort(-pl%ntpenc(1:npl), ind(1:npl)) + end if + case("levelg") + if (ascending) then + call util_sort(pl%levelg(1:npl), ind(1:npl)) + else + call util_sort(-pl%levelg(1:npl), ind(1:npl)) + end if + case("levelm") + if (ascending) then + call util_sort(pl%levelm(1:npl), ind(1:npl)) + else + call util_sort(-pl%levelm(1:npl), ind(1:npl)) + end if + case("peri") + if (ascending) then + call util_sort(pl%peri(1:npl), ind(1:npl)) + else + call util_sort(-pl%peri(1:npl), ind(1:npl)) + end if + case("atp") + if (ascending) then + call util_sort(pl%atp(1:npl), ind(1:npl)) + else + call util_sort(-pl%atp(1:npl), ind(1:npl)) + end if + case default + call util_sort_pl(pl, sortby, ascending) + end select + call pl%rearrange(ind) + end associate + + return + end subroutine symba_util_sort_pl + + module subroutine symba_util_sort_tp(self, sortby, ascending) + !! author: David A. Minton + !! + !! Sort a Swiftest test particle object in-place. + !! sortby is a string indicating which array component to sort. + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! Swiftest test particle object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + ! Internals + integer(I4B), dimension(self%nbody) :: ind + + associate(tp => self, ntp => self%nbody) + select case(sortby) + case("nplenc") + if (ascending) then + call util_sort(tp%nplenc(1:ntp), ind(1:ntp)) + else + call util_sort(-tp%nplenc(1:ntp), ind(1:ntp)) + end if + case("levelg") + if (ascending) then + call util_sort(tp%levelg(1:ntp), ind(1:ntp)) + else + call util_sort(-tp%levelg(1:ntp), ind(1:ntp)) + end if + case("levelm") + if (ascending) then + call util_sort(tp%levelm(1:ntp), ind(1:ntp)) + else + call util_sort(-tp%levelm(1:ntp), ind(1:ntp)) + end if + case default + call util_sort_tp(tp, sortby, ascending) + end select + call tp%rearrange(ind) + end associate + + return + end subroutine symba_util_sort_tp + + module subroutine symba_util_sort_rearrange_pl(self, ind) + !! author: David A. Minton + !! + !! Rearrange SyMBA massive body structure in-place from an index list. + !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! Symba massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + ! Internals + class(symba_pl), allocatable :: pl_sorted !! Temporary holder for sorted body + integer(I4B) :: i, j + + associate(pl => self, npl => self%nbody) + call util_sort_rearrange_pl(pl,ind) + allocate(pl_sorted, source=self) + pl%lcollision(1:npl) = pl_sorted%lcollision(ind(1:npl)) + pl%lencounter(1:npl) = pl_sorted%lencounter(ind(1:npl)) + pl%nplenc(1:npl) = pl_sorted%nplenc(ind(1:npl)) + pl%ntpenc(1:npl) = pl_sorted%ntpenc(ind(1:npl)) + pl%levelg(1:npl) = pl_sorted%levelg(ind(1:npl)) + pl%levelm(1:npl) = pl_sorted%levelm(ind(1:npl)) + pl%isperi(1:npl) = pl_sorted%isperi(ind(1:npl)) + pl%peri(1:npl) = pl_sorted%peri(ind(1:npl)) + pl%atp(1:npl) = pl_sorted%atp(ind(1:npl)) + pl%info(1:npl) = pl_sorted%info(ind(1:npl)) + pl%kin(1:npl) = pl_sorted%kin(ind(1:npl)) + do i = 1, npl + do j = 1, pl%kin(i)%nchild + pl%kin(i)%child(j) = ind(pl%kin(i)%child(j)) + end do + end do + deallocate(pl_sorted) + end associate + return + end subroutine symba_util_sort_rearrange_pl + + module subroutine symba_util_sort_rearrange_tp(self, ind) + !! author: David A. Minton + !! + !! Rearrange SyMBA test particle object in-place from an index list. + !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! Symba massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + ! Internals + class(symba_tp), allocatable :: tp_sorted !! Temporary holder for sorted body + + associate(tp => self, ntp => self%nbody) + call util_sort_rearrange_tp(tp,ind) + allocate(tp_sorted, source=self) + tp%nplenc(1:ntp) = tp_sorted%nplenc(ind(1:ntp)) + tp%levelg(1:ntp) = tp_sorted%levelg(ind(1:ntp)) + tp%levelm(1:ntp) = tp_sorted%levelm(ind(1:ntp)) + deallocate(tp_sorted) + end associate + return + end subroutine symba_util_sort_rearrange_tp end submodule s_symba_util \ No newline at end of file diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index fdaded124..044e56428 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -12,7 +12,7 @@ module subroutine util_sort_body(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(:), allocatable :: ind + integer(I4B), dimension(self%nbody) :: ind associate(body => self, n => self%nbody) select case(sortby) @@ -73,7 +73,7 @@ module subroutine util_sort_pl(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(:), allocatable :: ind + integer(I4B), dimension(self%nbody) :: ind associate(pl => self, npl => self%nbody) select case(sortby) @@ -140,7 +140,7 @@ module subroutine util_sort_tp(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(:), allocatable :: ind + integer(I4B), dimension(self%nbody) :: ind associate(tp => self, ntp => self%nbody) select case(sortby)