From 788e59c4ef4746502b7d451e2c9294c0a5b18ef2 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 31 Jul 2023 16:52:31 +0200 Subject: [PATCH] Add docu and validation files --- doc/ADI_scheme.org | 1 - doc/ValidationHetDiff.org | 114 ++++++++++++ doc/images/deSolve_AlphaHet1.png | Bin 0 -> 19543 bytes examples/CMakeLists.txt | 5 +- examples/FTCS_2D_proto_closed_mdl.cpp | 77 ++++++++ examples/FTCS_2D_proto_example_mdl.cpp | 77 ++++++++ proto/FTCS.ipynb | 2 +- scripts/Adi2D_Reference.R | 2 +- scripts/HetDiff.R | 237 +++++++++++++++++++++++++ scripts/gold/HetDiff1.csv | 121 +++++++++++++ 10 files changed, 631 insertions(+), 5 deletions(-) create mode 100644 doc/ValidationHetDiff.org create mode 100644 doc/images/deSolve_AlphaHet1.png create mode 100644 examples/FTCS_2D_proto_closed_mdl.cpp create mode 100644 examples/FTCS_2D_proto_example_mdl.cpp create mode 100644 scripts/HetDiff.R create mode 100644 scripts/gold/HetDiff1.csv diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 1124a92..2de70da 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -536,4 +536,3 @@ left derivative becomes zero and we are left with: &\displaystyle =\frac{\alpha_{i-1/2}}{\Delta x^2} (C_{i-1} - C_i) \end{aligned} \end{equation} - diff --git a/doc/ValidationHetDiff.org b/doc/ValidationHetDiff.org new file mode 100644 index 0000000..c9b5613 --- /dev/null +++ b/doc/ValidationHetDiff.org @@ -0,0 +1,114 @@ +#+TITLE: 2D Validation Examples +#+AUTHOR: MDL +#+DATE: 2023-07-31 +#+STARTUP: inlineimages +#+LATEX_CLASS_OPTIONS: [a4paper,9pt] +#+LATEX_HEADER: \usepackage{fullpage} +#+LATEX_HEADER: \usepackage{amsmath, systeme} +#+LATEX_HEADER: \usepackage{graphicx} +#+LATEX_HEADER: \usepackage{} +#+OPTIONS: toc:nil + + +* Simple setup using deSolve/ReacTran + +- Square of side 10 +- Discretization: 11 \times 11 cells +- All boundaries closed +- Initial state: 0 everywhere, 1 in the center (6,6) +- Time step: 1 s, 10 iterations + +The matrix of spatially variable diffusion coefficients \alpha is +constant in 4 quadrants: + +\alpha_{x,y} = + +| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | +| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | +| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | +| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | +| 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | +| 1 | 1 | 1 | 1 | 1 | 1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | + +The relevant part of the R script used to produce these results is +presented in listing 1; the whole script is at [[file:scripts/HetDiff.R]]. +A visualization of the output of the reference simulation is given in +figure [[#fig:1][1]]. + +Note: all results from this script are stored in the =outc= matrix by +the =deSolve= function. I stored a different version into +[[file:../scripts/gold/HetDiff1.csv]]: this file contains 11 columns (one +for each time step including initial conditions) and 121 rows, one for +each domain element, with no headers. + +#+caption: Result of ReacTran/deSolve solution of the above problem at 4 +[[./images/deSolve_AlphaHet1.png]] + + +#+name: lst:1 +#+begin_src R :language R :frame single :caption Listing 1, generate reference simulation using R packages deSolve/ReacTran :captionpos b :label lst:1 +library(ReacTran) +library(deSolve) + +## harmonic mean +harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & y have different lengths") + 2/(1/x+1/y) +} + +N <- 11 # number of grid cells +ini <- 1 # initial value at x=0 +N2 <- ceiling(N/2) +L <- 10 # domain side + +## Define diff.coeff per cell, in 4 quadrants +alphas <- matrix(0, N, N) +alphas[1:N2, 1:N2] <- 1 +alphas[1:N2, seq(N2+1,N)] <- 0.1 +alphas[seq(N2+1,N), 1:N2] <- 0.01 +alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001 + +cmpharm <- function(x) { + y <- c(0, x, 0) + ret <- numeric(length(x)+1) + for (i in seq(2, length(y))) { + ret[i-1] <- harm(y[i], y[i-1]) + } + ret +} + +## Construction of the 2D grid +x.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +y.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +grid2D <- setup.grid.2D(x.grid, y.grid) +dx <- dy <- L/N + +D.grid <- list() +## Diffusion coefs on x-interfaces +D.grid$x.int <- apply(alphas, 1, cmpharm) +## Diffusion coefs on y-interfaces +D.grid$y.int <- t(apply(alphas, 2, cmpharm)) + +# The model +Diff2Dc <- function(t, y, parms) { + CONC <- matrix(nrow = N, ncol = N, data = y) + dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + return(list(dCONC)) +} + +## initial condition: 0 everywhere, except in central point +y <- matrix(nrow = N, ncol = N, data = 0) +y[N2, N2] <- ini # initial concentration in the central point... + +## solve for 10 time units +times <- 0:10 +outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL, + dim = c(N, N), lrw = 1860000) +#+end_src + diff --git a/doc/images/deSolve_AlphaHet1.png b/doc/images/deSolve_AlphaHet1.png new file mode 100644 index 0000000000000000000000000000000000000000..796e2187ebbd8d921e0d4af8ebcfee5a51b7c77c GIT binary patch literal 19543 zcmdUXc|6p6`}b&_lSCUO%b}tyaZq;Cq0%BM6=fMBDzay1R1>F!BMDh13CWg}osp6y zONEeS$THcNvCWv}xxUkKx1Rev&+mEuczWHZ`&DG-`~6;@>w0h3Ts@?#&9#hw8488s z+P8Pt5fo}+De@nB0sPH{#{o6)f4@2GJ>iT(EkB6-$ML2vT^EJ=9kp-QpGV!_4R<&? zwDA}*ubC5?#7`s|>$vaI3f3My;ktyO7l6dz?8}+e4mb?2(Jofw!4t6}Jk7@_Xuk4JH{Vnlr;O#& z?c28>J9f;Tm}u9QGs4Xp=N>DWA1xd!G4kNB(#6$zPQBKaDI9);qb1{xm|G>OZ8r07 z%eJl;rau*?_h6V7=x)Xsj&(Z1h(L=%PpwhBAKK8aHE2pw?LT z*M{juw#vy~KCQPeB_pGbIn(Sr(yWQzxN#$nHH#D0oO>>=5n5$wS2Xd0r);&J+8=6a z0zrFlBh9|jVO8y(L+Mm1Ra`I5d#*cxTpK+*n3@?ZND&-!*G6|IMws%iUth)Cx^CUN zL+AXQI|&5u$trCjk4<4!mMrt}DC@Yx)>dMyfC?$s*Sqhvih*hzUNn0x4vvKfrnbXSavmy zqOYeu-Keh@#F**y8e@cEsA1Vor$l4=ge{=9W?OkQRbqwKt!p9TaAIMJj;cCJ)axh#fO~acxgRz4~l>oVu5(&vX-(8jht#6RFW~La?l9Ebc>mfT{OH z`D$@p`iwnm2HrbXNZHBq_1)cY)7x{-8>z%HV@cH54faKj%^4;^K|!jn!|Q~EN(Pdm zUj$QMLPbIiV`#5KR8&8mN=V$iy?)CL z2I{GM6(QVQT+)UQcfwITid(|Lk3#jY@w&r- z+Ud6hg+l(g=nM+Q^Z75n`NP*={qb8reYNfylgWJh_N`->y1IIZ+GHg+Zcb~O1AgKh z(bQ|?i5`yDK%&;U>Bj0LH*DBpYnW?iA{W(OL@ec^avyYvVT?*|wpAvz?=Hp073m*5tR|$Fc%;=)BK*OF2e;L)G!ZUF zwGp(?0o=kGL!@kJ0TfxuCHVUFv=iA@W1Yn!cxaS2ZrqrhbUL;mGHcNpzq1$S#=B!x z-Ix>%%h|xd03Dx|HIdV1S2$WQ^k!fx!=$ieGCb3^_0jHdl+~lkHmYL9SsL%gm8@Op zqJzm)@v^FSUF%1(3g`Jc%HTZv`};3jwoDUUTwDzGVG$i~{S+%glb3ayr=0DjwxH_H zFR4u;^IVGEfXcpv;*w2?T!s2jzq)1N&t`-3%~S`J&bd(bw-{dht`f|4IIa>Yz1Lp7 zVJ9xGa@h}q1^%Mua7G!2?AF}LC9eqLj~0CQ8l|@HH);Cq%MVxdW2%L&{kp(D*HmPJ z@s)*UN9_8z?qE1!Y>o5!2LwdlzCB5y7-X8796BVsgVfN_@a7E^HBuhK*Wr$%c=-AG zZ%LgvVq#)~#bS0u=_V`Mwmh`!D6o{v_n`Filuh(hEQXh>a}yL4?C9vYy6Z0dA`*$| z>C&pi^KCha2O{<)ynbTS4FBfBnKK_YkmmehdX`c%ii_2SLaZy`Up#qo!M5blBRTxv z5fOZW<{guAzG+a8JVlO^8N7H^4g{^-mH50>eQE<%}Q6sX{NG_5W<4k+Uf>FYN% zHH8bhxw>W;=gZ+`<>Ve&)!r`lnzfR5+qZ8Yj3bBjswmVJNy*gI)V2~IF9rVe^mLOV zcSj2g3uk9%3yYp$K~>k0cTb$jb<##TF-nf7rU#osb;RVKk@PlAz9m_1dQ=s_rE^#; zd*Hf?&dhZh3=OyJT&VM4?{BCn@^5e-eu6?iR$8YhpAI#CU(*65 zEG_$odazn8zH))iwjF)iDpz^BiU#VBe{M&gs|D##X2DKtSNK&)-=5saEx}t#eM3Or z(>K^hHd^T_kJ(8sZf?cx+xy7mOw;0mJ9o@%nU$d%y=I2DZryq!->I+2}5$jlUuR&yWkg5H`5$237Z#CPBT_K3URz*s$k1R7tze%;X# zyLXd>gamDVzU|$+M-(VUA0Hn|f9-)yt5S~ck4IU_87ZNX^t`;h)P(ATS6|&_V>9&V z^h-_j%Lo}^=<`rA@Xwi3daR?!H%q=kH6&+dWs>)JMHF!?~cCL)CpsEt%j;^WR{X z*AuLIa}Bi4$k=O1Q#R#^BKW_f?%cnBrhVevxpM^t1&?|FL4`{i`19=9v*%q~8=sUR zk7gdFP8I7tvHQwWB}seZn?|NK)d%wyz9Gq04PERIP(K{5D?xx6NwskR#`A<1ap zzN!=)m!*C3s;j%;)~#7;)}Ouf)8wk;g1aqdkhBOzqqzT2H<1UC)0r>28TE=#k^Oy} zV%N^c0eZ`AjvX}P7(PaRLPA1e zp=#(}JkQ4pyCrSJtJ#usz6~&I!s`$3p4cGMAJDI5a&s$sCWe(l_dd?a z39dCbe%zmDa&q!iR!fP`yr-upa#ga5>I^5S9??)o%EE;UckkYvTsJj6y_T2Psv^r! zDoRIp2kZH2J>>10*`1t7MJmS3*mH%*^PY*L65)c$M^a!8=Efk9DaX2A2u`LS#pK|@ zgD}}#+gn@VTvhXwRaREs-flB9+Oa_sAmo;)<$$Rc_lqv}t3pC(S*bvu zga9hcS<)U}V*&CS4ol$KIJ_;lVjWo1$CNIT%niUgT&N>G<#WQ`tM}bpA>&XZY~aiB zO<~t5Bk%Nrt?$#D!qLaK`BCm+jJ_M;bu!PMJquBvG1S-BS5fyhG8GG>KIQg(=rwD1 zYymV=t3-?VGuMbozBK?VHXpck)WD#*xtXBCq77w|yeYMPeSMM!8Am<-sqqcC&^8KF zEBQ*sudv32nU^IvY8;?RPMx*452@&BbrcOfbN+m7aCCGuw}2ARAG+Bye=ds@e!r-y z3wJm=Gpq?3D>eVxCEC-!bjgxy*O&yH`tZ;Yz}A{zi!4Lv48hRM2BwFaaWmNjpZ+Lo z4SF{>kv5q6`pi*t4l6OZfA-DCB=3sgYZmpqaYl;>g+P%GD}8E-aE_qbMDVE@1u)sB zN-r*aunZmo&sPm_n0g6V6I^I}p_{!c6R0db5lh?r1n`?}H!AWdc1!f6MTkz{pT8#M%F&#^|N|#$B1&UDmzc=InT9eS8^;>LqeAfCOit}Dw z@9`{Dj0~gom#FQj&TP4ZTL$+1==kmn{nEukyj#m;-Ukp~PD`$#&t%02vFgSXjfq!e zTsm&G@y)TEn93pEJIFGjLm}Oxd~(V-?lBLE3Yl?k8l_cEO|8AL@g{*l^qd~Nbm>wS zH7+jh)~#D_-n==Mdi)BHZk*a#_Sw_YdUj*;fd{FnLqkIqVhsaR{#+XX;Ya9*DcW}m zXu5Y@T#7J{e+}D`qU~1w_twWNeR31K^%c`y?9n1}PrJIh^0a5Z?s{=no%`^pF#hF# z>_$ZvGpKLY0?p9rCt-{fghE2jx9v?lf{g8GpaK>Y7`0DN52H|bl9G~AQzcbLhlh>x zolc|8EiEl;Z*TXS9pw^PzI-`gq?-{DOOeNP==kyD1_rr4^D}1@P8t~%x=*}+@!~~) ze?R;VaBUd6Ik_&@D983Uk#?t-(ujtGZ?h+=|EHTh^kvYZ?MpxaC|$@sH8ZmfTcTFD zo(2W925!Lbc%J?H_rG`nm~+@wqBVBI2^c{$NUK@@Z~<_DX00euU;#b9g7{0cLVx~s zgdO@RS097f+IPGYX6vqRTF_^j_z{!s6Xe`|4Rr?-Vu3vD&tGHzh}P%(@HH-O0WuN` z$omOK9;{O{Wcu1SCPEhg4>-TCzsyt{ z`C%e|)kreZ18foBMN4QIhKbWu^CNX5)O z9OD7{dSSR;TNL^5ZS2!p1{S69PGY{7QqJx2ak@9?h-u=xRKjACSo2KA{IRuWFCJZ4*fnx3BToXEip&n}w~dYMZX};pM{8>` zG%L1hgk=mGXz$shi3Ze?+u_=_Il}nv14K@YPfJV7$gl%!`T7h13)hMF<)Bv7*4D!B ztxMPN@`|mx$|+NHC>;T(t>S#S*NPc4ZTcyVA3l62yTYl39%{}cQv0JTO5VMBb9Wd; zzK1Y9c>{oRElR7du359@{%@sHz!bmIACPJMkq$_x^XPwhxa>KJ+6gj=-MMqs_4U|d z&*>Dtbt3}_pR<)Jz7y(&v5 z!094*uhIMp=#SQI&e=>swcB%NGO%~0ruwGlhFFwbCkAilG(pvgv1*LF(0P8rbbu49 zp)q-L&9cvwQWXD+{or$eyYrL*GP||SggfAuNm~dB9 z)!noT=$ZxPKPrYRhy!aSkHf)?j0xZ9M@hiUe%F_f4(@!ktw97|x#;rc%RNHhJa%1q z8KkvoOB=1!tkgMl^r*t~OO>2eWhBh}dyL0c}y`m!eht~-( zO;8RcAQAmLE-8TVr2{VDLvC&RPbZoP9>mc*xw&rC$-dUwR~Y6*$l^18Rf75p{3!P^ z5jSsM584=7*)npp`sGVKpGb8!^7*c;k>& z{wPX7|NA8!t%JVyHQLWN=UfWqfazs!cA{tI4)DC3R-GFRB2*1b5Tc1v&GsU^OPm?2 zM%<)oCnVKk=7rx{oGfc+ce^^jwQ!{;lox;w@=q^0df`V3!ER-2sicZa;^&|exk(?T z5+}A)qG$H8#rj2^b9bN1st~d{cW!@jLvkGpVPh~-aZ?GxTUP?JDP=W%%R2agFxX({ z-vODo&+ofs?=xidV3BMXuozh%6%`d=o4{4k@ocf=Jm}LR_{&SVeWz2gTp~3D3^6~S z+9S-$JIPUF2#4%jerx0E4eZToDPS&`mIhn=qH3V4^)5-CdDsi}3z){v6AJMD{SS>% zB_$Zr^rk?N4D{EUgJd9G_8V!pr)4 zgTo(zlg3Z#UIae|TD>s6YbjzpfV5l0EZ{t}u*MicQUs9bp*f7SP7_@0clDG4pe?ART z?7XjpBEK3RcL1pE>&qN$O66R(>_xwqr{~F&C-)yXV32+?JT?ooE~iXAYI78DqAHN- zG&ikMesXZ3Lk5jvWSWfmszu3p-Iy&cMW|^h)oxgY0kV1P(3CYf7cSgo;LRC@0#JF< zjdE9Tan%wqB;ZwbIegBnRTWmaaGI0XT*u2!0wtH{rB8THaLs7Q#a4^=7iigv1?TMT z>*K~_KTBG(8E(Zd*)68At3I;LB1AX@B!Ts& zguT7##HziPVZuUY*W{uOS~r?M!Y&I;ZWrOMyX30k!=RS66J6CRDzLZbbN?QSFuXON zE6XLm(lmT13@u#A@N%$t-P%adE7k=?e?32cwdvW%xw*M{dG5CIkC9#rqZ#B^eG6I( zT&}9(3%2lS%?2UzW}Lv=f`7r=zB*!6E4L#;LDsKZN6d>mh-UgLKy|h#0}&#W{sAmC zd_XAX3BM)3SNL=H(Kq|%RFXni!DsVc7})Z6?$Hho!PO9e@*z$7i?|iz44^K?pAq^~ zmE%9^Sc3zbTmHyj*?>rii`^RkNeGGYXSd&;6Wn>3@ zy|?DN4c;q+OXhVjbChyf$p>^n{O_S$2yFBw#8=l7d;;xyVA)Lh@8 z+#@NCQs`IUOiX4kt_kXg4&hAlnLS+*8yg$Fumd=ze=1?Gc(0D4<;!&yO4WF zq{NOD4L^D*8m0510wsGs$FAd~u`&2~Ap*)y6@qRG?0yAZ3#=tG-Xb{ZOajx16<7E#oC}u1_lOJ>}cbVl}C6g3T0gi1M&Zw zBlL?Ogg`r_QjsJ1qnq>t8uPT09Je(WLVd~te04GZM^|>c)zaPF z-OOyKAj;}5u+cyk&A#yX>C>l<@}3%qw~c6g!wB7c@Bo?jBLX9w*g4Jat?wu1P7E_q zWUVtA4qf1Ag1#ML>Nx<@Xl@60WR+!{`}?aynM7bN(-vWXK^5B_o*{%H#CIxTz~y8s zSoh0^TeHbvN4*Rtv6!>Xjg7rzGBI|c#{6(LAJ*FuBoFwdz|W}E)^K{dnD4%P>6&Mw zIGP@7fbCAN;p1bywP(D&wpOyU{|}^G>GdoshBXWHBT7e1RbEzhJf7Rf#L^+?TipsW62Z=QF4e z^!mfdh&p)lp&PZCNO%DF2tZl)9BHsb1wrR0Fb4=AQAEN0f(!wW^dNOcts<#C+geXg zuRiY?vvZ8W{$H#L@2 zFe(TDFz-HkH2K!ZzT=_!v!%c^L5vz?&W{nN-x@_o8_o2|;FR$W5E}qx76^E1`OdeI zO2DEZ&OFJCRr8nxbKPDGO@O#VM&A2tw{PE0%`Bm1LvajX?*kRFjTw0Zo%aaNDj4!6 zGX)>+25O>_bcmw@fAM%3CnOrs=FlaQH9W+p4d#O3vF?O=i5L?HiQiA$YBzmna2<;oRd=6C>@_XK7xYlZ~$06ZRHO?Sxh zWLLQZ*+py!`Fge0n(p+83M3~2eXALZV+Tba$ zac;}IbLUQ%4nYZ+O0{Gr_#uwX?U6NAVd5}IW%7GhZP!M`vJ-5J}qK@{6Ka z(t7}ztp0mr1f(6{PKnuQU}Gj7*(nBD#(E{pymM#TeKRo9Qs12iyrGjHWLtf9k(&41 zIQ66_^E)eT)kHK)KoA`J)K<scb)lud}w zZvXGfFtx=ZL=fX!l#`Q_m!~uzi4ZpM6rcec!{A!XXM86TMaBcWKKmV!{mwoM;DKP) zMJk0n)=`)nA1^W=31S;$mMr0L!>sdMq7(pcTZiyjxLc|DVQkFa#s*Z^6KHb?og9eR z7PKe9ra7aSKKWA7uHC)ipdykctGclnjL%kbCyvEvHgW4fh7AM(QZ|(;!FbLIQc0Qkvp{8}BQ zcn^q$B}b7UDhCoU^>v8vW00Mp??6*^>$uR&m?hx=3&Q)DcR?|z#=u*~U5mxPqzZVx z)(LgxvHJRyuFKpdWc4u*=ouBkU`xs>DBSA@v9Hu4{@;k&A?k^6IPTs(HJl4=#A~*Y zlzsz6{u}D5A^9KFLKwv@kh!RfRn4_+U8R`^7I+4nTo@hEb_FSEX~S3eOk46Cd&D&+ zt^!6S(o%^>_W3z)xlv_lMay0H$^N4({AB-$%_qj{YcR*~veT}&*_@njoAEl$N=rE{ zEv;5g&Doh5UOql+MXlt90R&!vaA{HR3F8N{fVpN*7a)$Ib0Z-cqAjKOoLo{{5aP(< zUso7gouaA;4vy|Yi2Mc0l(dL1_`3mIXhczGb1}A1gu+8C67CkAPH##z2sxj@i8Ri~ zj~~O3$~gH%0`D-{yAP41|7%(WVq_cPM>~`Zfn+*53MwJ6eFlSZ{`~n?D)4=2lL7_) zAjA~~X-7LmsJQ>O4dP5Sf}yR1p*ceX`X*!Wh|(Hu3f^7$sjO!K+?j@1Wo(Au;4FgY z5na7}#fo4rODq;UBS#m(vlA?OA5}y}X|68`fsH*A?h@hE!Jsf7JJxiz|6`oyKN5?z zYuk$-Y+jgn7|*Huh_Dwk1{4hN48*jAg7(hau2pk+@se@^8YiTleOpQNJ}a;qWv3Wyc$P*8~JzX})M zI{3?E%TMxnb_D}T+il)@_sbKMVX&NyLzhmi{$33nYVii@fEJOFk>L#lk&{GQc^!+s@8jTQs|U=y_sEzDG_z=Of?Zs>c6n-K6keF(ggG! z3?ZHAy~p^)3NzZx@bIRlrt_o{b^aijE89$q^0KnB^7He5)VgxP!UL#oa{Ra}Ob|fk zS9o5iAO3@*E*!}*43JKJQo!oDBk$^1HLSl7Nx%(k)E_;nTGrTh7P4W2f+l$m=4g?y zs^wfSD6^OQ9p(Xc(^Waf5RSA(1SZg$|P z=onjROkYLCNPg3@N_(qsqITkh>VJ^Uv#%oOz<_-YevHrjM_8)*nB}|pm6P_7#DB_e zB`rb@h|pL5%mq+KTfBit9XN2{?Ad__$5J)XfVCeZX9#+{IZRi$vD4E)PzPK*AsPTJ z4#pkie^-dKJJ+itUMRWn^U(g{w%v*)eVv8Lwd(+a_u3{S{{`=g&ZL$Y_p8l5f0sUX zi!(O5IUf<}FFBz>whkS5-PtJ=Qjy)WOB4j(0x)id>1>Mxf^}Jj!{go006HS_Jp=*S zr8rJzEfho)>V_KdtWY%v_1s`gTV6uAdC2vOo|K~!s=A3sz!6T40=ZNy7YYp(>A~AsqJ;$xrc*asRlbk_%%`iFI4HoBs_ za2h>Kuna-L{@Jt)Iq%PIZzA|ZM~*n1(-AI=-S6Q+P2pxIt?Jm*&$^Ts0(kjFOTbh8 zr%L#g#Ohj(xLd~|i2^Yp2%V=IWC|5l97D$^CnrNZ>hWWjrU%DF@L}S58Z#{x0D6Sy zTPhwYhI5H{L(Kj_as$GxXwyS7d0h+27b%G&C@%Huyh_Nyq=o?q*QXifZi~{f>ns8? z4bf}?HFpPa)!}E%L%K~94bB^qO9t50Fz|S4xx)(NQ6ko?OS=?YjCa*hg=}= z>=KCOgfIFCQ%7w{5;@;=!zl{fJ}m7`glOnqbaw_4ABNvxk3e-aFQK@&BLM6&@-@K8 zxkCdpJeq&$$SNr*jfR5l7{32>MaofNtwHantt~gnMHms%%vfR!}H1#j<_6({W#5VYY!=-M63+%3O#HGFc?# zPGC7?ik$zZ-RV4k&R@6sn49AoN~XvI!}an zyKznsjQml<{FIl-WmBaIgesM;k8F~BbpYI%-O3uwLMP^mjJd>iAy$tvZs}OjP`NSj zx!y)qXZkX`#xiQZyfW?zJgWhy&}X4N4FCHYD{ei@8CQc>4sw;z%bKTYBO#|Js$6K`7`qJGyHMc(%- z(6cYLD*Y74fBQk!PN{~dx(KH2=U@msu8I16`S%IklCOJ0>-v9%p4l#$Z2H%E4xd4o zIQ1LhxIYtS%FeeTy#=yaz>`?76SVbMyHK4a|6-)5K!O0tp%z`~d0mqOPpPmhFkhWn zs^@H)#0u-_$?au^`C{v^=K8@QtWKTt8KknHv3WZ?tXciDeAYeea{2lp^VcUZR&daO zsjUs#1Gy(*&kym#mabPeMu$N0(Y0|?+Osg|K<1(&D^uXiHX~wwZ(m<=-)ifm8tk@6 z|9Y`pniDo@j664ot*5H4Q;wuW%$;(V%Vn?vlyRJ)7E$dbrPBdcEbgJ}AL_pe8vwq3 zS(k}~$-l(dKY4zo*NwWBU}|oMm;%s8HTw9oD$7X-$c3q|^lNaZ4s?E+bQ7%R_OEs) zJZ^gcZlmcKcdv`iOiNqN$|-aFm}0>50JAbUd6GY*LflIfl{}FuS*`U8qy~Osq=xn5 zvXz==&YgpmE-^3gAa?&JOd8dZ^T(G2{8ya^LIv4ogREWddn5(EYG{;QsJd%k zJY5{C?Lo(VyIPm{%nCN*IN{V^SnRL_W>C=>#4P+lu zszF;@o3MuW9O9PQ+79!91YJ{81HO_kK#xA4xL-&soG?bShQP!7A* z`w8wCzsak@#dBBCtjxu*zIOr z@sbGgyg9pCQfOkrplk9T`7J%yb+O#8p|On_v6=M4z2tH0!)^{Tv{~~bKYJbIO{M7z z^>M*9hmy8%_&qTo)PDbtAKuFQi^ptVSSUu8|C?`k{Sd#iM5*RgoVe3hT34C4jm7MVp0+;^ljg<1N3onV~RdVX=mXugO3|an!U`e zq2$zCX_Hwz1ACn7_5s=J+BbC>BauTSVxn~Vb3j>u&5~_9XJ==j-jwi%3=JXLG}fMv zq}LYDA3u5&6fFWM{1C1^I9iSxq&$Fo(uEl!+0lK_@n)e_RaJ{K1q1~l-CcJcBsdTR zz2|d_!5cu-J8+V~JUVEx0^&7k=2;Qk{yhDB$m+*A#MsWx%fUf5ca?>!YfgN8H!S~%;ICf2x=S<*{2^metT#R-B?Y{is=W)+ zL~ZVH9M=Vu4EBqvs;Wg+3-k!sS#51un@>x-J--90l zX=~oI{tn8yZ1ftppyYZ1s|>$C*`+7ni)#99J}v%!2}jq~lzMr2 zd05&4Bnq^U`m#_94cr0NDL_pCAvT{Eel&^)o*Bf}Qeh`AI+}$&Qg3B+@+6o)Zw6f0 z?#8K8nxMvaomG^T1#R=@{zM3)Y}e4}Y@Q7$#Jah-z*6e-Q&|wd1y5W@PtP9K@7Ukp zOd=6<#GKmzFX!eel1vXDUa7f7N(!hRJkZ~4Jh`TJw(woxONjSQ$MK(Av1@+n&YW{A ztOidCB{7Hk_lB@S$QxN{{!zq0dK)TaEdN7cq+B)W#GQmO(e)nO^{g%LdWb?yBWNCu zQ)kO(N0zRyc>S7w*$m`1l~~iH4y0WTHIh2fB5S1Nl=H{;b-^l7SXFP1e3Z`V(|vLy zqtHYZ%6M4}h7~K5mY%MTV|pV+I}V1f2)?zo6)2J(hM|C${o{|&y_-OX z+Xl^N=N6FY!2GnfwwB^W7M4MbEi5R2ZL>E}t&rOT{oAD2^9k@5{P=jc7Gk594KyUn zN7b~q8$-$x-3>dqAo_BO7{k^X?C1hu$S;2dA{5f^KywD8(jGLcUg6{$M){#+QHXWU z2|*i?sx5wlY9#Oer>0HWVvm~kScyDfBoxPP*X%aBd7J1W!`s>(3-8wl{xCS${=|i9 z6D8EG{==QnCq%74>H8rTU>57dx56)Rw?p`BO^uG5n;VD;K0XY9ER>ON3hWlhIuKwt zF!eZ16@)pn*6P1>Q9Jda4JZ`n1oF86-*nWON!0~N`>}>q8_=X-2LSFa_^d#3^Jf5n zLctnJy=`f6hlRU=fdN?9Q`o-Uha8MbE81R!?>T<|YC#~m?DF0a%-Z*G)Xr=PnVPF= zWZc&iiMxHhAP}>Hwej5w-zN*Oe|fDG#_i9i|2Z`~3)W=2t3>$V$VdnD=?n&_8X|be z+PSOw0v(3-rCt}Tci(c#WEUg!x;AdvVYi7g&$;S$oQ>OKY!z$0$1OB*UgG@=1x*1Ku&T))F|}Ai)E0(L%#oE*(}epta&aNrD?(uF9{2%$^hn+SCB- zG%6}02vt)Fw*Fz&ND~dZw7~lH%8;M1OLum3oQRUQu8UF5ZVW~cg47ti*$MY8NWnz*+o!qnm^Q=;Q`}a=!cK@X`$|D$XG6cL-aYv@}(@_ zZguNM{hONJaTX2d=E|3IgCGMql>c1AKQ>89P^)bY=q<3N7OR{4YH!B>&qf4@D#o9lW$5+0ASR z*GS;5wUa>+5!634#s%@`cIY8E(r|kK2 zc=j!YZeySrnZQyjTrlu|ULgVT2xUY$P`%>siZ1L0^HD7@ji7T)>ttuVxGBfK*<}d-i0gcq*?0S&GR|u&gHjMF0af0`aS}P=>4hsTp?9I%~&Yam2 zz8|#kii&v}4N?rkn&5WbPf592ehM}TKF`aV!bH&$i(#RhwTIYbcG_N+2nySApj9y5SB z^)vxj$ydg>U+BB-`W{a@Us*AnOYEG}LGHt7bBX5(qwnikG?|2vt_qd_gRc*_+nc`B zrnsfHp1?3&D>h5e+>3GLuGjEYuCj1mAHA;ReP5F-?*4+faFBRBY3jju$^k#mQ(wRL+v-@q_EEja=Z;M zyam|Nb5Z2K(MgL&9$=(H;yYO{Sso&T+md!X_fd%8(4*eluq3w+1TpP$%b3uh-i&qAS?sR-xyze1@Ng){2wH}&PLegrJt zD^LBy09hU~xdf=N-1mT)4U9fg;b&|r{GaLYBNcubNqZWiBro64(P47*=us?|2zey$ zc}VVpD|6z+37DZ7ecAOI@Tm(bu=Ea{=~nAXV9wn%+y8=l0P7gTCw0JDCGvPZCVRo( zJ)*kc;p8;rnhllz9q54%Dqu_{V>-b1e|mW-MEGvKTj>{`iP(|V17l-)Vqy6sMhCU~P_r?BI;Z(+x=Cyb+F6-!(6G?M)SfD7M9XMs4 zcsnu@QuBHV2vey$Zgf<9`sd>%VqRl^Y$p$@4ZPOGmNtiqxg>E5mAmHPYDD|*9 z-eiF6ZAf<>_F2+MK9*9yf$+D)@^OW~DtdsgLu8~6Y*YxW0epmklL<>1fK})-xkOlM z8bZ={6atG5847nK3#qEAc>Y~uTiZOu0Hr&TrheRKTO@&?Xo-^C8dgfmmPT>Uppn^F?L(0?iTG<(ofz1Khv) v!)Jwt{q5iU;p?J>#`*un%cK>sIO{rIoSHgZ#0?+8h1$1UcUOj{*`@yl%@`U- literal 0 HcmV?d00001 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 61fac3d..eb2dda4 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -2,13 +2,14 @@ add_executable(first_example first_example.cpp) add_executable(second_example second_example.cpp) add_executable(boundary_example1D boundary_example1D.cpp) add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp) +add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp) add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) - target_link_libraries(first_example tug) target_link_libraries(second_example tug) target_link_libraries(boundary_example1D tug) target_link_libraries(FTCS_2D_proto_example tug) +target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(reference-FTCS_2D_closed tug) -# target_link_libraries(FTCS_2D_proto_example easy_profiler) \ No newline at end of file +# target_link_libraries(FTCS_2D_proto_example easy_profiler) diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp new file mode 100644 index 0000000..442b64a --- /dev/null +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -0,0 +1,77 @@ +/** + * @file FTCS_2D_proto_example.cpp + * @author Hannes Signer, Philipp Ungrund + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach + * and constant boundary condition + * + */ + +#include + +int main(int argc, char *argv[]) { + + // ************** + // **** GRID **** + // ************** + + // create a grid with a 20 x 20 field + int row = 64; + int col = 64; + int n2 = row/2-1; + Grid grid = Grid(row,col); + + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); + + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row,col,0); + concentrations(n2,n2) = 1; + concentrations(n2,n2+1) = 1; + concentrations(n2+1,n2) = 1; + concentrations(n2+1,n2+1) = 1; + grid.setConcentrations(concentrations); + + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); + + + // ****************** + // **** BOUNDARY **** + // ****************** + + // create a boundary with constant values + Boundary bc = Boundary(grid); + + // (optional) set boundary condition values for one side, e.g.: + bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + + + // ************************ + // **** SIMULATION ENV **** + // ************************ + + // set up a simulation environment + Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + + // (optional) set the timestep of the simulation + simulation.setTimestep(1000); // timestep + + // (optional) set the number of iterations + simulation.setIterations(5); + + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); + + // **** RUN SIMULATION **** + + // run the simulation + simulation.run(); + + return 0; +} diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp new file mode 100644 index 0000000..80d4112 --- /dev/null +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -0,0 +1,77 @@ +/** + * @file FTCS_2D_proto_example.cpp + * @author Hannes Signer, Philipp Ungrund + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach + * and constant boundary condition + * + */ + +#include + +int main(int argc, char *argv[]) { + + // ************** + // **** GRID **** + // ************** + + // create a grid with a 20 x 20 field + int row = 64; + int col = 64; + int n2 = row/2-1; + Grid grid = Grid(row,col); + + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); + + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row,col,0); + concentrations(n2,n2) = 1; + concentrations(n2,n2+1) = 1; + concentrations(n2+1,n2) = 1; + concentrations(n2+1,n2+1) = 1; + grid.setConcentrations(concentrations); + + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); + + + // ****************** + // **** BOUNDARY **** + // ****************** + + // create a boundary with constant values + Boundary bc = Boundary(grid); + + // (optional) set boundary condition values for one side, e.g.: + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + + + // ************************ + // **** SIMULATION ENV **** + // ************************ + + // set up a simulation environment + Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + + // (optional) set the timestep of the simulation + simulation.setTimestep(1000); // timestep + + // (optional) set the number of iterations + simulation.setIterations(5); + + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); + + // **** RUN SIMULATION **** + + // run the simulation + simulation.run(); + + return 0; +} diff --git a/proto/FTCS.ipynb b/proto/FTCS.ipynb index 4708428..b64c934 100644 --- a/proto/FTCS.ipynb +++ b/proto/FTCS.ipynb @@ -458,7 +458,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.9.7" }, "orig_nbformat": 4 }, diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index 947ddc8..6210cd8 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-05-11 17:31:41 delucia" +## Time-stamp: "Last modified 2023-07-31 14:26:49 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 diff --git a/scripts/HetDiff.R b/scripts/HetDiff.R new file mode 100644 index 0000000..417e606 --- /dev/null +++ b/scripts/HetDiff.R @@ -0,0 +1,237 @@ +## Time-stamp: "Last modified 2023-07-31 16:28:48 delucia" + +library(ReacTran) +library(deSolve) +options(width=114) + +## harmonic mean +harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & y have different lengths") + 2/(1/x+1/y) +} + +## harm(0, 1) ## 0 +## harm(1, 2) ## 0 + + +############# Providing coeffs on the interfaces +N <- 11 # number of grid cells +ini <- 1 # initial value at x=0 +N2 <- ceiling(N/2) +L <- 10 # domain side + +## Define diff.coeff per cell, in 4 quadrants +alphas <- matrix(0, N, N) +alphas[1:N2, 1:N2] <- 1 +alphas[1:N2, seq(N2+1,N)] <- 0.1 +alphas[seq(N2+1,N), 1:N2] <- 0.01 +alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001 + +image(log10(alphas), col=heat.colors(4)) + +r180 <- function(x) { + xx <- rev(x) + dim(xx) <- dim(x) + xx +} +mirror <- function (x) { + xx <- as.data.frame(x) + xx <- rev(xx) + xx <- as.matrix(xx) + xx +} + +array_to_LaTeX <- function(arr) { + rows <- apply(arr, MARGIN=1, paste, collapse = " & ") + matrix_string <- paste(rows, collapse = " \\\\ ") + return(paste("\\begin{bmatrix}", matrix_string, "\\end{bmatrix}")) +} + + +cat(array_to_LaTeX(mirror(r180(alphas)))) + + + +r180(alphas) + +filled.contour(log10(alphas), col=terrain.colors(4), nlevels=4) + +cmpharm <- function(x) { + y <- c(0, x, 0) + ret <- numeric(length(x)+1) + for (i in seq(2, length(y))) { + ret[i-1] <- harm(y[i], y[i-1]) + } + ret +} + +## Construction of the 2D grid +x.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +y.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +grid2D <- setup.grid.2D(x.grid, y.grid) + +D.grid <- list() + +# Diffusion on x-interfaces +D.grid$x.int <- apply(alphas, 1, cmpharm) + +# Diffusion on y-interfaces +## matrix(nrow = N, ncol = N+1, data = rep(c(rep(1E-1, 50),5.E-1,rep(1., 50)), 100) ) +D.grid$y.int <- t(apply(alphas, 2, cmpharm)) + +dx <- L/N +dy <- L/N + +# The model equations +Diff2Dc <- function(t, y, parms) { + CONC <- matrix(nrow = N, ncol = N, data = y) + dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + return(list(dCONC)) +} + +## initial condition: 0 everywhere, except in central point +y <- matrix(nrow = N, ncol = N, data = 0) +y[N2, N2] <- ini # initial concentration in the central point... + +## solve for 10 time units +times <- 0:10 +outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL, + dim = c(N, N), lrw = 1860000) + +outtimes <- c(0, 4, 7, 10) + +## NB: assuming current working dir is "tug" +cairo_pdf("doc/images/deSolve_AlphaHet1.pdf", family="serif", width=12, height=12) +image(outc, ask = FALSE, mfrow = c(2, 2), main = paste("time", outtimes), + legend = TRUE, add.contour = FALSE, subset = time %in% outtimes, + xlab="",ylab="", axes=FALSE, asp=1) +dev.off() + +## outc is a matrix with 11 rows and 122 columns (first column is +## simulation time); +str(outc) + +## extract only the results and transpose the matrix for storage +ret <- data.matrix(t(outc[ , -1])) +rownames(ret) <- NULL + +## NB: assuming current working dir is "tug" +data.table::fwrite(ret, file="scripts/gold/HetDiff1.csv", col.names=FALSE) + + + + +#################### 2D visualization + +## Version of Rmufits::PlotCartCellData with the ability to fix the +## "breaks" for color coding of 2D simulations +Plot2DCellData <- function(data, grid, nx, ny, contour = TRUE, + nlevels = 12, breaks, palette = "heat.colors", + rev.palette = TRUE, scale = TRUE, plot.axes=TRUE, ...) { + if (!missing(grid)) { + xc <- unique(sort(grid$cell$XCOORD)) + yc <- unique(sort(grid$cell$YCOORD)) + nx <- length(xc) + ny <- length(yc) + if (!length(data) == nx * ny) + stop("Wrong nx, ny or grid") + } else { + xc <- seq(1, nx) + yc <- seq(1, ny) + } + z <- matrix(round(data, 6), ncol = nx, nrow = ny, byrow = TRUE) + pp <- t(z[rev(seq(1, nrow(z))), ]) + + if (missing(breaks)) { + breaks <- pretty(data, n = nlevels) + } + + breakslen <- length(breaks) + colors <- do.call(palette, list(n = breakslen - 1)) + if (rev.palette) + colors <- rev(colors) + if (scale) { + par(mfrow = c(1, 2)) + nf <- layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(4, + 1)) + } + par(las = 1, mar = c(5, 5, 3, 1)) + image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1, asp = 1, + breaks = breaks, col = colors, axes = FALSE, ann=plot.axes, + ...) + + if (plot.axes) { + axis(1) + axis(2) + } + if (contour) + contour(unique(sort(xc)), unique(sort(yc)), pp, breaks = breaks, + add = TRUE) + if (scale) { + par(las = 1, mar = c(5, 1, 5, 5)) + PlotImageScale(data, breaks = breaks, add.axis = FALSE, + axis.pos = 4, col = colors) + axis(4, at = breaks) + } + invisible(pp) +} + +PlotImageScale <- function(z, zlim, col = grDevices::heat.colors(12), breaks, + axis.pos = 1, add.axis = TRUE, ...) { + if (!missing(breaks)) { + if (length(breaks) != (length(col) + 1)) { + stop("must have one more break than colour") + } + } + if (missing(breaks) & !missing(zlim)) { + breaks <- seq(zlim[1], zlim[2], length.out = (length(col) + 1)) + } + if (missing(breaks) & missing(zlim)) { + zlim <- range(z, na.rm = TRUE) + zlim[2] <- zlim[2] + c(zlim[2] - zlim[1]) * (0.001) + zlim[1] <- zlim[1] - c(zlim[2] - zlim[1]) * (0.001) + breaks <- seq(zlim[1], zlim[2], length.out = (length(col) + 1)) + } + + poly <- vector(mode = "list", length(col)) + for (i in seq(poly)) { + poly[[i]] <- c(breaks[i], breaks[i + 1], breaks[i + 1], + breaks[i]) + } + if (axis.pos %in% c(1, 3)) { + ylim <- c(0, 1) + xlim <- range(breaks) + } + if (axis.pos %in% c(2, 4)) { + ylim <- range(breaks) + xlim <- c(0, 1) + } + plot(1, 1, t = "n", ylim = ylim, xlim = xlim, axes = FALSE, + xlab = "", ylab = "", xaxs = "i", yaxs = "i", ...) + for (i in seq(poly)) { + if (axis.pos %in% c(1, 3)) { + polygon(poly[[i]], c(0, 0, 1, 1), col = col[i], border = NA) + } + if (axis.pos %in% c(2, 4)) { + polygon(c(0, 0, 1, 1), poly[[i]], col = col[i], border = NA) + } + } + box() + if (add.axis) { + axis(axis.pos) + } +} + +cairo_pdf("AlphaHet1.pdf", family="serif", width=8) +par(mar = c(1,1,1,1)) +Plot2DCellData(log10(mirror(alphas)), nx=N, ny=N, nlevels=5, palette = terrain.colors, contour=FALSE, plot.axes=FALSE, + scale = F, + main=expression(log[10](alpha))) +text(3,8,"1") +text(8,8,"0.1") +text(3,3,"0.01") +text(8,3,"0.001") +# title("Diff. Coefficients (log10)") +dev.off() + diff --git a/scripts/gold/HetDiff1.csv b/scripts/gold/HetDiff1.csv new file mode 100644 index 0000000..2c4a25f --- /dev/null +++ b/scripts/gold/HetDiff1.csv @@ -0,0 +1,121 @@ +0,1.15723009391899e-05,0.000573422585977779,0.00271483227696389,0.00592017488450261,0.00919024372055321,0.0119992110381186,0.01421611640169,0.0158871549414673,0.0171115136326332,0.0179901833258047 +0,4.32496365970161e-05,0.00110003761455522,0.0038588402926536,0.00723956931728921,0.0103665385368182,0.0129275133140941,0.0149018727052112,0.0163722917415511,0.0174425542133075,0.0182067473313329 +0,0.000156410551955808,0.0022501860978233,0.00596846775043427,0.00951979828485079,0.0123439077799255,0.0144650737976258,0.0160242190667615,0.0171550743320483,0.0179654704632021,0.0185369538227271 +0,0.000449900131743341,0.00397802650798475,0.0085351147388942,0.0120552559726064,0.0144482796517419,0.0160547847343675,0.0171511707942764,0.0179099936829702,0.0184375120234885,0.0188001645647094 +0,0.00096247784605253,0.00575161551352107,0.0106416083625011,0.0139078265679269,0.0158706072121663,0.0170493964826529,0.0177834857154376,0.0182584391833927,0.018573228563985,0.0187811039627322 +0,0.00139148060908994,0.00653728609085073,0.0112226963812925,0.0141895931925161,0.0159026822227563,0.0168925530607027,0.0174883152516615,0.017865325174865,0.0181135647391488,0.0182789062492789 +0,7.88755181655377e-05,0.000890855516524495,0.0025489601065516,0.00451888907818989,0.00640234215109739,0.00803400181298965,0.00938380805959316,0.0104782803692741,0.0113602330285538,0.0120715601783272 +0,1.96835046894113e-06,5.20032828635902e-05,0.000248131836356188,0.000630746867714403,0.00117274437938983,0.00182064447664627,0.00252211668991391,0.00323643607872941,0.00393564702867779,0.00460249317726503 +0,4.05263620284678e-08,2.41553600349545e-06,1.87665958722685e-05,6.74081501157623e-05,0.000163264089823103,0.000313023273063137,0.000515631508985307,0.000764816320228696,0.00105167332970099,0.00136658044925295 +0,7.12888460576857e-10,9.36408593803959e-08,1.16671799220012e-06,5.86484773654181e-06,1.84017858631692e-05,4.34716444970184e-05,8.51739231661631e-05,0.000146375097555955,0.000228524220562364,0.000331765311847696 +0,1.11266213313442e-11,3.21569209939452e-09,6.45997582539074e-08,4.58674913815976e-07,1.88637915450085e-06,5.56437226690733e-06,1.31531030829188e-05,2.6577288329037e-05,4.78110315422556e-05,7.86848274922097e-05 +0,4.35985871981009e-05,0.0011173013784839,0.00393465640259938,0.00738648425859891,0.0105643563174128,0.0131485926879987,0.0151246465122267,0.016583690239707,0.0176362908990393,0.0183809352240671 +0,0.000162940275870536,0.00214415472789055,0.0055971341202626,0.00904266600518215,0.0119317484324307,0.01418453738493,0.0158744170951516,0.0171103507405317,0.0179967639568099,0.0186199630282087 +0,0.000589272728343116,0.00438867743664915,0.00866893626959719,0.0119151946783314,0.0142428604635392,0.0159134847420158,0.0171142649822806,0.0179721981758363,0.0185775454158148,0.0189955063245207 +0,0.00169526769334991,0.00776694628989879,0.0124247154337811,0.0151382671118746,0.0167367107498905,0.0177366662921586,0.0183940071716407,0.0188368236363926,0.0191342578330509,0.0193277227804878 +0,0.00362839452566188,0.0112526737744286,0.0155521246084951,0.0175590393037948,0.0184985843185778,0.0189571544442947,0.0191923705608016,0.0193168273546316,0.0193791769704849,0.0194019806796609 +0,0.00525177970919126,0.0128429513741672,0.0165189505188843,0.0180778870279157,0.0187197075826292,0.018968965691842,0.0190518322037241,0.0190652772723224,0.0190482483782045,0.019016041284815 +0,0.000385153961868949,0.00229577940499986,0.00480460165448024,0.00711204718132809,0.0089909576548431,0.0104472162269959,0.0115548772926463,0.0123946614477843,0.0130350658282232,0.013528343285464 +0,1.16858258297074e-05,0.000164584714174317,0.000570116218846185,0.00119271144355875,0.00194723110061389,0.00275621955267251,0.00356441433332025,0.00433765719695555,0.00505760032168857,0.00571654514201105 +0,2.81208380577636e-07,8.97423573740115e-06,5.03422264992399e-05,0.000147427646174569,0.000310274427544454,0.000536832367026033,0.000817602280027168,0.00113989443932459,0.00149069663261676,0.00185828031564749 +0,5.63399964199174e-09,3.96964588417007e-07,3.55564115271452e-06,1.44716592579598e-05,3.91537025779409e-05,8.28348829677685e-05,0.000148977185800308,0.000239043494747069,0.00035275053097983,0.000488513445484224 +0,9.84398346768558e-11,1.5293568296117e-08,2.20334810416761e-07,1.26108145384204e-06,4.44865566281344e-06,1.16871019830672e-05,2.52195633384727e-05,4.73253788486059e-05,8.00591501276681e-05,0.000125072519296284 +0,0.00015939225040812,0.00233357034425003,0.00625562530069953,0.0100135877171516,0.0129676882444301,0.0151374162693059,0.0166874294020243,0.0177764984876221,0.0185308441084237,0.0190433426894296 +0,0.000595680647430986,0.0044801608514742,0.00890743169921897,0.0122772025886856,0.0146734676860174,0.0163633416618241,0.0175504470842378,0.0183769265178609,0.018943746718196,0.019322547574371 +0,0.00215448195385171,0.0091767592756406,0.0138199838635887,0.0162232592499362,0.0175799863177681,0.0184333811027871,0.0190005234819769,0.0193808295698239,0.0196287328140203,0.01977995018463 +0,0.00619992143295609,0.0162622253524985,0.0198667054266907,0.0207089785376823,0.0207820327860086,0.0206824193920932,0.0205605370395911,0.0204471756023412,0.0203412141522796,0.0202389042951594 +0,0.0132779136133337,0.0236231555894036,0.0250044710521133,0.0242144479160946,0.0231938485225262,0.0223380815152074,0.0216793603038312,0.0211797341820616,0.020794392809717,0.0204895366674834 +0,0.0192489696423083,0.0271184370488408,0.026840217635384,0.0252838506050641,0.023849962885927,0.0227245048675836,0.0218698389880889,0.0212219274277141,0.0207240679736698,0.020334308600127 +0,0.00188379770767293,0.00657190925731914,0.0104892667525962,0.0130700680836906,0.0146349533533907,0.015535980568487,0.0160269602381333,0.0162728756742906,0.0163767130755721,0.016400710298742 +0,7.01829175866307e-05,0.000581093706531354,0.00153370940121776,0.00269094081762158,0.00386803514667116,0.00496316154582596,0.00593261862345271,0.00676642929553236,0.0074719703393928,0.00806413175855109 +0,1.97847089603411e-06,3.69925955384929e-05,0.000157406497985584,0.000384860592255597,0.000709846476804252,0.0011079382186538,0.00155170108970297,0.00201681216519113,0.00248426641823041,0.00294058581228933 +0,4.51078087575571e-08,1.85104501980116e-06,1.25057013282593e-05,4.22892094350572e-05,9.98395815944757e-05,0.000189787308802236,0.000312690644970973,0.000466034643916367,0.000645463690913993,0.000845828528229199 +0,8.80252798763183e-10,7.92215439962928e-08,8.56765592551646e-07,4.05803525996784e-06,1.2449180396882e-05,2.9295973763362e-05,5.7746988001425e-05,0.0001003733295813,0.000158928915289496,0.00023429575744269 +0,0.000466368704740309,0.00424791180746615,0.00926944539909274,0.0131713216643424,0.0157625616947689,0.0174133722552654,0.0184577152890545,0.0191156877417616,0.0195250042740468,0.0197700479556243 +0,0.00174300015594418,0.00815930380362134,0.0132116161354557,0.016172754442502,0.0178696343095053,0.0188636958720877,0.0194550603754561,0.0198042367684176,0.0200009359043926,0.0200979478510297 +0,0.00630533087106852,0.0167264504542955,0.0205352907522532,0.021433296617797,0.0214908287288404,0.021342901494693,0.0211590163834043,0.0209808421912509,0.020813382900871,0.0206559438303523 +0,0.0181516683021334,0.0296877506877057,0.0296188711707772,0.0274991387847508,0.0255697170579367,0.0241218925670229,0.0230708035427124,0.0223012794968431,0.0217225548818433,0.0212745536462067 +0,0.0389067518217964,0.0432743493327086,0.037525495598965,0.0324527348417642,0.0288543692365208,0.0263663031569143,0.0246230583185604,0.0233728213888174,0.0224524193870472,0.0217568092848378 +0,0.0565396239160925,0.0500897261990562,0.0408387347867424,0.0344846493850116,0.0302543240032215,0.0273643802632421,0.0253287379325492,0.0238539136707014,0.0227584403765214,0.02192583980509 +0,0.00794755170946868,0.0172628249510397,0.0219912124473273,0.0237906919799591,0.0240591222673087,0.0235982544067714,0.0228369786001951,0.021996032718532,0.0211823304839058,0.0204423856134907 +0,0.000373684524467729,0.00190283496658901,0.00395283626547245,0.00595433814031928,0.00766198835968331,0.0090150942497423,0.0100370905098692,0.0107818312810106,0.0113083452714233,0.0116697843722264 +0,1.24794246503553e-05,0.000141270114233274,0.000466648605475989,0.000969170133300241,0.00158727378005724,0.00225743530947124,0.00292973311195242,0.00357075743939607,0.00416127927601125,0.00469265749237697 +0,3.25056776407433e-07,7.95461921640221e-06,4.12089493515841e-05,0.000117250753288648,0.00024402042254145,0.000420314286347587,0.000638908995931823,0.000889809192557268,0.00116261453581859,0.00144792147140295 +0,7.0854379238238e-09,3.75524392933622e-07,3.08321489667331e-06,1.21925554449326e-05,3.27749649029249e-05,6.95563221553643e-05,0.000126018794389096,0.000204052811537285,0.000304018569398211,0.000425019521621274 +0,0.00102836224521133,0.0064274670666241,0.0121491816508084,0.015968397765568,0.0181492003638804,0.019314880428558,0.0199102032969297,0.0201928751968812,0.0203044557297043,0.020320071184116 +0,0.00384385388740204,0.0123515070854559,0.0173313777131353,0.0196324133235912,0.0206081184457373,0.0209612687545928,0.0210258040457677,0.0209599611581841,0.0208374482363417,0.0206927731773186 +0,0.0139084634226096,0.0253420008838306,0.0269858208570455,0.0260863431776468,0.0248662406090429,0.023805495300228,0.0229587141516884,0.0222941985053244,0.0217679373130175,0.021344966099159 +0,0.0400577489259572,0.0450598366852064,0.0390556903095381,0.0336304680408348,0.0297597604093353,0.0270807237205804,0.0252035706273644,0.0238573134643687,0.0228662562144913,0.0221177566182138 +0,0.0859596622787531,0.0659645822190159,0.0498440895724295,0.0400638374763904,0.0339433683852068,0.0299353070092985,0.0272039578032914,0.0252772077773997,0.0238777953326208,0.0228341669415018 +0,0.125398084679699,0.0772427043755313,0.0551767906604444,0.0434193651923798,0.0363278914606691,0.0317025734216055,0.0285266548704347,0.0262622056865752,0.0246012510386757,0.0233529211513904 +0,0.0279719026476077,0.0399912836159061,0.041824217737992,0.0399960375533372,0.0370306079722165,0.0339334669014817,0.0310904350363235,0.0286238301128732,0.0265419722346029,0.0248094228214593 +0,0.0016878543861307,0.00546802879678612,0.00908467643641864,0.0118685267101036,0.0137867836918316,0.014992473171871,0.0156668769769553,0.0159679687584234,0.0160191464240754,0.0159109388085067 +0,6.66210210366172e-05,0.000466311941721689,0.0012074539008488,0.00214399617569523,0.00313678410870533,0.0040909928004654,0.0049535861902635,0.00570178710489221,0.00633199242057907,0.00685150908457396 +0,1.96581789972705e-06,2.90655193631781e-05,0.000116181708737399,0.000279461298722284,0.000515297156287373,0.00080891264649635,0.00114182790936584,0.00149639884138984,0.0018579271762422,0.00221530254986257 +0,4.73511582533572e-08,1.48929389156945e-06,9.32003158799015e-06,3.08850626482769e-05,7.30754947692817e-05,0.000140629698309681,0.000235669226477032,0.00035796975388455,0.000505546298234415,0.000675276938000184 +0,0.0015836697501133,0.00789553113761796,0.0138353209275783,0.0174900423844659,0.0193829171959464,0.0202511267409907,0.0205839676229735,0.0206522039767607,0.0205960670144616,0.0204840229205406 +0,0.00592024795026366,0.015178056357401,0.0197491887189647,0.0215219374447417,0.0220321907378727,0.0220033162511111,0.0217646356223314,0.0214641136669779,0.0211632832330552,0.0208849932589192 +0,0.0214261698026196,0.0311627748380295,0.030789710824671,0.0286485315817415,0.0266434190236869,0.0250513323715184,0.023828703885917,0.0228921858496036,0.0221672538845208,0.0215984913078681 +0,0.0617381544584604,0.0554953966664163,0.0446769294750762,0.0370604778930011,0.0320153173963507,0.0286236282907804,0.0262786287048439,0.0246098514264461,0.0233900725628177,0.0224762767353041 +0,0.132662095447386,0.0815715457115294,0.0573620962224714,0.0444656436513204,0.0367970084501995,0.0318900160341415,0.0285856373696197,0.0262711290074283,0.0245996244789251,0.0233597987792184 +1,0.194600594053136,0.0967188378926377,0.0645033312825157,0.0489815255136541,0.0400071807557426,0.0342746856436233,0.0303863403322649,0.0276362424095596,0.0256323186910777,0.0241347191227446 +0,0.072774737081717,0.0731608818750184,0.0649513569830106,0.0562393519385766,0.048657668424965,0.0424321620057059,0.0374285709142598,0.0334382577022264,0.0302593069735787,0.0277213581035032 +0,0.00546229967617301,0.0117769944796152,0.0160960402630076,0.0186595966553889,0.0199729770538671,0.020458608892436,0.0204223397121343,0.0200740130964532,0.0195534451842866,0.0189514347157146 +0,0.000246612705267533,0.00110781620007196,0.00231152127659001,0.00359442519681962,0.00480227348450595,0.00586145140404475,0.00674755083035105,0.00746338966694999,0.00802538295172145,0.00845539933096027 +0,8.01024172201437e-06,7.40436516662321e-05,0.000234915899100426,0.000490064251651924,0.000819680663206465,0.00119845837778001,0.00160264998868245,0.00201295787243647,0.00241514785548912,0.00279964947569617 +0,2.07955309779492e-07,4.0104486961948e-06,1.96923787162901e-05,5.61707039451707e-05,0.00011992118154755,0.000214121267381263,0.000338952769617425,0.000492335832038914,0.000670725912883758,0.000869802668989746 +0,9.15481631874651e-06,0.000117693874573294,0.000377019455101549,0.000741577814721888,0.00115788487944848,0.00159137882061179,0.00202342112731787,0.00244488663055897,0.00285171452271575,0.00324239719313011 +0,4.46105775530757e-05,0.000303984477552903,0.000717806961735601,0.00118873254936734,0.00166857393355206,0.00213667568842936,0.00258503807646796,0.00301133151201256,0.00341571030049681,0.00379926644315538 +0,0.000217726325331875,0.000873770408536327,0.00159096186837063,0.00224649712918517,0.00282908355295164,0.00334893383349626,0.00381769528011674,0.00424489367607894,0.00463786438443504,0.00500213301923337 +0,0.000918933831995098,0.00232435924877404,0.00343741337091249,0.00429315823548186,0.00497145196445075,0.00552787579177755,0.00599805857236328,0.00640506369246164,0.00676421367084299,0.00708600072229698 +0,0.00331122058413453,0.00568447123669356,0.00712540944284069,0.00809397304876091,0.00879084024979697,0.0093162979589995,0.00972657995173228,0.0100559006360597,0.010326180081839,0.0105521545670149 +0,0.0100043670600405,0.0127097661456827,0.0139937515689219,0.0147091809677894,0.0151257626062355,0.0153636070507909,0.0154869837094123,0.0155343002312855,0.0155299222027964,0.0154901124905561 +0,0.000136637116798085,0.000338905367123057,0.000531352383440882,0.000703792701127588,0.000856724436487964,0.000992875271716,0.00111515974854185,0.00122613935565325,0.00132792886232481,0.00142222359205319 +0,5.41197942669333e-06,2.66179151421015e-05,6.06837139232639e-05,0.00010289135389313,0.00014958545120669,0.00019826296968678,0.000247306756529671,0.000295722537183644,0.000342937360699378,0.000388655525576309 +0,1.72186371215188e-07,1.70696816511905e-06,5.79640506738781e-06,1.29195485396102e-05,2.30515604056451e-05,3.58992272238309e-05,5.10570165270148e-05,6.80982556623907e-05,8.66227559703461e-05,0.000106279145950884 +0,4.32364439009495e-09,8.64589599304522e-08,4.3955011662034e-07,1.29702112369114e-06,2.8633412989769e-06,5.28573284839916e-06,8.65071350621842e-06,1.29921291564899e-05,1.83031949023414e-05,2.45484723579094e-05 +0,9.12457420765116e-11,3.74315585927427e-09,2.9038454685416e-08,1.15682654773554e-07,3.22205377431422e-07,7.18623792678699e-07,1.37865133447057e-06,2.37328126021504e-06,3.76631554161963e-06,5.61170605601181e-06 +0,2.08011169409904e-08,6.36518490608328e-07,3.48383663648875e-06,1.01031069547314e-05,2.13412103139757e-05,3.74893709216357e-05,5.85171433417359e-05,8.42328798186352e-05,0.000114372709309901,0.000148647128283624 +0,1.23164013985779e-07,2.01619957543725e-06,8.06746622712326e-06,1.9408951787608e-05,3.62807655737012e-05,5.85218420212134e-05,8.58109894680849e-05,0.000117776723517294,0.000154044674846032,0.00019425785338382 +0,7.38055239000408e-07,7.13661815273176e-06,2.1957201331878e-05,4.48051934222451e-05,7.46425998373469e-05,0.000110486321507166,0.000151514246571192,0.000197051618892773,0.000246540311467432,0.00029951243350813 +0,3.94931915614762e-06,2.37664477362685e-05,5.85078377488893e-05,0.000104356634872865,0.000158457684921854,0.000218839976363591,0.000284118381194624,0.000353286447661602,0.00042558888273829,0.000500441217237315 +0,1.87536920881212e-05,7.42297836614388e-05,0.000150953066786234,0.000240330072584709,0.000337698410142427,0.000440229993907094,0.00054607757879758,0.000653971107402446,0.000763006862111349,0.000872522857886958 +0,7.67060564087528e-05,0.000211583511011325,0.000364290214247424,0.000523711941940329,0.000684842914248833,0.000844968149576811,0.00100249413331233,0.00115645811127351,0.00130628569191617,0.00145165008061455 +0,1.29198147797745e-07,7.26337646539414e-07,1.8789774469963e-06,3.58894576036159e-06,5.83743184885366e-06,8.59992729697426e-06,1.1850973027563e-05,1.55658536023551e-05,1.97211620171109e-05,2.42949568717564e-05 +0,1.89355863862761e-09,2.0149347443802e-08,7.32356468455249e-08,1.74272463156743e-07,3.31196104192811e-07,5.481872642976e-07,8.2692218882057e-07,1.16749109225162e-06,1.56902754100022e-06,2.03013177037545e-06 +0,4.63537374431423e-11,9.81049734780514e-10,5.2662674195959e-09,1.63744082880102e-08,3.80213424268806e-08,7.37013316252193e-08,1.26462911129735e-07,1.98832947513683e-07,2.92825213625396e-07,4.09988785684554e-07 +0,9.52762309519183e-13,4.03064693534399e-11,3.21954008431541e-10,1.31921085909251e-09,3.77652149951806e-09,8.65283546407388e-09,1.70473483728495e-08,3.01291907541136e-08,4.90801983841489e-08,7.50526597099419e-08 +0,1.69833113269482e-14,1.45989004470562e-12,1.76611013922202e-11,9.70786436099117e-11,3.4867494868619e-10,9.60502789076322e-10,2.20870660779088e-09,4.45777595611339e-09,8.15402734448277e-09,1.38148719622932e-08 +0,3.8915196963406e-11,2.71147134577767e-09,2.44901018232562e-08,1.01911652264688e-07,2.85170970090303e-07,6.29937137214719e-07,1.19197488155735e-06,2.02429502148315e-06,3.17588061786517e-06,4.69121110232959e-06 +0,2.69185075542407e-10,1.00640042548647e-08,6.5998653123948e-08,2.25483644613014e-07,5.52286794898488e-07,1.10885707763952e-06,1.95314930874714e-06,3.13753072469243e-06,4.70863852449393e-06,6.70761014036735e-06 +0,1.89155326006871e-09,4.15998208034506e-08,2.08604035348676e-07,6.01200476492204e-07,1.30548633543253e-06,2.39315754080736e-06,3.9233691392589e-06,5.9449170992954e-06,8.49806448732054e-06,1.16159780561656e-05 +0,1.20839005460873e-08,1.63063271190061e-07,6.45781047622635e-07,1.61034741128995e-06,3.16276730648069e-06,5.37778778885908e-06,8.30838167561027e-06,1.19919722977472e-05,1.64544997795969e-05,2.17132051523398e-05 +0,6.98490407129936e-08,6.0448889346693e-07,1.94055108474347e-06,4.2597649069398e-06,7.66233240335793e-06,1.22016422946016e-05,1.79016972891742e-05,2.476703497099e-05,3.27887852389906e-05,4.19486886713683e-05 +0,3.50437914406735e-07,2.02949943667367e-06,5.38075813382185e-06,1.04873870799375e-05,1.73426483579829e-05,2.58972599259919e-05,3.60796128969124e-05,4.78063686693464e-05,6.09885319751288e-05,7.55352438512354e-05 +0,2.52830836713648e-10,3.03681967057117e-09,1.23485455469166e-08,3.26088143094103e-08,6.82703234078585e-08,1.23666594902017e-07,2.02940810635938e-07,3.10010651321603e-07,4.48552679684649e-07,6.21997431145188e-07 +0,5.62803310698392e-13,1.27640446069762e-11,7.32343235940447e-11,2.42833323537223e-10,6.0008786533641e-10,1.23564347077632e-09,2.24827721329559e-09,3.74216806129839e-09,5.82509155188699e-09,8.60722506603327e-09 +0,1.01387223716053e-14,4.50063177569099e-13,3.76367908055613e-12,1.6114550692278e-11,4.81228068266938e-11,1.1484654272007e-10,2.3535141406549e-10,4.32106599374873e-10,7.30344360289206e-10,1.15746018360108e-09 +0,1.76594647095188e-16,1.5591711552107e-14,1.93446159551681e-13,1.08935554130944e-12,4.00524857658289e-12,1.12878863819327e-11,2.65435340079881e-11,5.47630083633411e-11,1.02366893517743e-10,1.77190911910303e-10 +0,2.72938929463367e-18,4.86658533444932e-16,9.10403238648642e-15,6.85236595138709e-14,3.15126909211129e-13,1.06497145870696e-12,2.91648358432319e-12,6.85836462119616e-12,1.4373187057492e-11,2.75299294530364e-11 +0,6.21517803140922e-14,9.59223157644613e-12,1.40071205391153e-10,8.23460517316941e-10,3.01588393968386e-09,8.29925719783238e-09,1.88974373890384e-08,3.76435597362723e-08,6.79236252804073e-08,1.13611867015149e-07 +0,4.89626154105836e-13,4.05815676062707e-11,4.27904064325323e-10,2.04932427478431e-09,6.51652855169281e-09,1.61738279532313e-08,3.40432642093828e-08,6.37469204528233e-08,1.09424110270049e-07,1.75651555512571e-07 +0,3.92394931909311e-12,1.90205223022698e-10,1.52475376523412e-09,6.12909774497683e-09,1.7202341627476e-08,3.88279099539024e-08,7.57896022082591e-08,1.33417356513312e-07,2.17461256827561e-07,3.33987690067168e-07 +0,2.8933367514824e-11,8.49634390328452e-10,5.32093344063423e-09,1.83496200775791e-08,4.62765938242325e-08,9.63901823877104e-08,1.76587279811171e-07,2.95134919493546e-07,4.60499444520641e-07,6.81219273509691e-07 +0,1.95394555097804e-10,3.60599662518081e-09,1.80378231091857e-08,5.41614780599373e-08,1.24057163023099e-07,2.4042528230335e-07,4.1611636054975e-07,6.638317374418e-07,9.95928543382885e-07,1.42428849426674e-06 +0,1.14722473747572e-09,1.37465796551839e-08,5.57323434828026e-08,1.46670400626854e-07,3.05907787720263e-07,5.51856057412464e-07,9.01663873101398e-07,1.37107045554111e-06,1.97435903337969e-06,2.72436530869665e-06 +0,6.01331068747146e-13,1.48666863208301e-11,9.22909151658988e-11,3.28914376451997e-10,8.6846657139908e-10,1.90065440926915e-09,3.65824807559149e-09,6.41362850871594e-09,1.04751195906183e-08,1.61832994397799e-08 +0,2.47209629421286e-16,1.20861421473939e-14,1.10842575020953e-13,5.18453630267307e-13,1.68538118789818e-12,4.36405036350932e-12,9.67329876973588e-12,1.91546745620659e-11,3.48213414272979e-11,5.92006053141907e-11 +0,1.89201707214526e-18,1.74461832628427e-16,2.25617452264583e-15,1.32210769025309e-14,5.05114773691188e-14,1.47735746305447e-13,3.60121598845708e-13,7.69378961925591e-13,1.48783448922978e-12,2.66185680931199e-12 +0,2.82938853517197e-20,5.166160471115e-18,9.88324205334377e-17,7.5994903425817e-16,3.56758223375618e-15,1.23000563514818e-14,3.43481298564478e-14,8.23311838199757e-14,1.75814766445246e-13,3.43039404178284e-13 +0,3.8648935971775e-22,1.41879457944149e-19,4.08084894427787e-18,4.18430606959238e-17,2.4523734952935e-16,1.01226704556591e-15,3.28784556445063e-15,8.97275887121034e-15,2.14634799098414e-14,4.63096283150802e-14 +0,8.69726640954885e-17,2.92087981800005e-14,6.81319014247032e-13,5.60977017733479e-12,2.67241290013728e-11,9.12044259353017e-11,2.49116450278525e-10,5.80809267726238e-10,1.20365994453422e-09,2.27799582449115e-09 +0,7.66381761027786e-16,1.38259280480535e-13,2.31866103567873e-12,1.54608384966719e-11,6.35500210965569e-11,1.94448402919345e-10,4.88264062451473e-10,1.06479889696641e-09,2.08990883233503e-09,3.78080528306937e-09 +0,6.87383174408679e-15,7.21096117919272e-13,9.14517774138443e-12,5.09570157318565e-11,1.84185808276166e-10,5.10862690038548e-10,1.1861963221612e-09,2.42559699182021e-09,4.50988933263573e-09,7.78908854017501e-09 +0,5.72300830893883e-14,3.59756333500574e-12,3.53239285399593e-11,1.6766501693894e-10,5.4153652527201e-10,1.38002614943133e-09,2.99692285950551e-09,5.80167590670959e-09,1.03017611203301e-08,1.71032265873292e-08 +0,4.4015162593732e-13,1.71041885299831e-11,1.32555658504776e-10,5.42962213883109e-10,1.58183574169237e-09,3.73024467366084e-09,7.61934524058875e-09,1.40282403312057e-08,2.38787854849993e-08,3.82284985005463e-08 +0,2.94015789946406e-12,7.24491356536544e-11,4.48309024211327e-10,1.59267035200818e-09,4.19217068372523e-09,9.14631653977701e-09,1.75503380261564e-08,3.06757992440402e-08,4.99505561914961e-08,7.69389888117878e-08 +0,1.25520662515885e-15,6.34202350482921e-14,5.98701785935928e-13,2.87279812512801e-12,9.55263541481702e-12,2.52380922965864e-11,5.69558750757137e-11,1.14607040894637e-10,2.11363987866348e-10,3.6401786096932e-10 +0,2.69293069884163e-19,2.76242819230984e-17,3.94910273521441e-16,2.5438021672209e-15,1.0629681164291e-14,3.38496697614216e-14,8.94675031954185e-14,2.06465258826618e-13,4.29763614969243e-13,8.24924625699432e-13 +0,3.40393634256781e-22,6.53632617070714e-20,1.31310268914671e-18,1.05893998255463e-17,5.20829887103865e-17,1.87956724398352e-16,5.48951609054929e-16,1.3751099157628e-15,3.06669713938502e-15,6.2447989892164e-15 +0,4.00278935798164e-24,1.50166583042231e-21,4.40866620105182e-20,4.60970828253607e-19,2.75310436445269e-18,1.15735040672427e-17,3.82669184362989e-17,1.06268988633135e-16,2.58592708871164e-16,5.67426085674954e-16 +0,4.89860423986417e-26,3.68055944030794e-23,1.6212071379999e-21,2.25639400694407e-20,1.6800782168471e-19,8.44489572108405e-19,3.24448865148392e-18,1.02486142861079e-17,2.79143092738667e-17,6.76923776023426e-17