diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8cdd9d41..33a56bd4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,4 +28,4 @@ jobs: sudo apt-get install -y python3-pip python3 -m pip install --upgrade pip python3 -m pip install "$PWD"[test] - python3 -m pytest + python3 -m pytest -rN diff --git a/src/pydrex/cli.py b/src/pydrex/cli.py index 8b8d3055..911d709c 100644 --- a/src/pydrex/cli.py +++ b/src/pydrex/cli.py @@ -7,6 +7,7 @@ import numpy as np from pydrex import exceptions as _err +from pydrex import io as _io from pydrex import logger as _log from pydrex import minerals as _minerals from pydrex import stats as _stats @@ -31,7 +32,9 @@ def __call__(self): if args.range is None: i_range = None else: - i_range = range(*(int(s) for s in args.range.split(":"))) + start, stop_ex, step = (int(s) for s in args.range.split(":")) + # Make command line start:stop:step stop-inclusive, it's more intuitive. + i_range = range(start, stop_ex + step, step) density_kwargs = {"kernel": args.kernel} if args.smoothing is not None: @@ -58,6 +61,7 @@ def __call__(self): i_range=i_range, density=args.density, savefile=args.out, + strains=_io.read_scsv(args.scsv).strain, **density_kwargs, ) except (argparse.ArgumentError, ValueError, _err.Error) as e: @@ -73,6 +77,15 @@ def _get_args(self) -> argparse.Namespace: help="range of strain indices to be plotted, in the format start:stop:step", default=None, ) + parser.add_argument( + "-f", + "--scsv", + help=( + "path to SCSV file with a column named 'strain'" + + " that lists shear strain percentages for each strain index" + ), + default=None, + ) parser.add_argument( "-p", "--postfix", diff --git a/src/pydrex/data/rng/seeds.scsv b/src/pydrex/data/rng/seeds.scsv new file mode 100644 index 00000000..8c030997 --- /dev/null +++ b/src/pydrex/data/rng/seeds.scsv @@ -0,0 +1,1010 @@ +--- +schema: + delimiter: ',' + missing: '-' + fields: + - name: seeds + type: integer + fill: 0 +--- +seeds +101 +121 +1218 +3052 +4232 +4557 +5675 +8236 +8711 +10113 +10325 +12198 +12467 +12762 +13178 +13515 +14756 +15181 +16535 +16608 +16913 +19182 +20342 +20442 +20519 +20692 +20893 +20909 +21232 +21400 +21590 +24094 +26169 +26974 +28635 +31173 +31719 +33355 +34472 +34851 +36199 +43410 +43879 +44099 +44409 +44886 +45928 +46708 +46908 +47611 +48561 +49754 +50018 +50036 +50296 +50569 +50669 +51406 +51623 +52535 +53478 +54061 +54882 +55276 +56346 +56394 +57961 +59776 +60397 +60411 +60955 +63519 +63566 +64142 +64206 +67640 +70280 +70618 +71666 +71838 +72127 +74830 +75257 +75269 +75509 +76437 +76509 +77002 +77468 +79002 +79073 +80782 +80962 +83466 +84510 +85192 +86036 +86446 +87099 +87126 +90646 +90949 +91002 +91827 +92409 +92420 +94664 +95791 +97951 +98873 +100710 +101899 +104164 +104963 +105813 +106577 +107646 +107784 +108242 +109640 +110784 +111304 +111478 +112182 +113363 +113535 +113663 +114704 +115470 +116042 +116335 +117661 +118109 +120051 +120765 +121729 +124533 +126017 +126657 +129519 +129878 +130272 +130950 +131319 +132147 +132333 +133142 +134203 +134561 +134679 +136585 +136866 +137027 +138114 +138160 +139623 +141670 +142614 +142630 +142677 +143678 +144881 +144972 +145400 +145478 +145637 +146239 +147742 +148686 +149713 +150432 +152078 +152272 +152827 +153859 +154469 +154522 +154661 +157064 +159107 +161425 +161575 +162130 +162194 +163472 +163863 +165201 +165590 +166203 +167019 +168486 +170180 +171159 +172981 +173861 +174256 +175802 +176530 +176625 +176661 +176762 +177884 +177902 +178679 +179030 +179426 +179772 +180748 +182886 +183850 +184169 +185118 +185422 +185920 +186416 +187608 +188228 +191526 +191953 +193226 +194848 +195211 +199081 +199740 +200064 +200969 +201228 +201897 +202641 +203371 +205585 +206823 +207006 +207651 +207694 +209018 +212606 +213511 +213662 +214205 +215936 +216028 +217889 +218486 +218898 +219503 +220767 +221923 +224888 +227310 +227658 +227774 +229028 +229228 +229427 +232214 +233523 +234135 +234590 +235239 +236275 +236990 +237639 +240855 +241246 +241308 +241519 +242014 +244336 +245199 +247561 +247885 +250933 +251730 +254091 +254312 +254421 +254663 +255256 +255659 +255767 +255813 +255850 +257203 +257382 +259837 +260294 +260583 +261521 +261725 +261746 +262517 +262591 +262643 +263132 +263585 +263664 +264112 +266552 +267007 +267096 +267664 +269002 +271054 +271783 +271879 +272072 +272103 +273326 +275228 +275445 +276563 +278368 +278855 +278877 +280582 +280807 +283301 +284463 +284933 +287994 +289416 +291901 +292819 +295332 +295940 +296427 +297562 +298877 +300906 +301643 +302043 +302942 +304463 +305273 +306387 +310707 +311713 +313520 +314883 +315141 +317174 +320123 +320152 +320651 +322318 +323683 +324677 +325364 +326698 +327262 +327333 +327578 +329023 +329785 +332494 +332962 +334043 +335832 +336802 +336919 +338138 +338987 +343049 +345528 +346444 +346585 +347563 +349501 +352775 +355260 +356715 +357839 +358841 +358865 +358948 +359893 +360175 +360396 +361291 +362095 +362290 +363809 +364865 +365604 +370152 +371380 +372149 +373716 +373928 +374838 +374909 +378247 +379464 +379554 +381361 +382202 +382449 +382599 +382866 +383357 +383533 +383853 +385205 +385469 +386029 +388625 +389027 +390437 +391075 +391332 +391520 +393397 +393906 +395389 +398151 +398261 +398460 +399329 +399601 +402816 +403551 +403990 +404522 +405987 +412442 +412630 +412656 +413271 +413467 +415170 +415344 +415773 +415936 +418239 +418388 +418701 +418778 +419633 +420532 +421043 +421263 +423522 +423915 +425898 +427400 +427946 +428117 +429464 +429593 +433933 +433983 +434484 +435212 +435348 +435604 +437148 +437283 +438641 +439831 +441009 +441581 +441979 +442138 +442804 +444584 +447019 +447986 +449235 +449525 +450820 +450822 +452296 +452971 +454420 +454661 +455744 +455876 +456950 +458972 +459948 +462615 +463811 +464004 +465753 +466349 +467139 +467531 +467584 +467710 +468484 +470558 +470682 +470934 +471145 +471826 +471955 +472864 +473529 +475202 +476465 +476753 +477252 +477384 +478311 +479818 +481452 +482138 +483000 +483973 +484116 +484217 +484872 +485865 +487222 +491172 +491999 +493057 +493696 +493824 +494459 +498480 +499571 +500049 +500346 +500445 +502323 +502950 +504119 +504713 +504941 +506016 +507663 +511945 +512633 +513989 +516960 +517294 +517801 +517834 +522321 +523270 +523637 +526204 +526431 +526893 +526906 +526973 +528611 +530070 +530178 +530407 +530639 +530947 +531128 +531275 +531590 +533010 +533197 +533276 +534678 +537643 +538678 +539939 +540030 +542130 +542633 +544457 +546530 +548646 +548717 +548805 +552227 +553607 +557117 +558482 +558904 +559903 +560888 +561368 +563245 +563294 +563520 +564397 +564819 +565582 +566792 +567812 +569815 +572072 +572878 +576465 +576587 +577557 +579867 +580260 +580682 +581478 +582940 +583940 +584145 +586959 +588493 +588532 +589990 +590776 +593155 +593430 +593632 +594885 +595898 +598620 +599420 +601715 +601872 +602178 +602591 +602611 +602888 +605984 +606798 +607131 +607931 +608085 +610722 +613000 +613276 +613315 +614909 +615460 +615779 +617764 +617963 +619214 +620203 +622414 +623978 +624052 +624286 +624650 +625527 +626657 +626716 +626825 +630264 +630760 +631079 +631264 +631575 +631733 +633212 +634957 +635967 +636516 +641406 +641606 +642197 +643329 +643397 +643632 +644194 +644308 +644785 +645815 +647888 +648320 +649396 +650704 +652522 +654155 +654313 +654704 +655820 +658800 +658832 +659286 +660245 +660659 +661129 +664010 +664975 +667304 +667599 +668089 +668210 +668477 +668915 +669689 +670024 +673294 +673510 +676008 +678502 +679143 +679305 +679762 +680229 +682241 +682719 +684165 +685298 +686964 +690741 +691328 +691402 +692105 +693619 +696799 +700232 +700421 +700540 +701239 +701297 +701498 +702498 +702840 +702962 +702971 +703363 +703402 +704854 +706082 +706428 +707107 +707957 +711738 +712288 +712397 +712757 +713260 +714453 +714743 +715834 +716660 +716975 +718168 +718628 +719645 +720116 +721325 +721367 +722486 +724389 +724722 +725961 +726835 +732226 +733141 +736511 +737244 +741208 +741645 +742280 +743755 +745060 +745957 +746212 +746979 +747513 +749077 +749419 +749989 +751348 +755830 +756270 +757168 +757597 +759115 +760318 +761602 +762269 +763085 +763215 +763560 +764092 +765464 +766537 +766679 +767822 +768656 +769114 +771565 +771988 +775845 +775887 +776480 +777256 +777739 +779012 +779863 +779897 +779912 +780362 +780999 +781265 +781550 +781796 +782246 +783556 +787035 +787654 +788261 +788442 +788591 +789326 +791016 +791183 +792031 +793039 +794748 +795684 +795955 +796421 +797516 +799256 +799777 +801081 +803588 +807751 +807914 +808120 +809399 +814791 +814922 +815033 +815346 +817360 +818647 +818951 +820052 +820449 +826627 +827412 +830216 +830482 +830580 +831032 +832097 +832245 +839261 +839381 +839789 +841781 +843469 +844115 +844932 +845283 +845317 +846331 +847041 +847373 +847699 +848528 +851434 +851465 +851466 +851774 +852268 +852789 +853654 +855316 +856162 +856547 +858407 +859504 +860219 +860654 +862215 +862307 +864936 +865908 +866458 +866753 +867911 +870193 +871661 +873418 +877441 +878009 +878255 +880540 +881472 +881988 +882393 +883930 +884255 +885117 +885118 +887174 +887953 +888765 +889304 +889352 +890564 +893397 +894275 +894633 +896104 +897298 +898800 +899154 +900886 +901107 +901859 +901993 +902375 +907005 +907141 +909813 +911359 +912036 +912699 +912766 +913740 +914224 +915912 +916088 +916528 +919320 +922633 +924463 +926970 +927283 +927983 +929035 +932962 +933185 +933715 +933810 +934681 +934725 +934727 +940356 +941238 +942410 +942480 +943924 +945225 +945769 +946655 +947763 +948215 +948388 +949448 +951693 +952402 +952635 +952872 +953753 +959636 +960626 +963786 +964389 +965439 +965750 +966015 +967207 +967844 +968522 +968848 +968990 +969958 +969985 +971247 +971737 +972461 +972976 +974358 +974404 +974628 +977057 +977819 +977995 +978427 +978718 +980127 +980330 +981505 +982017 +982332 +983276 +983368 +984405 +984453 +985561 +985918 +986581 +986698 +986809 +988366 +989661 +989866 +990988 +992476 +992937 +993197 +993759 +994253 +994818 +995242 diff --git a/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv new file mode 100644 index 00000000..ff5180d2 --- /dev/null +++ b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv @@ -0,0 +1,124 @@ +--- +# Digitized data from Figure 5 in Kaminski & Ribe 2001, https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9 +# This data records the mean angle of the a-axes of an olivine aggregate from the horizontal shear direction. +# The flow field is a static simple shear and the aggregate starts with random (i.e. isotropically distributed) grain orientations. + +schema: + delimiter: ',' + missing: '-' + fields: + - name: equivalent_strain_M0 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M0 + type: float + fill: NaN + unit: degree + - name: equivalent_strain_M50 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M50 + type: float + fill: NaN + unit: degree + - name: equivalent_strain_M200 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M200 + type: float + fill: NaN + unit: degree +--- +equivalent_strain_M0, angle_M0,equivalent_strain_M50, angle_M50,equivalent_strain_M200, angle_M200 + 0.3076923076923084, 44.93959731543627, 1.098613804496157, 44.77430481378372, 1.5111685699920994, 44.72085812156877 + 1.9237233354880399,44.718921206880175, 3.367665014723837, 44.37236387335087, 3.7802197802197792, 44.29219383502844 + 4.19277454571572,44.382124404232606, 5.5748330101271275, 44.00933456156514, 5.911752735282144, 43.91091652503244 + 6.4618257559434,44.026404780609305, 7.390073978309269, 43.76427140580321, 7.390073978309269, 43.64428903552475 + 8.730876966171083, 43.71627845769182, 10.381096028154849, 43.30748138181449, 10.277957336780863, 43.0062852881645 + 10.999928176398763, 43.3868987470251, 13.243194713782946, 42.80741200283247, 13.475256769374413, 42.10306095185686 + 13.330862601450834,43.052174367136885, 15.385908209437625, 42.32208498611974, 15.744307979602093, 41.497695356360985 + 15.538030596854123, 42.7215155332013, 16.98197227608992, 42.024527036765505, 17.291388350211875, 41.01667512642642 + 17.78645406880701,42.402471503142664, 19.972994325935495, 41.529599759366846, 24.40795805501687, 38.67374665953426 + 20.07613301730948,42.073610733697606, 23.17029375852905, 40.951884646476046, 26.573870573870565, 37.89593367548361 + 22.34518422753716, 41.76929181271859, 26.47073188249658, 40.37596926913943, 27.811534870358393, 37.43520137361432 + 24.61423543776484, 41.4584283988153, 29.564892623716148, 39.7928549495861, 34.30927242691948, 35.035008593634714 + 26.88328664799252,41.134475999063454, 32.65905336493571, 39.21693957224949, 36.47518494577317, 34.24606997161278 + 29.1523378582202, 40.80725135284947, 35.85635279752926, 38.62302683937109, 37.60971055088701, 33.81773290971867 + 31.421389068447887, 40.4898434460219, 39.05365223012281, 38.00031833762587, 43.79803203332614, 31.350095494324965 + 33.69044027867557, 40.16916329273219, 42.147812971342375, 37.417204018072546, 45.55138978668389, 30.569210234429306 + 35.95949148890325, 39.85502763236676, 48.542411836529475, 36.157389130148694, 46.685915391797735, 30.050286482974958 + 38.22854269913093, 39.54089197200133, 51.636572577749035, 35.545479041728534, 52.15226603461896, 27.35218293133803 + 40.497593909358606, 39.22021181871162, 54.7307333189686, 34.89037530000813, 53.69934640522875, 26.551943372427214 + 42.766645119586286, 38.91262065127047, 57.721755368814186, 34.2766654760338, 54.83387201034259, 25.970885893507234 + 45.035696329813966,38.611573976753604, 60.712777418659755, 33.63775935430098, 59.99080657904186, 23.145301073449442 + 47.304747540041646,38.294166069926035, 63.80693815987932, 32.91786513263021, 61.537886949651636, 22.319993769319733 + 49.573798750269326, 37.98984714894703, 66.7979602097249, 32.24596385907082, 62.5692738633915, 21.759504696733202 + 51.842849960497006, 37.68552822796802, 69.78898225957047, 31.53566822702232, 67.9030176173033, 18.90563831796692 + 54.11190117072469, 37.38120930698901, 72.78000430941606, 30.860167482354576, 68.49974861739567, 18.304012432713485 + 56.40387209014658, 37.08125338129285, 75.77102635926164, 30.16067026363114, 69.99525964231844, 17.81551278229403 + 58.65000359118005, 36.78893269734169, 78.65890971773322, 29.46237286861049, 75.77102635926164, 14.781672848110055 + 61.00156575450691, 36.47414258768383, 81.75307045895279, 28.728080762506302, 77.31810672987142, 14.002644458230606 + 63.18810601163541, 36.19665608769437, 84.64095381742439, 28.022584425268942, 78.45263233498525, 13.396733488324372 + 65.45715722186308, 35.90869839902606, 87.52883717589599, 27.302690203598168, 83.81584428643251, 10.6422810735031 + 67.72620843209077, 35.61092397097133, 90.41672053436757, 26.58279598192739, 85.36292465704229, 9.827257972540117 + 69.99525964231844,35.316421789378744, 93.30460389283917, 25.891697529123448, 86.49745026215612, 9.24534347668957 + 72.26431085254613, 35.03828084009685, 96.29562594268475, 25.11601150527319, 91.86066221360338, 6.348626251395267 + 74.5333620627738, 34.75359539789068, 99.28664799253032, 24.387118605831528, 93.40774258421315, 5.530389336942683 + 76.80241327300149, 34.46890995568451, 102.07139265962793, 23.6402283508481, 94.439129497953, 5.00396668734593 + 79.07146448322918,34.177680020554064, 105.06241470947352, 22.868141798106194, 100.52431228901816, 2.361627669167966 + 81.34051569345685,33.909355810658596, 107.9502980679451, 22.097854980918466, 102.79336349924583, 1.6744559121185958 + 83.60956690368454, 33.64103160076312, 110.79692594986709, 21.33476710594745, 104.13416648710763, 1.3505035123667497 + 85.87861811391221, 33.36289065148123, 113.62292609351431, 20.56268055320554, 113.7260647848883, 0.08850712680145989 + 88.1476693241399, 33.07820520927506, 116.40767076061192, 19.794193471571987, 115.99511599511597, 0.06560140156648941 + 90.41672053436757,32.816425492303864, 119.08927673633553, 19.04550348103438, 117.74847374847371, 0.02469832078973866 + 92.68577174459526, 32.55137352887054, 121.97716009480712, 18.232023010546406, 120.7394957983193,-0.08928493097479873 + 94.99607843137254, 32.28075874645157, 124.76190476190473, 17.44913804447944, 123.00854700854697,-0.08928493097479873 + 97.22387416505062, 32.01799735554174, 127.54664942900234, 16.641056780653997, 125.27759821877467,-0.08928493097479873 + 99.49292537527829, 31.75948988503269, 130.12511671335199, 15.901365467887274, 127.54664942900234,-0.10564616328550613 + 101.76197658550598,31.520615893296476, 133.01300007182357, 15.037492401882346, 129.81570063923002, -0.1612743531418772 + 104.03102779573365,31.232658204628166, 135.79774473892118, 14.207814311406779, 132.08475184945772, -0.1449131208311769 + 106.30007900596134,30.987239719967675, 138.4793507146448, 13.41593066756893, 133.32241614594554,-0.08928493097479873 + 108.56913021618902,30.725460002996485, 141.16095669036844, 12.559256543780705, 144.255117431588, -0.1612743531418772 + 110.8381814264167,30.473497025411714, 143.73942397471805, 11.728978541453742, 146.52416864181566,-0.14818536729332266 + 113.0247216835452,30.236259156906574, 146.4210299504417, 10.89810062727539, 148.79321985204336,-0.15472986021759993 + 115.37628384687206, 29.96302657731789, 149.1026359261653, 10.097218305666658, 151.06227106227104,-0.08274043805052145 + 117.64533505709974, 29.74051381789238, 151.88738059326292, 9.283737835178684, 153.3313222724987,-0.08928493097479873 + 119.91438626732742,29.495095333231887, 154.67212526036053, 8.48345542542134, 155.60037348272638,-0.08928493097479873 + 122.18343747755509,29.266038080882094, 157.45686992745814, 7.692771605286275, 157.86942469295408,-0.08928493097479873 + 124.45248868778279,29.017347349759465, 160.24161459455573, 6.920685052544364, 160.13847590318176,-0.08928493097479873 + 126.72153989801046,28.807923576182514, 163.12949795302734, 6.119802730935625, 162.38460740421522,-0.08928493097479873 + 128.99059110823814,28.549416105673462, 165.91424262012495, 5.381911153723081, 164.6765783236371,-0.08928493097479873 + 131.2596423184658,28.326903346247953, 168.9052646699705, 4.608024865426998, 166.9456295338648,-0.08928493097479873 + 133.5286935286935, 28.09784609389816, 171.69000933706812, 3.9277248259481183, 169.21468074409248,-0.08928493097479873 + 135.79774473892118,27.878605580934785, 174.7841700782877, 3.222228488710762, 171.48373195432015,-0.08928493097479873 + 138.06679594914885,27.646276082122853, 177.87833081950726, 2.567124746990352, 173.75278316454782,-0.08928493097479873 + 140.33584715937656, 27.4204910762352, 180.9724915607268, 1.9696125430036133, 176.02183437477552,-0.09910167036122175 + 142.60489836960423, 27.21106730265825, 184.16979099332036, 1.3685008679085158, 178.2908855850032,-0.08928493097479873 + 144.8739495798319, 26.98201005030846, 187.26395173453994, 0.8735735905098565, 180.55993679523087,-0.08928493097479873 + 147.14300079005957,26.779130769655787, 191.3158288956608, 0.2655200782772269, 182.82898800545857,-0.08928493097479873 + 149.41205200028728,26.569706996078832, 194.20002872944045, 0.06969170964415383, 185.09803921568624,-0.08928493097479873 + 151.57796451914095,26.381225599859572, 197.40592305298182,-0.02796060838803527, 187.36709042591391,-0.10891840974763767 + 154.1564318034906, 26.12795372368995, 199.1936603701309, -0.0652884569191059, 189.6361416361416,-0.11546290267192205 + 156.42548301371826, 25.94143567534797, -, -, 191.9051928463693,-0.08928493097479873 + 158.69453422394596, 25.72546740884674, -, -, 194.17424405659696,-0.05329021989126659 + 160.96358543417364,25.525860374656208, -, -, 196.44329526682463,-0.07292369866409842 + 163.2326366444013,25.336070079852096, -, -, -, - + 165.501687854629,25.136463045661564, -, -, -, - + 167.77073906485668, 24.95321724378173, -, -, -, - + 170.03979027508436, 24.77651593482618, -, -, -, - + 172.30884148531203,24.576908900635644, -, -, -, - + 174.57789269553973,24.396935345217948, -, -, -, - + 176.8469439057674,24.203872803951693, -, -, -, - + 179.0128564246211,24.030770966104498, -, -, -, - + 181.38504632622275, 23.86028692542701, -, -, -, - + 183.65409753645045,23.686857862933593, -, -, -, - + 185.92314874667812, 23.50361206105376, -, -, -, - + 188.1921999569058,23.349816477333185, -, -, -, - + 190.46125116713347, 23.17311516837763, -, -, -, - + 192.73030237736117,23.002958352346358, -, -, -, - + 194.99935358758884,22.836073782777223, -, -, -, - + 197.26840479781652, 22.66591696674595, -, -, -, - + 198.8154851684263,22.551388340571055, -, -, -, - diff --git a/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv b/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv new file mode 100644 index 00000000..57721214 --- /dev/null +++ b/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv @@ -0,0 +1,61 @@ +--- +# Digitized data from Figure 4 in Kaminski & Ribe 2004, https://doi.org/10.1111%2Fj.1365-246x.2004.02308.x +# This data records the mean angle of the a-axes of an olivine aggregate in the shear plane. +# The flow field is a static simple shear and the aggregate starts with random (i.e. isotropically distributed) grain orientations. + +schema: + delimiter: ',' + missing: '-' + fields: + - name: dimensionless_time_X0 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0 + type: float + fill: NaN + unit: degree + - name: dimensionless_time_X0d2 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0d2 + type: float + fill: NaN + unit: degree + - name: dimensionless_time_X0d4 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0d4 + type: float + fill: NaN + unit: degree +--- +dimensionless_time_X0, angle_X0,dimensionless_time_X0d2, angle_X0d2,dimensionless_time_X0d4, angle_X0d4 + 0,0.13513513513513775, 0,0.06756756756757198, -0.004098360655737765,3.552713678800501e-15 + 0.09836065573770492, 0.2702702702702737, 0.09836065573770492, 0.3378378378378404, 0.09631147540983603, 0.27027027027027195 + 0.19672131147540978,-0.2702702702702666, 0.19672131147540978,-0.2702702702702684, 0.194672131147541, -0.4054054054054035 + 0.3032786885245901, -2.297297297297293, 0.3012295081967214,-2.2972972972972965, 0.3012295081967213, -2.3648648648648614 + 0.40163934426229503, -9.189189189189186, 0.3995901639344262, -9.054054054054053, 0.39754098360655715, -8.108108108108098 + 0.49999999999999994,-20.945945945945944, 0.49999999999999994,-19.256756756756758, 0.4959016393442623, -16.486486486486484 + 0.598360655737705, -35, 0.598360655737705,-31.486486486486484, 0.5963114754098366, -26.081081081081045 + 0.6967213114754098, -42.70270270270271, 0.6987704918032793, -38.04054054054055, 0.6946721311475408, -31.351351351351347 + 0.7950819672131146, -45.00000000000001, 0.7971311475409832, -39.93243243243243, 0.7930327868852456, -33.04054054054057 + 0.8975409836065573, -45.4054054054054, 0.8954918032786883, -40.74324324324325, 0.8934426229508197, -33.44594594594596 + 1, -45.4054054054054, 1.0020491803278693, -41.14864864864866, 0.9979508196721318, -34.05405405405407 + 1.0983606557377048, -45.4054054054054, 1.100409836065574, -41.2837837837838, 1.0963114754098358, -34.662162162162176 + 1.1987704918032789, -45.20270270270271, 1.1987704918032789, -41.35135135135137, 1.196721311475411, -35.40540540540541 + 1.2991803278688523, -45.13513513513514, 1.2991803278688525, -41.35135135135137, 1.2950819672131149, -35.810810810810814 + 1.3975409836065575, -45.00000000000001, 1.3975409836065575, -41.55405405405406, 1.395491803278689, -36.418918918918926 + 1.4959016393442621, -45.00000000000001, 1.4959016393442617, -41.75675675675676, 1.4938524590163944, -36.6891891891892 + 1.6024590163934427, -45.00000000000001, 1.6024590163934427, -42.1621621621622, 1.5983606557377048, -36.89189189189188 + 1.7008196721311473, -45.00000000000001, 1.7008196721311477,-42.162162162162126, 1.6967213114754094, -36.891891891891895 + 1.7991803278688525, -45.00000000000001, 1.7991803278688525,-42.567567567567586, 1.797131147540984, -37.2972972972973 + 1.8975409836065575, -45.00000000000001, 1.897540983606557, -42.77027027027028, 1.8954918032786876, -37.50000000000001 + 1.9979508196721312, -45.2027027027027, 1.9979508196721303, -42.97297297297298, 1.9959016393442626, -37.905405405405396 + 2.102459016393443, -45.13513513513514, 2.096311475409836, -43.37837837837841, 2.0942622950819674, -38.31081081081081 + 2.1967213114754096, -45.20270270270271, 2.1946721311475406,-43.378378378378365, 2.192622950819672, -38.51351351351353 + 2.3012295081967222, -45.2027027027027, 2.3012295081967222, -43.37837837837841, 2.2971311475409837, -38.918918918918926 + 2.3995901639344264, -45.2027027027027, 2.3995901639344273, -43.58108108108105, 2.397540983606557, -39.12162162162162 + 2.500000000000001, -45.2027027027027, 2.500000000000001,-43.783783783783804, 2.495901639344262, -39.52702702702706 diff --git a/src/pydrex/diagnostics.py b/src/pydrex/diagnostics.py index 7331010a..458dcad5 100644 --- a/src/pydrex/diagnostics.py +++ b/src/pydrex/diagnostics.py @@ -240,7 +240,7 @@ def coaxial_index(orientations, axis1="b", axis2="a"): The BA index of [Mainprice et al. 2015](https://doi.org/10.1144/SP409.8) is derived from the eigenvalue `symmetry` diagnostics and measures point vs girdle - symmetry in the aggregate. $BA \in [0, 1]$ where $BA = 0$ corresponds to a perfect + symmetry in the aggregate. $BA ∈ [0, 1]$ where $BA = 0$ corresponds to a perfect axial girdle texture and $BA = 1$ represents a point symmetry texture assuming that the random component $R$ is negligible. May be less susceptible to random fluctuations compared to the raw eigenvalue diagnostics. diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index e098ff28..8276a8bd 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -202,6 +202,8 @@ class Mineral: - `n_grains` (int) — number of grains in the aggregate - `fractions` (list of arrays) — grain volume fractions for each texture snapshot - `orientations` (list of arrays) — grain orientation matrices for each texture snapshot + - `seed` (`int`) — seed used by the random number generator to set up the isotropic + initial condition when `fractions_init` or `orientations_init` are not provided """ @@ -214,6 +216,7 @@ class Mineral: orientations_init: np.ndarray = None fractions: list = field(default_factory=list) orientations: list = field(default_factory=list) + seed: int = None def __str__(self): # String output, used for str(self) and f"{self}", etc. @@ -249,7 +252,7 @@ def __post_init__(self): self.fractions_init = np.full(self.n_grains, 1.0 / self.n_grains) if self.orientations_init is None: self.orientations_init = Rotation.random( - self.n_grains, random_state=1 + self.n_grains, random_state=self.seed ).as_matrix() # Copy the initial values to the storage lists. @@ -358,8 +361,7 @@ def apply_gbs(orientations, fractions, config): len(np.nonzero(mask)) / len(fractions), ) # No rotation: carry over previous orientations. - # TODO: Should we really be resetting to initial orientations here? - orientations[mask, :, :] = self.orientations[0][mask, :, :] + orientations[mask, :, :] = self.orientations[-1][mask, :, :] fractions[mask] = config["gbs_threshold"] / self.n_grains fractions /= fractions.sum() _log.debug( diff --git a/src/pydrex/stats.py b/src/pydrex/stats.py index dd9edb55..6ac42545 100644 --- a/src/pydrex/stats.py +++ b/src/pydrex/stats.py @@ -2,60 +2,17 @@ import numpy as np from pydrex import geometry as _geo -from pydrex import minerals as _minerals -from pydrex import tensors as _tensors -_RNG = np.random.default_rng(seed=8845) - -def average_stiffness(minerals, config): - """Calculate average elastic tensor from a list of `minerals`. - - The `config` dictionary must contain volume fractions of all occurring mineral phases, - indexed by keys of the format `"_fraction"`. - - """ - n_grains = minerals[0].n_grains - assert np.all( - [m.n_grains == n_grains for m in minerals[1:]] - ), "cannot average minerals with varying grain counts" - elastic_tensors = [] - for phase in _minerals.MineralPhase: - if phase == _minerals.MineralPhase.olivine: - elastic_tensors.append( - _tensors.Voigt_to_elastic_tensor(_minerals.OLIVINE_STIFFNESS) - ) - elif phase == _minerals.MineralPhase.enstatite: - elastic_tensors.append( - _tensors.Voigt_to_elastic_tensor(_minerals.ENSTATITE_STIFFNESS) - ) - - average_tensor = np.zeros((3, 3, 3, 3)) - for n in n_grains: - for mineral in minerals: - if mineral.phase == _minerals.MineralPhase.olivine: - average_tensor += ( - _tensors.rotated_tensor(mineral.orientations[n, ...].transpose()) - * mineral.fractions[n] - * config["olivine_fraction"] - ) - elif mineral.phase == _minerals.MineralPhase.enstatite: - average_tensor += ( - _tensors.rotated_tensor(minerals.orientations[n, ...].transpose()) - * mineral.fractions[n] - * config["enstatite_fraction"] - ) - return _tensors.elastic_tensor_to_Voigt(average_tensor) - - -def resample_orientations(orientations, fractions, n_samples=None, rng=_RNG): +def resample_orientations(orientations, fractions, n_samples=None, seed=None): """Generate new samples from `orientations` weighed by the volume distribution. If the optional number of samples `n_samples` is not specified, it will be set to the number of original "grains" (length of `fractions`). - The argument `rng` can be used to specify a custom random number generator. + The argument `seed` can be used to seed the random number generator. """ + rng = np.random.default_rng(seed=seed) if n_samples is None: n_samples = len(fractions) sort_ascending = np.argsort(fractions) diff --git a/src/pydrex/utils.py b/src/pydrex/utils.py new file mode 100644 index 00000000..52b23a6e --- /dev/null +++ b/src/pydrex/utils.py @@ -0,0 +1,13 @@ +"""> PyDRex: Miscellaneous utility methods.""" +import numpy as np + + +def remove_nans(a): + """Remove NaN values from array.""" + a = np.asarray(a) + return a[~np.isnan(a)] + + +def angle_fse_simpleshear(strain): + """Get angle of FSE long axis anticlockwise from the X axis in simple shear.""" + return np.rad2deg(np.arctan(np.sqrt(strain**2 + 1) + strain)) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index b51bce1a..b67e97fd 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -1,6 +1,4 @@ """> PyDRex: Visualisation functions for test outputs and examples.""" -import functools as ft - import numpy as np from matplotlib import projections as mproj from matplotlib import pyplot as plt @@ -23,7 +21,13 @@ def polefigures( - orientations, ref_axes, i_range, density=False, savefile="polefigures.png", **kwargs + orientations, + ref_axes, + i_range, + density=False, + savefile="polefigures.png", + strains=None, + **kwargs, ): """Plot pole figures of a series of (Nx3x3) orientation matrix stacks. @@ -42,18 +46,32 @@ def polefigures( grid = fig.add_gridspec( 4, n_orientations, height_ratios=((1, 3, 3, 3)), hspace=0, wspace=0.2 ) - fig_time = fig.add_subfigure(grid[0, :]) + fig_strain = fig.add_subfigure(grid[0, :]) first_row = 1 - fig_time.suptitle( - f"N ⋅ (max strain) / {i_range.stop}", x=0.5, y=0.85, fontsize="small" - ) - ax_time = fig_time.add_subplot(111) - ax_time.set_frame_on(False) - ax_time.grid(False) - ax_time.yaxis.set_visible(False) - ax_time.xaxis.set_tick_params(labelsize="x-small", length=0) - ax_time.set_xticks(list(i_range)) - ax_time.set_xlim((-i_range.step / 2, i_range.stop - i_range.step / 2)) + ax_strain = fig_strain.add_subplot(111) + + if strains is None: + fig_strain.suptitle( + f"N ⋅ (max strain) / {i_range.stop}", x=0.5, y=0.85, fontsize="small" + ) + ax_strain.set_xlim( + (i_range.start - i_range.step / 2, i_range.stop - i_range.step / 2) + ) + ax_strain.set_xticks(list(i_range)) + else: + fig_strain.suptitle("strain (%)", x=0.5, y=0.85, fontsize="small") + ax_strain.set_xticks(strains[i_range.start : i_range.stop : i_range.step]) + ax_strain.set_xlim( + ( + strains[i_range.start] - strains[i_range.step] / 2, + strains[i_range.stop - i_range.step] + strains[i_range.step] / 2, + ) + ) + + ax_strain.set_frame_on(False) + ax_strain.grid(False) + ax_strain.yaxis.set_visible(False) + ax_strain.xaxis.set_tick_params(labelsize="x-small", length=0) fig100 = fig.add_subfigure( grid[first_row, :], edgecolor=plt.rcParams["grid.color"], linewidth=1 @@ -112,27 +130,6 @@ def polefigures( fig.savefig(_io.resolve_path(savefile)) -def check_marker_seq(func): - """Raises a `ValueError` if number of markers and data series don't match. - - The decorated function is expected to take the data as the first positional - argument, and a keyword argument called `markers`. - - """ - - @ft.wraps(func) - def wrapper(data, *args, **kwargs): - markers = kwargs["markers"] - if len(data) % len(markers) != 0: - raise ValueError( - "Number of data series must be divisible by number of markers." - + f" You've supplied {len(data)} data series and {len(markers)} markers." - ) - func(data, *args, **kwargs) - - return wrapper - - def _get_marker_and_label(data, seq_index, markers, labels=None): marker = markers[int(seq_index / (len(data) / len(markers)))] label = None @@ -141,81 +138,58 @@ def _get_marker_and_label(data, seq_index, markers, labels=None): return marker, label -@check_marker_seq def simple_shear_stationary_2d( + strains, + target_angles, angles, - indices, - timestop, + point100_symmetry, + angles_err=None, savefile="pydrex_simple_shear_stationary_2d.png", markers=("."), - labels=None, θ_fse=None, + labels=None, ): """Plot diagnostics for stationary A-type olivine 2D simple shear box tests.""" fig = plt.figure(figsize=(5, 8), dpi=300) grid = fig.add_gridspec(2, 1, hspace=0.05) ax_mean = fig.add_subplot(grid[0]) ax_mean.set_ylabel("Mean angle ∈ [0, 90]°") - ax_mean.axhline(0, color=plt.rcParams["axes.edgecolor"]) ax_mean.tick_params(labelbottom=False) - ax_strength = fig.add_subplot(grid[1], sharex=ax_mean) - ax_strength.set_ylabel("Texture strength (M-index)") - ax_strength.set_xlabel(r"Strain ($\dot{ε}_0 t$)") - - colors = [] - data_Kaminski2001 = _io.read_scsv( - _io.data("thirdparty") / "Kaminski2001_GBMshear.scsv" - ) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M0) / 200, - data_Kaminski2001.angle_M0, - alpha=0.33, - label=r"$M^{\ast}=0$", - ) - colors.append(lines[0].get_color()) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M50) / 200, - data_Kaminski2001.angle_M50, - alpha=0.33, - label=r"$M^{\ast}=50$", - ) - colors.append(lines[0].get_color()) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M200) / 200, - data_Kaminski2001.angle_M200, - alpha=0.33, - label=r"$M^{\ast}=200$", - ) - colors.append(lines[0].get_color()) - - for i, (misorient_angles, misorient_indices) in enumerate(zip(angles, indices)): - timestamps = np.linspace(0, timestop, len(misorient_angles)) + ax_mean.set_xlim((strains[0], strains[-1])) + ax_mean.set_ylim((0, 60)) + ax_symmetry = fig.add_subplot(grid[1], sharex=ax_mean) + ax_symmetry.set_xlim((strains[0], strains[-1])) + ax_symmetry.set_ylim((0, 1)) + ax_symmetry.set_ylabel(r"Texture symmetry ($P_{[100]}$)") + ax_symmetry.set_xlabel(r"Strain ($\dot{D}_0 t = γ/2$)") + + for i, (θ_target, θ, point100) in enumerate( + zip(target_angles, angles, point100_symmetry) + ): marker, label = _get_marker_and_label(angles, i, markers, labels) - ax_mean.plot( - timestamps, - misorient_angles, - marker, - markersize=5, - alpha=0.33, - label=label, - color=colors[i], - ) - ax_strength.plot( - timestamps, - misorient_indices, + lines = ax_mean.plot(strains, θ_target, alpha=0.66, label=label) + color = lines[0].get_color() + ax_mean.plot(strains, θ, marker, markersize=5, alpha=0.33, color=color) + if angles_err is not None: + ax_mean.fill_between( + strains, θ - angles_err[i], θ + angles_err[i], alpha=0.22, color=color + ) + ax_symmetry.plot( + strains, + point100, marker, markersize=5, alpha=0.33, + color=color, label=label, - color=colors[i], ) if θ_fse is not None: - ax_mean.plot(timestamps, θ_fse, "r--", label="FSE") - + ax_mean.plot(strains, θ_fse, linestyle=(0, (5, 5)), alpha=0.66, label="FSE") if labels is not None: ax_mean.legend() + ax_symmetry.legend() fig.savefig(_io.resolve_path(savefile)) @@ -228,7 +202,6 @@ def _lag_2d_corner_flow(θ): ) -@check_marker_seq def corner_flow_2d( x_paths, z_paths, diff --git a/tests/conftest.py b/tests/conftest.py index bd232377..7519e103 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,8 +2,8 @@ import matplotlib import pytest from _pytest.logging import LoggingPlugin, _LiveLoggingStreamHandler -from numpy import random as rn +from pydrex import io as _io from pydrex import logger as _log from pydrex import mock as _mock @@ -11,6 +11,22 @@ _log.quiet_aliens() # Stop imported modules from spamming the logs. +# Set up custom pytest CLI arguments. +def pytest_addoption(parser): + parser.addoption( + "--outdir", + metavar="DIR", + default=None, + help="output directory in which to store PyDRex figures/logs", + ) + parser.addoption( + "--runslow", + action="store_true", + default=False, + help="run slow tests (HPC cluster recommended, large memory requirement)", + ) + + # The default pytest logging plugin always creates its own handlers... class PytestConsoleLogger(LoggingPlugin): """Pytest plugin that allows linking up a custom console logger.""" @@ -34,10 +50,12 @@ def pytest_runtest_teardown(self, item): yield from self._runtest_for(item, "teardown") -# Hook up our logging plugin last, -# it relies on terminalreporter and capturemanager. @pytest.hookimpl(trylast=True) def pytest_configure(config): + config.addinivalue_line("markers", "slow: mark test as slow to run") + + # Hook up our logging plugin last, + # it relies on terminalreporter and capturemanager. if config.option.verbose > 0: terminal_reporter = config.pluginmanager.get_plugin("terminalreporter") capture_manager = config.pluginmanager.get_plugin("capturemanager") @@ -50,13 +68,13 @@ def pytest_configure(config): ) -def pytest_addoption(parser): - parser.addoption( - "--outdir", - metavar="DIR", - default=None, - help="output directory in which to store PyDRex figures/logs", - ) +def pytest_collection_modifyitems(config, items): + if config.getoption("--runslow"): + return # Don't skip slow tests. + skip_slow = pytest.mark.skip(reason="need --runslow option to run") + for item in items: + if "slow" in item.keywords: + item.add_marker(skip_slow) @pytest.fixture(scope="session") @@ -66,7 +84,7 @@ def outdir(request): @pytest.fixture(scope="function") def console_handler(request): - if request.config.option.verbose > 0: + if request.config.option.verbose > 0: # Show console logs if -v/--verbose given. return request.config.pluginmanager.get_plugin( "pytest-console-logger" ).log_cli_handler @@ -124,6 +142,12 @@ def ref_axes(request): @pytest.fixture -def rng(): - """A seeded RNG for tests to have (more) reproducible results.""" - return rn.default_rng(seed=8816) +def seeds(): + """1000 unique seeds for ensemble runs that need an RNG seed.""" + return _io.read_scsv(_io.data("rng") / "seeds.scsv").seeds + + +@pytest.fixture +def seed(): + """Default seed for test RNG.""" + return 8816 diff --git a/tests/test_corner_flow_2d.py b/tests/test_corner_flow_2d.py index 70d1b5ac..7f94d5c2 100644 --- a/tests/test_corner_flow_2d.py +++ b/tests/test_corner_flow_2d.py @@ -81,7 +81,7 @@ class TestOlivineA: def test_corner_prescribed_init_isotropic( self, params_Kaminski2001_fig5_shortdash, - rng, + seed, outdir, ): """Test CPO evolution in prescribed 2D corner flow. @@ -96,7 +96,7 @@ def test_corner_prescribed_init_isotropic( domain_height = 2.0e5 # Represents the depth of olivine-spinel transition. domain_width = 1.0e6 n_grains = 2000 - orientations_init = Rotation.random(n_grains, random_state=rng).as_matrix() + orientations_init = Rotation.random(n_grains, random_state=seed).as_matrix() n_timesteps = 20 # Number of places along the pathline to compute CPO. # Z-values at the end of each pathline. z_ends = list(map(lambda x: x * domain_height, (-0.1, -0.3, -0.54, -0.78))) @@ -191,7 +191,7 @@ def _get_velocity_gradient(point): # Loop over first dimension (time steps) of orientations. for idx, matrices in enumerate(mineral.orientations): orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx], rng=rng + matrices, mineral.fractions[idx], seed=seed ) direction_mean = _diagnostics.bingham_average( orientations_resampled, diff --git a/tests/test_geometry.py b/tests/test_geometry.py index e6c24cf6..24189396 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -26,7 +26,7 @@ def test_poles_example(hkl, ref_axes): np.testing.assert_allclose(ref_data["zvals"], zvals, atol=1e-16, rtol=0) -def test_lambert_equal_area(rng): +def test_lambert_equal_area(seed): """Test Lambert equal area projection.""" x, y = np.mgrid[-1:1:11j, -1:1:11j] x_flat = [j for i in x for j in i] @@ -35,6 +35,7 @@ def test_lambert_equal_area(rng): x_disk, y_disk = _geo.shirley_concentric_squaredisk(x_flat, y_flat) # Project onto the unit sphere by adding z = ± (1 - r). # Then project back onto the disk using Lambert equal-area, should be the same. + rng = np.random.default_rng(seed=seed) sign = rng.integers(low=0, high=2, size=len(x_disk)) x_laea, y_laea = _geo.lambert_equal_area( x_disk, y_disk, (-1) ** sign * (1 - (x_disk**2 + y_disk**2)) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 7aee7ef6..d82ecb89 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -1,356 +1,457 @@ -"""> PyDRex: 2D simple shear tests. - -NOTE: In scipy, rotations are represented as a matrix -that transforms [1, 0, 0] to the new a-axis vector. -We need the inverse of that rotation, -which represents the change of coordinates -from the grain-local to the global (Eulerian) frame. - -""" +"""> PyDRex: 2D simple shear tests.""" import contextlib as cl import numpy as np -from scipy.spatial.transform import Rotation +import pytest + +# from numpy import testing as nt +from scipy.interpolate import PchipInterpolator -from pydrex import core as _core from pydrex import diagnostics as _diagnostics +from pydrex import io as _io from pydrex import logger as _log from pydrex import minerals as _minerals from pydrex import stats as _stats +from pydrex import utils as _utils from pydrex import visualisation as _vis +# Subdirectory of `outdir` used to store outputs from these tests. +SUBDIR = "2d_simple_shear" + -class TestSinglePolycrystalOlivineA: - """Tests for a single A-type olivine polycrystal in 2D simple shear.""" +class TestOlivineA: + """Tests for A-type olivine polycrystals in 2D simple shear.""" + + class_id = "olivineA" def get_position(self, t): return np.zeros(3) - def test_shearYZ_initQ1( + @pytest.mark.slow + def test_dvdx_GBM_ensemble( self, params_Kaminski2001_fig5_solid, # GBM = 0 params_Kaminski2001_fig5_shortdash, # GBM = 50 params_Kaminski2001_fig5_longdash, # GBM = 200 - rng, + seeds, outdir, ): - """Test clockwise a-axis rotation around X. - - Initial condition: randomised a-axis orientation within the first - quadrant (between +Y and +Z axes) and steady flow with velocity gradient - - 0 0 0 - L = 0 0 2 - 0 0 0 + r"""Test a-axis alignment to shear in the Y direction. - Orientations set up for slip on (010)[100]. - Tests the effect of grain boundary migration, - similar to Fig. 5 in [Kaminski 2001]. + Velocity gradient: + $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ - [Kaminski 2001]: https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9 + Tests the effect of grain boundary migration, similar to Fig. 5 in + [Kaminski 2001](https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9). """ - strain_rate_scale = 2e-5 - timescale = 1 / (strain_rate_scale / 2) - n_grains = 1000 + strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. + timestamps = np.linspace(0, 2e5, 201) # Solve until D₀t=1 ('shear strain' γ=2). + n_timestamps = len(timestamps) # Number of steps + 1. + # i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. + # i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. def get_velocity_gradient(x): - Δv = np.zeros((3, 3)) - Δv[1, 2] = strain_rate_scale - return Δv - - # Initial orientations with a-axis in first quadrant of the YZ plane, - # and c-axis along the +X direction (important!) - # This means that the starting average is the same, - # which is convenient for comparing the texture evolution. - orientations_init = ( - Rotation.from_euler( - "zxz", - [[x * np.pi / 2, np.pi / 2, np.pi / 2] for x in rng.random(n_grains)], - ) - .inv() - .as_matrix() - ) - # Uncomment to check a-axis vectors, should be near [0, a, a]. - # assert False, f"{orientations_init[0:10, 0, :]}" - # Uncomment to check c-axis vectors, should be near [1, 0, 0]. - # assert False, f"{orientations_init[0:10, 2, :]}" - - # One mineral to test each value of grain boundary mobility. - minerals = [ - _minerals.Mineral( - _core.MineralPhase.olivine, - _core.MineralFabric.olivine_A, - _core.DeformationRegime.dislocation, - n_grains=n_grains, - fractions_init=np.full(n_grains, 1 / n_grains), - orientations_init=orientations_init, - ) - for i in range(3) - ] + # It is independent of time or position in this test. + grad_v = np.zeros((3, 3)) + grad_v[1, 0] = 2 * strain_rate + return grad_v + + shear_direction = [0, 1, 0] # Used to calculate the angular diagnostics. - # Optional plotting setup. + # Output setup with optional logging and data series labels. + θ_fse = np.empty((len(seeds), n_timestamps)) + angles = np.empty((3, len(seeds), n_timestamps)) + point100_symmetry = np.empty_like(angles) + optional_logging = cl.nullcontext() if outdir is not None: + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM_ensemble" + optional_logging = _log.logfile_enable(f"{out_basepath}.log") labels = [] - angles = [] - indices = [] - - for mineral, params in zip( - minerals, - ( - params_Kaminski2001_fig5_solid, # GBM = 0 - params_Kaminski2001_fig5_shortdash, # GBM = 50 - params_Kaminski2001_fig5_longdash, # GBM = 200 - ), - ): - time = 0 - timestep = 0.025 - timestop = 1 - deformation_gradient = np.eye(3) # Undeformed initial state. - while time < timestop * timescale: - deformation_gradient = mineral.update_orientations( - params, - deformation_gradient, - get_velocity_gradient, - pathline=(time, time + timestep * timescale, self.get_position), - ) - time += timestep * timescale - - n_timesteps = len(mineral.orientations) - misorient_angles = np.zeros(n_timesteps) - misorient_indices = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, [0, 1, 0] - ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled - ) - - # Check for uncorrupted record of initial condition. - assert np.isclose(misorient_angles[0], 45.0, atol=5.0) - assert misorient_indices[0] < 0.71 - # Check for mostly smoothly decreasing misalignment. - angles_diff = np.diff(misorient_angles) - assert np.max(angles_diff) < 3.6 - assert np.min(angles_diff) > -7.5 - assert np.sum(angles_diff) < -24.0 - - # Check alignment and texture strength (half way and final value). - halfway = round(n_timesteps / 2) - match params["gbm_mobility"]: - case 0: - np.testing.assert_allclose( - misorient_angles, - misorient_angles[0] - * np.exp( - np.linspace(0, timestop, n_timesteps) - * (np.cos(2 * np.deg2rad(misorient_angles[0])) - 1) - ), - atol=5.2, - rtol=0.1, + + with optional_logging: + for p, params in enumerate( + ( + params_Kaminski2001_fig5_solid, # GBM = 0 + params_Kaminski2001_fig5_shortdash, # GBM = 50 + params_Kaminski2001_fig5_longdash, # GBM = 200 + ), + ): + for s, seed in enumerate(seeds): + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], seed=seed + ) + deformation_gradient = np.eye(3) # Undeformed initial state. + θ_fse[s, 0] = 45 + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "$M^∗=%s$; #%s/%s; step %s/%s (t=%s)", + params["gbm_mobility"], + s + 1, + len(seeds), + t, + n_timestamps - 1, + time, + ) + deformation_gradient = mineral.update_orientations( + params, + deformation_gradient, + get_velocity_gradient, + pathline=(time, timestamps[t], self.get_position), + ) + _, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info( + "› velocity gradient = %s\n› strain D₀t=%s", + get_velocity_gradient(None).flatten(), + strain_rate * time, + ) + if p == 0: + θ_fse[s, t] = _diagnostics.smallest_angle( + fse_v, shear_direction + ) + + texture_symmetry = np.zeros(n_timestamps) + # mean_angles = np.zeros(n_timestamps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx], seed=seed + ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # mean_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, shear_direction + # ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] + + # Use SCCS axis (hexagonal symmetry) for the angle instead (opt). + mean_angles = np.array( + [ + _diagnostics.smallest_angle( + _diagnostics.anisotropy(v)[1][2, :], shear_direction + ) + for v in _minerals.voigt_averages([mineral], params) + ] ) - assert np.isclose(misorient_angles[halfway], 25, atol=2.0) - assert np.isclose(misorient_angles[-1], 17.0, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.925, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.975, atol=0.05) - case 50: - assert np.isclose(misorient_angles[halfway], 15, atol=2.0) - assert np.isclose(misorient_angles[-1], 10, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.925, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.97, atol=0.03) - case 200: - assert np.isclose(misorient_angles[halfway], 9, atol=1.5) - assert np.isclose(misorient_angles[-1], 7, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.975, atol=0.05) - assert np.isclose(misorient_indices[-1], 0.99, atol=0.03) - - # Optionally store plotting metadata. - if outdir is not None: - labels.append(f"$M^∗$ = {params['gbm_mobility']}") - angles.append(misorient_angles) - indices.append(misorient_indices) - - # Optionally plot figure. + + # Store data series and optional plotting metadata. + angles[p, s, :] = mean_angles + point100_symmetry[p, s, :] = texture_symmetry + + # Update labels and store the last mineral of the ensemble for polefigs. + if outdir is not None: + labels.append(f"$M^∗$ = {params['gbm_mobility']}") + mineral.save( + f"{out_basepath}.npz", + postfix=f"M{params['gbm_mobility']}", + ) + + # Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`. + _log.info("interpolating target CPO angles...") + strains = timestamps * strain_rate + data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv") + cs_M0 = PchipInterpolator( + _utils.remove_nans(data.equivalent_strain_M0) / 200, + _utils.remove_nans(data.angle_M0), + ) + cs_M50 = PchipInterpolator( + _utils.remove_nans(data.equivalent_strain_M50) / 200, + _utils.remove_nans(data.angle_M50), + ) + cs_M200 = PchipInterpolator( + _utils.remove_nans(data.equivalent_strain_M200) / 200, + _utils.remove_nans(data.angle_M200), + ) + target_angles = [cs_M0(strains), cs_M50(strains), cs_M200(strains)] + + # Take ensemble means and optionally plot figure. + result_angles = angles.mean(axis=1) + result_angles_err = angles.std(axis=1) + result_point100_symmetry = point100_symmetry.mean(axis=1) + result_θ_fse = θ_fse.mean(axis=0) if outdir is not None: + schema = { + "delimiter": ",", + "missing": "-", + "fields": [ + { + "name": "strain", + "type": "integer", + "unit": "percent", + "fill": 999999, + } + ], + } + _io.save_scsv( + f"{out_basepath}_strains.scsv", + schema, + [[int(γ * 200) for γ in strains]], + ) _vis.simple_shear_stationary_2d( - angles, - indices, - timestop=timestop, - savefile=f"{outdir}/simple_shearYZ_stationary_olivineA_initQ1.png", + strains, + target_angles, + result_angles, + result_point100_symmetry, + angles_err=result_angles_err, + savefile=f"{out_basepath}.png", markers=("o", "v", "s"), + θ_fse=result_θ_fse, labels=labels, ) - def test_shearXZ_initQ1( + # # Check that FSE is correct. + # # First, get theoretical FSE axis for simple shear. + # # We want the angle from the Y axis (shear direction), so subtract from 90. + # θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + # nt.assert_allclose(result_θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # # Check Bingham angles, ignoring the first portion. + # # Average orientations of near-isotropic distributions are unstable. + # nt.assert_allclose( + # result_θ_fse[i_first_cpo:], + # result_angles[0][i_first_cpo:], + # rtol=0.11, + # atol=0, + # ) + # nt.assert_allclose( + # target_angles[0][i_first_cpo:], + # result_angles[0][i_first_cpo:], + # rtol=0.11, + # atol=0, + # ) + # nt.assert_allclose( + # target_angles[1][i_first_cpo:], + # result_angles[1][i_first_cpo:], + # rtol=0, + # atol=5.7, + # ) + # nt.assert_allclose( + # target_angles[2][i_first_cpo:], + # result_angles[2][i_first_cpo:], + # rtol=0, + # atol=5.5, + # ) + + # # Check point symmetry of [100] at strains of 0%, 50%, 100%, 150% & 200%. + # nt.assert_allclose( + # [0.015, 0.11, 0.19, 0.27, 0.34], + # result_point100_symmetry[0].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.15, 0.33, 0.57, 0.72], + # result_point100_symmetry[1].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.22, 0.64, 0.86, 0.91], + # result_point100_symmetry[2].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + + @pytest.mark.slow + def test_dudz_GBS_ensemble( self, - params_Kaminski2004_fig4_triangles, # GBS = 0.4 - params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_circles, # GBS = 0 - rng, + params_Kaminski2004_fig4_squares, # GBS = 0.2 + params_Kaminski2004_fig4_triangles, # GBS = 0.4 + seeds, outdir, ): - """Test clockwise a-axis rotation around Y. - - Initial condition: randomised a-axis orientation within the first - quadrant (between +X and +Z axes) and steady flow with velocity gradient - - 0 0 2 - L = 0 0 0 - 0 0 0 + r"""Test a-axis alignment to shear in X direction. - Orientations set up for slip on (010)[100]. + Velocity gradient: + $$\bm{L} = \begin{bmatrix} 0 & 0 & 2 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ """ - strain_rate_scale = 2e-5 - timescale = 1 / (strain_rate_scale / 2) - n_grains = 1000 + strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. + timestamps = np.linspace(0, 5e5, 251) # Solve until D₀t=2.5 ('shear' γ=5). + n_timestamps = len(timestamps) + # i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. def get_velocity_gradient(x): - Δv = np.zeros((3, 3)) - Δv[0, 2] = strain_rate_scale - return Δv - - # Initial orientations with a-axis in first quadrant of the XZ plane, - # and c-axis along the -Y direction (important!) - # This means that the starting average is the same, - # which is convenient for comparing the texture evolution. - orientations_init = ( - Rotation.from_euler( - "zxz", - [[x * np.pi / 2, np.pi / 2, 0] for x in rng.random(n_grains)], - ) - .inv() - .as_matrix() - ) - # Uncomment to check a-axis vectors, should be near [a, 0, a]. - # assert False, f"{orientations_init[0:10, 0, :]}" - # Uncomment to check c-axis vectors, should be near [0, -1, 0]. - # assert False, f"{orientations_init[0:10, 2, :]}" - - # One mineral to test each grain boundary sliding threshold. - minerals = [ - _minerals.Mineral( - _core.MineralPhase.olivine, - _core.MineralFabric.olivine_A, - _core.DeformationRegime.dislocation, - n_grains=n_grains, - fractions_init=np.full(n_grains, 1 / n_grains), - orientations_init=orientations_init, - ) - for i in range(3) - ] + # It is independent of time or position in this test. + grad_v = np.zeros((3, 3)) + grad_v[0, 2] = 2 * strain_rate + return grad_v + + shear_direction = [1, 0, 0] # Used to calculate the angular diagnostics. - # Optional plotting and logging setup. + # Output setup with optional logging and data series labels. + θ_fse = np.empty((len(seeds), n_timestamps)) + angles = np.empty((3, len(seeds), n_timestamps)) + point100_symmetry = np.empty_like(angles) optional_logging = cl.nullcontext() if outdir is not None: - optional_logging = _log.logfile_enable( - f"{outdir}/simple_shearXZ_initQ1.log" - ) + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS_ensemble" + optional_logging = _log.logfile_enable(f"{out_basepath}.log") labels = [] - angles = [] - indices = [] with optional_logging: - for mineral, params in zip( - minerals, + for p, params in enumerate( ( - params_Kaminski2004_fig4_triangles, # GBS = 0.4 - params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_circles, # GBS = 0 + params_Kaminski2004_fig4_squares, # GBS = 0.2 + params_Kaminski2004_fig4_triangles, # GBS = 0.4 ), ): - time = 0 - timestep = 0.025 - timestop = 1 - deformation_gradient = np.eye(3) # Undeformed initial state. - while time < timestop * timescale: - deformation_gradient = mineral.update_orientations( - params, - deformation_gradient, - get_velocity_gradient, - pathline=(time, time + timestep * timescale, self.get_position), - # atol=1e-99, - ) - time += timestep * timescale - - n_timesteps = len(mineral.orientations) - misorient_angles = np.zeros(n_timesteps) - misorient_indices = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + for s, seed in enumerate(seeds): + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], seed=seed ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, [1, 0, 0] - ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled + deformation_gradient = np.eye(3) # Undeformed initial state. + θ_fse[s, 0] = 45 + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "X = %s; # %s/%s; step %s/%s (t = %s)", + params["gbs_threshold"], + s + 1, + len(seeds), + t, + n_timestamps - 1, + time, + ) + deformation_gradient = mineral.update_orientations( + params, + deformation_gradient, + get_velocity_gradient, + pathline=(time, timestamps[t], self.get_position), + ) + _, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info( + "› velocity gradient = %s", + get_velocity_gradient(None).flatten(), + ) + _log.info("› strain D₀t = %s", strain_rate * time) + if p == 0: + θ_fse[s, t] = _diagnostics.smallest_angle( + fse_v, shear_direction + ) + + texture_symmetry = np.zeros(n_timestamps) + # mean_angles = np.zeros(n_timestamps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx], seed=seed + ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # mean_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, shear_direction + # ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] + + # Use SCCS axis (hexagonal symmetry) for the angle instead (opt). + mean_angles = np.array( + [ + _diagnostics.smallest_angle( + _diagnostics.anisotropy(v)[1][2, :], shear_direction + ) + for v in _minerals.voigt_averages([mineral], params) + ] ) - # Check for uncorrupted record of initial condition. - assert np.isclose(misorient_angles[0], 45.0, atol=5.0) - assert misorient_indices[0] < 0.71 - # Check for mostly smoothly decreasing misalignment. - angles_diff = np.diff(misorient_angles) - assert np.max(angles_diff) < 3.6 - assert np.min(angles_diff) > -7.5 - assert np.sum(angles_diff) < -24.0 - # Check that recrystallization is causing faster alignment. - np.testing.assert_array_less( - misorient_angles - 3.8, # Tolerance for GBM onset latency. - misorient_angles[0] - * np.exp( - np.linspace(0, timestop, n_timesteps) - * (np.cos(2 * np.deg2rad(misorient_angles[0])) - 1) - ), - ) - - # Check alignment and texture strength (half way and final value). - halfway = round(n_timesteps / 2) - match params["gbs_threshold"]: - case 0: - assert np.isclose(misorient_angles[halfway], 11, atol=1.5) - assert np.isclose(misorient_angles[-1], 8, atol=1.25) - assert np.isclose(misorient_indices[halfway], 0.975, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.99, atol=0.03) - case 0.2: - assert np.isclose(misorient_angles[halfway], 13, atol=2.075) - assert np.isclose(misorient_angles[-1], 11, atol=1.5) - assert np.isclose(misorient_indices[halfway], 0.755, atol=0.08) - assert np.isclose(misorient_indices[-1], 0.755, atol=0.075) - case 0.4: - assert np.isclose(misorient_angles[halfway], 19, atol=2.0) - assert np.isclose(misorient_angles[-1], 16, atol=2.5) - assert misorient_indices[halfway] < 0.75 - assert misorient_indices[-1] < 0.71 - - # Optionally store plotting metadata. + # Store data series and optional plotting metadata. + angles[p, s, :] = mean_angles + point100_symmetry[p, s, :] = texture_symmetry + + # Update labels and store the last mineral of the ensemble for polefigs. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") - angles.append(misorient_angles) - indices.append(misorient_indices) + mineral.save( + f"{out_basepath}.npz", + postfix=f"X{params['gbs_threshold']}", + ) - # Optionally plot figure. + # Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`. + _log.info("interpolating target CPO angles...") + strains = timestamps * strain_rate + data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv") + cs_X0 = PchipInterpolator( + _utils.remove_nans(data.dimensionless_time_X0), + 45 + _utils.remove_nans(data.angle_X0), + ) + cs_X0d2 = PchipInterpolator( + _utils.remove_nans(data.dimensionless_time_X0d2), + 45 + _utils.remove_nans(data.angle_X0d2), + ) + cs_X0d4 = PchipInterpolator( + _utils.remove_nans(data.dimensionless_time_X0d4), + 45 + _utils.remove_nans(data.angle_X0d4), + ) + target_angles = [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)] + + # Take ensemble means and optionally plot figure. + result_angles = angles.mean(axis=1) + result_angles_err = angles.std(axis=1) + result_point100_symmetry = point100_symmetry.mean(axis=1) + result_θ_fse = θ_fse.mean(axis=0) if outdir is not None: + schema = { + "delimiter": ",", + "missing": "-", + "fields": [ + { + "name": "strain", + "type": "integer", + "unit": "percent", + "fill": 999999, + } + ], + } + _io.save_scsv( + f"{out_basepath}_strains.scsv", + schema, + [[int(γ * 200) for γ in strains]], + ) _vis.simple_shear_stationary_2d( - angles, - indices, - timestop=timestop, - savefile=f"{outdir}/simple_shearXZ_stationary_olivineA_initQ1.png", + strains, + target_angles, + result_angles, + result_point100_symmetry, + angles_err=result_angles_err, + savefile=f"{out_basepath}.png", markers=("o", "v", "s"), + θ_fse=result_θ_fse, labels=labels, ) + + # # Check that FSE is correct. + # # First, get theoretical FSE axis for simple shear. + # # We want the angle from the Y axis (shear direction), so subtract from 90. + # θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + # nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # # Check point symmetry of [100] at strains of 0%, 100%, 200%, 300% & 400%. + # nt.assert_allclose( + # [0.015, 0.52, 0.86, 0.93, 0.94], + # point100_symmetry[0].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.42, 0.71, 0.77, 0.79], + # point100_symmetry[1].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.36, 0.57, 0.6, 0.62], + # point100_symmetry[2].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # )