From a86410fda7783dcd6c9ff2e58c7fa6cb8e17aeb7 Mon Sep 17 00:00:00 2001 From: Michelle Date: Mon, 29 Sep 2025 20:54:06 +0200 Subject: [PATCH] diffusion.h I'm done with the backbone of it, I haven't had feedback on it. --- bin/abc_lab | Bin 0 -> 22672 bytes include/core/global_config.h | 36 +++++ include/core/omp_config.h | 4 +- include/modules/fluids/diffusion1d.h | 129 +++++++++++++++++ include/modules/fluids/fluids.h | 5 + include/modules/mesh/mesh.h | 4 + include/modules/{grid1d.h => mesh/mesh1d.h} | 15 +- include/numerics/numerics.h | 8 +- include/utils/generators.h | 4 + include/utils/generators/linspace.h | 37 +++++ include/utils/utils.h | 1 + obj/main.d | 56 ++++++++ obj/main.o | Bin 0 -> 21184 bytes src/main.cpp | 148 ++++---------------- 14 files changed, 307 insertions(+), 140 deletions(-) create mode 100755 bin/abc_lab create mode 100644 include/core/global_config.h create mode 100644 include/modules/fluids/fluids.h create mode 100644 include/modules/mesh/mesh.h rename include/modules/{grid1d.h => mesh/mesh1d.h} (77%) create mode 100644 include/utils/generators.h create mode 100644 include/utils/generators/linspace.h create mode 100644 obj/main.d create mode 100644 obj/main.o diff --git a/bin/abc_lab b/bin/abc_lab new file mode 100755 index 0000000000000000000000000000000000000000..ac1a151b0e2aaeec911d8fa9b9c5f6ef40bf16ca GIT binary patch literal 22672 zcmeHveRxzwmiO(XJ82-KI|>>VhX#X7L}EG+AR?NOG~AY(5C}+Aply;4=`|lrf5f1p zM3T7Pr)}ijoe$=nFJ>OsadyX%omt=4`9Kqr4(OYO^~)TvWd=T=qstt}2^Nru6YDe;sk-6|Q7-+7frDChLrTAl$Y zmljA9@jXqlNm;;+jHZ>BYXGHMdM6_pnmL^+B)NQEb1Iy|`2-}JC=n#NbZP%Ytwcb{ zs^KJOqAbT#;h}G7d;*H{M193txs(D8OeyT*P%lTm8NNl68&;BGGjI1(qA}bqau>ZFV=4Hx>TY!WXg{{C%$=Q}y^#YeN&RTwWj-#+W; z&#FHl9uSSk|G=KaH@;MUGDA9b`?|HTLcmHqB0bj3j&Ey?9m8>Am z*Ds;yvN6hbf;JZYt2Fc_Y3Q{uek?t2rm6RfH1wO&=vkA7emISu#x!=8r{TXMjs82+ z@c%9i{k=5x=BD9qNyGm)Y3!MmM$aG9=%I0D!zVRmr_m4fWAfjahF$?Z(^HiL-;2*Q zDPNLAT;n{wfYUMhxATFV&FN88oQF}}RjU!KWC$gfYh!DBo68??`vNYPKeJas`@w?f1D_JZ&2T%`T76*Y2xb)*!j= zYHX?PXa%gPYcF$EyS#0Fk1ybCbU0SlwXZ^@chxMb3zRklHgEDc8ypK7+gLuFX96gw0Xj;4emx)gWDhIK$}U@-O|$D z;0{1(XM<-`z}w!&1WN4QcE4+b+wXB9cKO`C&9K+Yq&9YTx;A-y{`NL^i#M>@)lodA zeN}b`>fm*^`vX3YyH&Kg0T$UCn%zEEz~}Y`{AezmfEJ=cZ-a|-J6E|@6}udD^IWyd z>gKpwST*f)HU$HkOZ48d&w)poyO^DbCp5ZH=y#9)GYEJx$}G;a*om^S!Pnx3@)V^|ZFOcTfR*(!fVctJ@2gF<&-$$tCxD zh|%u}5NO!6Spu^^(AY3*76trygl=qrG+}na8Z`RS%1YOqqGD;O;#^$mDk++q$k|7; zb4IfzMWq@?eLdvn7R@P1Lhbe>v`nh3udi6_blHpM6zP#9u974nOpFXXGfen4;$y-y zA{)?*&p2tkGyyON(1;Qfo*NVKF)&U8N+=QDimZqawbyT+Ac3z}D`wme@N$}xle}K6 zxW-9gPW$fv`aT{@+0qX<{;v$rl77MQ;MLx7v{HE*I3)vk*Yh&M7mtu+j`Vv@b2EH` zw2$NJmO|vlOK)?0^J~w*KC|=@$KQD7UObXb(ia^6Hj~ekzUFwK@b9toffx9LTCB^& zdZ>%XpAo)Z6ZCRUH}g0o=mOuroUGTL$>Lc|_&YkfkRkd(9X*E#JVkZ%D+Cd6SVu?n zOr#MVoiEv0nOLt=xTJMQBJsFO>y~R0L@5zpiGHn)Ue5JXJQkssie()gu1%y`9bK%2 zsIXo~r*)1<_vq+Uw@A%8y0*^IN;c`};u%3DojN)OSR!@l=;ISqNqSUA#{f#C?K(OJ zb|UT2(Q^}2N!qESTXgiNbaa|aBJI-A_3NVNbo9#wGhm;NK3PZKr=wq~qhHEi^1y$_ z11BxtIxBa7VwOXuL$_gC_x1%c;``<9LGu7p6E7VPa(qS(zAaPBCE{NG1qnDvq74@^|ctAP4yB6+`UO>;w>IaH?Bzlsv z;Rm~-A5$i@`zg3P`aa7$_BG10j5B7U7lKp49{mAi!Mg*Rqkj}~zeVQwP;>%yYmsGY z54a_XjbA0as=zFFe|Ru*W=-GoGr_vnMn_znn-lb>q_WpRubyn2L`1u>Qo28vkQo`|9ob(ddwOB zi{;UUm{!*Y1t4mlhq1xn#dG+IHh@=pO2H( zx&3nJSyp#=i9zljF?1}EL#w;w?h6J>&v)T3S^dj%%*};w%iZyO%ino{S*B9w8TVj# zM9$S$?HP%jyR#n(js55$x$JcC&~}v1E~jeHa3%bY5-MgRkc=i93DQ&`cn}T3d_nW? z0fm~6>@bLwjl5xvK=so!45FvtdGtygD^a-cS@MJuewyUv@Uu+qIqPoX#6+_}m!ck0 z)Q|U?#@|e;NUzB=TBiYOzX=j!>4sN%j>4KrFSR!_~AjCg7&_)#-7&OpK>`hDe5<;}fOn+%&3dUc0 zI{}Puby)5{ahu#fk|`U8feJaGbkMWS?)h-Qo_cp?xVKAHHb21Ib3dYmh}Z}NdXD*D9?`Ki5`0hvXg@T9XC8gh8CSy6A0%lZR($XgZlZ!zlo14d-5k+WMcUzS<*jT?1Q z({{)oe;1u@Mr@8g{*^XB+0x_#6h*@<9vq9VMpG01)K)=G^;E^G)pfG^D)tNca%j0( zR!5xb0XejJrW{&fl|$nYaI*2aOzwTr((`+IexI}6jtFF}w<_U9d2=Pn8Cq{uLitMQ zit11`M(NyF$hmU3WxhFba_viz<7@Xvj;!y;=dTD{1Cg^2U3(VW;mDaa1DJZ5P*rY? zoCCwTyh<}o$z?#TOHs$#Y^<1Gik!Stwhr$3{AU@o6f&Uf#u82FjkS$0jW|PhkU6U* zr~0NG-=kQbIV4-2d95bgI4zF44?y+c;Ljh(Gadn#45Me9!6qGw_u%s~Eb8Z$7~f!( z+rZ|~;t4D;&R0}K8#z2XqDhUKV=N(ijCQQ5{t;!-?|q4()O{$+ z{D%*pL$`LnWK{3ShX=@G#y!!+P+PsQ&1~74k6_R0Rn((UjHO=mtCOg6;kPVXZ$z0A zf_XZI=-VhoP$(~6#wlP@ST?~Trd3enbo zKmqEFc2T7vbk4%n=D-ZH@?#?ZkdYl`=PU=FaD8KtokSW7*_pF z!m3Xo#jT3E$e8ZI{JYn|B(b81{Seed{L$9J>IpLBAPT4>oMjQqRzF-C8^Txe{QU!S zeYkW8afW1mO)})j>xloi7xUvmk8t@uc;$E``XH+UE(kfya^Zf~T&nFelC~gB(xbQ; zOxiNxRXmiFmkE!95j@%oVXKBj=O+Yya*@E7dVym=qE~_Y1qzc<*XK|cF3kr*aaSFP zE+R^3@{>1#5-xq=Mk@FqdL<+%aLZ~0D8&&E@nOj~5hbJHhdc!KDm#gmg-JlTSH zGO0ig&6h)KXozMYNLzY-kcc8@ADa9bR_X{O`%t?Qp1q#1(A@|;TV`UWu3#vl%Q!id zQ=o(layW;DZUmNn(1*ZrW(^iy>-%Xcj)oG#Dk&;Z6v3pYP%_x>zw>q$M+{Uv=)n1$(h&qhNnTY*UeZ|H%K&-N8b+O?~k;8)i zXO-G0cp1mBXkm`F4a1Gy3lS!E(6GgNRm>R|>Vb^Iek$rDrck9mK15KY7>1qkArTV3 zX?#G-czl=+waP+&a(tk9K;!hlXQ*@GEG<6pNyz5NpQs>Yno)$2ks2C?5|i+*l+f@! zoWMk1I{{IIhO3}LG_w;0s50tdYK<5kR(WBqIWQAW)S|-!C{aRl3)sY*m8B05eejvy zHk)RfI(rP;Zb;bn8F`xRRpZfH$)t+{L?@_XeZt~{7B{1P#Ep({VhQ>m#fK+x21N~f zjAY1>-z+|?13PPJ=n7e_!l{KWM#M0AV>vuVR54pnu>hr9Z;448@3(peDzN3gSbIhyfT544aN3=M(-?H^RgwtyE0!D|y zsea`wJf+Zb7SD{bJ(libsFeGU;Oymsv+#Yn`<&Ub?M581$m(|rie>e$iVAk@n)0d& z7DVU3DkW4^P)yUcV1YAKjz!lO*pVoOhn;FU&Wp}M{1rTQA%g*+iOGd;N6&o1wq7I@ zo;yV;J8J1Uk0>U0A4kmitK9!Fo?v@WyS1ZUR$rsl;XXw@DEG&(M?GY^J$%cw?nAb~ z70z%~fo-t=!MvOTTE^#V>w~2AA#GAmoZ5pu>a(nZ!{@8RH@=DZ>kK<+b(Dt(CT)iv zL@7vCM$&dje+v>rtuYs}3I7iyJcN4&P*yM06 zf@L0B@^;{IAP7%+a6d2;<*ZJ%)66!l+q%#Y1+xNOMqP;RqJF{CXyfsezU(wBhSMZi zU3SdU(}_WYT{JoYeJh837;|4Fx(TcK8qprn81&Kg)R@as`e+i?*ZTGQX#Y$`6>_)M zKck)FEZb}8Ig3#_+BuAL3?!{{V14WbV2S=wLRe&u$j@V%Q?e+9+W0OUlFQz*Y<-uF z|3lSg8x8Mjj&! zU-?coLQ}zZdw4GPEtT*Oo}<=rN_f(@5`G%7cHAkeZ&SbbvpuAC>ak%WJgx|ImQ~7t`O@b*|2}Kv*og3*|KyZ_3#TZZI?MZqi?`=P*hWc6b#=?+jJV&?!?K}M!$(rq+W!ne@503CW( z6<{b}e*$9=KJ2Gz!bVJ{y)=l-(MnL+kid$R9@HvMkRQNW)o&A9bVaSiw5znyi5ec{ z7ozt-gckmZjZRr4xV$=KET{-sdtVPus|b}Fd;2UsGR#0ehw|rP3XY>KJ-?%PjGE-c zZy}c6yIIs3tiOAm*xU>#VLa>TOhgVV1vAl7J3aH{y;TMD^NH4~VnWLajkY3opcRmD zh8|@4%blShVlF*Q?%x7ZFP;VU7-Q6cpM8WarLfo1+X`}QDZX_5Z+{WTS20>pXwS1$ zlb=G{u;Qv$)X#Au%{E5IM$M1y#vV_nt!aCzeDl9AKO0NHB>8H@H?uQThC_RU<&h_0 zSLEzNGaf;YJ44^LvV*|b6YzwXPc(AT@7Pm(N}B-Me?-n%rwlpO*R=zM12lB71Bg@` z-gX-HM2@Wex50<8^+zOmz`6>O$JhTEp97JT^)E#}U-L4J@*H-W0JXIHh3o*6c%wSB zVkWC^JUz_U8pl_Mt{@y+5eSOoxU0N}IG&(`acGQIdvNU~}|^i-^JFbhZSs=p2-3 z>%D>F^5buKGYiME5wtEr#{gzY z;2NA>uEAmaOM%is4rkEsMG}25BpEbZlesW}Joc4R@dxqvQ-JpzjK}u_ZUT$~Za)-{ z-+~BWJsgh@LvRzXRBX;PNuyP6VkbklT*~`WC7f z3CFE08SY;x89J?oYjU#9+rdk8{J<&YqmJu|qbi$@UZetia9%2Py&sR$zD3HlF3HWm z-7?`mbC-16)nzx#DVPo_lV1z?DB8jViHJmc1Nd}-ZX`hGgLzG`M4J%7G~ z!>b$UKPWZSo+T5roR4R2Glor%Wp!s~HW<*BfuWaQdG#-~Ia${3>_;*l%L<#e86VHg zT$M=;RVJ(;uzW%nRM7ZXiytfU$iBH3+qXEmy;Zr}4Vjl45)Fn=5)H1Hknt1>s+lMB z1uG}ir=g>M*g#`4ALHOBXp@l^FV39t(rRl~Mu&8}C9|VG%WBE!=(Of$bZoLt%;;#g z=45o-BU>{wI%+c?MU%!$)o4$}1R9U|=w=0DGKjqKPx#Ul5HC~6WHl{JA87`c6ptGw z-2Il}Uoxbj46NXQ=Dn5ibL6&}O}{crPnf4_#HTZU3gY2x)7RP3NVZ0tww2zlc|XJS zlWgfkrbCim%rZThEgj1;h0M~9Y(JB{sTMpZ4I#s8Inq|6Y1>5UfXVbpj`VJpY3oGk zv2n9NSI3(U=1AQ+PUK;MaTP=^rAr>Tw~@McN7c*AKw zK}qZx@YgApY`o!uj#rFW5_XHdhj=GL)X!hHPP7k)d#rt^j)~$RhT&G;kT}PMo%nH* zmG9?@vDIR@DA%m`@B0?-VW(;==KWpHVJ(OEaJY%XE)KVIxRb+O9QJXzpTmP34s(c? zAX&0J51RF=(^T!5eJFv*!gF zTJliI+@j*Sg`Au!nT%7tM(Oi-GL_yqqd94m9^)Nu;|%z-(3t!*pE4w?)P?`mP@??{ zrR(t-OTQg-;@>YOE-@p~p$E|=F<}`M+fGhTjVrV_H(-aqeTBv-{MgLsS4v{OYo(I( zI5N!66#YM9{8vaKehB_QATyTz>Du)o(=$nGq|KZzrPD(!W9fgJ$zCpr_xFXWkC8DT z+NI8ylZ<|KBEAUzac~Il5R5JGfzJF;plN64RLFjp(>K{PI-Q46 z%3*%c^7`*$(8qeOehKI_?^5T{e2h;k;@A!z$Aqpr(lf?<+z2|^Y2ynSI`5-&h|^R1 z>ny8x^gVx#8|xa(r?LD{l7@a~8hUFQdN>XJ3D9lA2p#N8!~Ztu)ZV@uG-f*Qqcj~0 zBOA*lf7#Q}mw`Uk{RmQBgD+q&YR7x)?hRg7z`aqT*Vv01Hf`dW_LfEoJpMosKX5fj zqjw>=03yT#RrGk?F!9qg1E+FLiZcmkeAyvjZ%2Ywc4^SazV zpL?^*(-!b;mYRHS+-lGmY;D~P5go?fVT&^W6E2@6t;y zPoq2FhC*JSt0~yl;D@=5xc^}}-u!R$wz-0S574w%Gw2l-;l+AzmIedf7XQ3E*iR(5 z`JuAJ(Fzqwmj`HKc>kZ?@7G-}aPdVA<)d%(vzral-aVnWgx^h&R13Pv-|oUsCyjK! zMWW3eZgeR%xoS<3&(l^V8yKouDkB4!F>-t&-}{vwkGdJKYYHR(bMD(wgmi+4u^dsY17{$ zanY?C^RyqBoQ-vEx~hZUYJtlnlH00X?xAa5Q5mVe=r#_LPnJx%WFlcctcNS;rV3m? zA_n+>=XDx-Gc;@Uj=1PTjU;wlF2P!t;vM_k zt}SWAkBL4{TZ6|{;%fIHPPuA4!B+fE=vXPHTeUmj^LA42t*oglag9ztN8JJn0b^c| zqHiwaVfTh4Q(Y7!v}-_;r;QlQch#h8(%k=4ld8NvZ$mSp>Njr`hF(hNRF%*@D6S2D zzc3jsPj}CYzM*up2=(?@x7H-_mC+QW8&WhXk6~kl4c$b6Nr@Y5lF0ww9Xv@FwUmpS zQ($Boe=l(R{hn3?m~ZUBF)jp3bLr0*ckbxDB1uL5&8>)TfB~Ndn?=sshN-YgD#EYY zo+3w~>=@hs$hwU^Dh+x4&5`YTVq0Dl{1IX(b)Fp@%N> zVQd_6q0HxLagzYgZE6WfMJ&z~AzQSu9Vlh|o(8EX;OPX)B65+hojnnYJk9(8-P{Om z0z)+o4g3lrKKPkd?b)sA5P8x=WH{A3BQ?@(ZxOrX--)PiWN9^EH9vCAcK{W-c#afnem2FY^*q zXNtTypAgWN;M6jy@@qiG|F7GH{o=e&K(TGda+4*ozvHOnMuv{yguFNp6mUN;ARm#R zg}gYI^MXcA7V_eJQ9#iTU}1@L=+n1>&{>_37w2;V-olxM{en-xhfq#ue*zche*zA3 zd9t7E6!vqn6h?-8D&)m^q<}VVu&7_ii}wG3%P-~%#rdXy;yje(Y0H@^|8vl&3ZWlc z0xfY!<6ECn`~Sx%)XUqroPg}VAUXvp>ZM<%$k*m+%mNC2p+evS?n;rrhsz5n{FhGt z_bKvSeDfyYij)Z~>J|E4NRb!k=K_jzb-|zNzX+F4^nWquDbXdv+|1vQt(-jCtbpR zac^OLxyDHETIy4|9zmYwH?b6_(1rZ06uMBShhpsLSsyLeFdKIq{#hf@Z$P5mqF?A- lUD&QkXx~pQ)nKb|JTDaW3SJ5qW64imu5py52&ACW{|6GMSZx3R literal 0 HcmV?d00001 diff --git a/include/core/global_config.h b/include/core/global_config.h new file mode 100644 index 0000000..c18e981 --- /dev/null +++ b/include/core/global_config.h @@ -0,0 +1,36 @@ +#pragma once + + +namespace core{ + + enum class GridKind { Uniform, NonUniform }; + enum class FDKind { Central, Forward, Backward }; + enum class BCKind { Dirichlet, Neumann /*, Robin*/ }; + enum class SolverKind { LU, Inverse /*, CG*/ }; + + template + struct BC { + FDKind fd{FDKind::Forward}; + BCKind kind{BCKind::Dirichlet}; + T value{T(0)}; + }; + + + // Global default config holder + template + struct Configs { + GridKind grid{GridKind::Uniform}; + FDKind fd{FDKind::Central}; + BC left{FDKind::Forward, BCKind::Dirichlet, T(0) }; + BC right{FDKind::Backward, BCKind::Dirichlet, T(0) }; + SolverKind solver{SolverKind::LU}; + + static Configs& defaults() { + static Configs g{}; // process-wide defaults + return g; + } + }; + + + +} // namespace core diff --git a/include/core/omp_config.h b/include/core/omp_config.h index 91b7e75..6449035 100644 --- a/include/core/omp_config.h +++ b/include/core/omp_config.h @@ -1,6 +1,6 @@ - - #pragma once + +#include #include diff --git a/include/modules/fluids/diffusion1d.h b/include/modules/fluids/diffusion1d.h index e69de29..221af99 100644 --- a/include/modules/fluids/diffusion1d.h +++ b/include/modules/fluids/diffusion1d.h @@ -0,0 +1,129 @@ +#pragma once + +#include "modules/mesh/mesh1d.h" +#include "core/global_config.h" +#include "utils/matrix.h" +#include "utils/vector.h" + + + + +namespace fluids { + + template + struct Diffusion1D{ + const core::Configs& cfg; + const mesh::Mesh1D& mesh; + T Gamma{1}; + + + + // Constructor + Diffusion1D(const core::Configs& configs, const mesh::Mesh1D& Mesh, T Gamma_const=T(1)): cfg(configs), mesh(Mesh), Gamma(Gamma_const) {} + + + + + + + void assemble(utils::Matrix& A, utils::Vector& b, utils::Vector& s){ + + uint64_t N = mesh.center_idx + 1; + + if (N < 3){ + throw std::runtime_error("Diffusion1D: need N>=3"); + } + if (A.rows() != N || A.cols() != N){ + A = utils::Matrix(N, N, T(0)); + } + if (b.size() != N){ + b = utils::Vector(N, T(0)); + } + + + if (cfg.grid == core::GridKind::Uniform){ + // Core of A + if (cfg.fd == core::FDKind::Central){ + uniform_central_finite_diffrence_2_order(A,b,s); + } + + + + // Left BC of A + if (cfg.left.kind == core::BCKind::Dirichlet){ + BC_uniform_backward_finite_diffrence_2_order_Dirichlet(A,b,s); + }else if (cfg.left.kind == core::BCKind::Neumann){ + BC_uniform_backward_finite_diffrence_2_order_Neumann(A,b,s); + } + } + } + + + void uniform_central_finite_diffrence_2_order(utils::Matrix& A, utils::Vector& b, utils::Vector& s){ + T xm, xc, xp; + + for (uint64_t i = 1; i < mesh.center_idx; ++i){ + xm = mesh.center(i-1); + xc = mesh.center(i); + xp = mesh.center(i+1); + + A(i, i-1) = Gamma/(xc - xm); + A(i, i) = -((Gamma/(xp - xc)) + (Gamma/(xc - xm))); + A(i, i+1) = Gamma/(xp - xc); + + b[i] = -s[i]*mesh.dx(i); + } + } + + void BC_uniform_backward_finite_diffrence_2_order_Dirichlet(utils::Matrix& A, utils::Vector& b, utils::Vector& s){ + T xm; + T xw = mesh.vertice(0); + T xc = mesh.center(0); + T xp = mesh.center(1); + T xe; + uint64_t N = mesh.center_idx; + + A(0, 0) = -((Gamma/(xp - xc)) + (Gamma/(xc - xw))); + A(0, 1) = Gamma/(xp - xc); + + b[0] = -s[0]*mesh.dx(0) - Gamma*(cfg.left.value/(xc - xw)); + + xm = mesh.center(N-1); + xc = mesh.center(N); + xe = mesh.vertice(N+1); + + A(N, N-1) = Gamma/(xc - xm); + A(N, N) = -((Gamma/(xe - xc)) + (Gamma/(xc - xm))); + + b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value/(xe - xc)); + A.print(); + b.print(); + } + + void BC_uniform_backward_finite_diffrence_2_order_Neumann(utils::Matrix& A, utils::Vector& b, utils::Vector& s){ + T xm; + T xc = mesh.center(0); + T xp = mesh.center(1); + + uint64_t N = mesh.center_idx; + + A(0, 0) = -Gamma/(xp - xc); + A(0, 1) = Gamma/(xp - xc); + + b[0] = -s[0]*mesh.dx(0) - (Gamma*cfg.left.value); + + xm = mesh.center(N-1); + xc = mesh.center(N); + + A(N, N-1) = Gamma/(xc - xm); + A(N, N) = -Gamma/(xc - xm); + + b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value); + + A.print(); + b.print(); + } + + }; + +} // end namespace fluids \ No newline at end of file diff --git a/include/modules/fluids/fluids.h b/include/modules/fluids/fluids.h new file mode 100644 index 0000000..066af5f --- /dev/null +++ b/include/modules/fluids/fluids.h @@ -0,0 +1,5 @@ +#pragma once + +#include "modules/fluids/diffusion1d.h" + + diff --git a/include/modules/mesh/mesh.h b/include/modules/mesh/mesh.h new file mode 100644 index 0000000..369ebd7 --- /dev/null +++ b/include/modules/mesh/mesh.h @@ -0,0 +1,4 @@ +#pragma once + +#include "modules/mesh/mesh1d.h" + diff --git a/include/modules/grid1d.h b/include/modules/mesh/mesh1d.h similarity index 77% rename from include/modules/grid1d.h rename to include/modules/mesh/mesh1d.h index 8fc35fc..5373580 100644 --- a/include/modules/grid1d.h +++ b/include/modules/mesh/mesh1d.h @@ -2,20 +2,18 @@ #include "utils/vector.h" - -namespace fvm { +namespace mesh { template - - struct Grid1D{ + struct Mesh1D{ uint64_t center_idx; // max cell index uint64_t vertices_idx; // max vertice index utils::Vector centers; // size N (unknowns at cell centers) utils::Vector vertices; // size N+1 (face positions) - Grid1D() = default; + Mesh1D() = default; - explicit Grid1D(const utils::Vector& midpoints){ + explicit Mesh1D(const utils::Vector& midpoints){ centers = midpoints; center_idx = centers.size()-1; vertices_idx = centers.size(); @@ -33,11 +31,12 @@ namespace fvm { T dx(uint64_t i) const { check(i); return vertices[i+1] - vertices[i]; } T center(uint64_t i) const { check(i); return centers[i]; } + T vertice(uint64_t i) const {; return vertices[i]; } void check(uint64_t i) const { - if (i > center_idx) throw std::runtime_error("Grid1D: cell index out of range"); + if (i > center_idx) throw std::runtime_error("Mesh1D: cell index out of range"); } }; -} \ No newline at end of file +} // end namespace mesh \ No newline at end of file diff --git a/include/numerics/numerics.h b/include/numerics/numerics.h index e3cd5ad..bb4760e 100644 --- a/include/numerics/numerics.h +++ b/include/numerics/numerics.h @@ -10,10 +10,6 @@ #include "./numerics/min.h" #include "./numerics/max.h" #include "./numerics/abs.h" -#include "./numerics/interpolation1d_base.h" // base -#include "./numerics/interpolation1d/interpolation1d_linear.h" // derived -#include "./numerics/interpolation1d/interpolation1d_polynomial.h" // derived -#include "./numerics/interpolation1d/interpolation1d_cubic_spline.h" // derived -#include "./numerics/interpolation1d/interpolation1d_rational.h" // derived -#include "./numerics/interpolation1d/interpolation1d_barycentric.h" // derived +#include "./numerics/interpolation1d.h" // base + diff --git a/include/utils/generators.h b/include/utils/generators.h new file mode 100644 index 0000000..3369d29 --- /dev/null +++ b/include/utils/generators.h @@ -0,0 +1,4 @@ +// "./utils/generators.h" +#pragma once + +#include "./utils/generators/linspace.h" diff --git a/include/utils/generators/linspace.h b/include/utils/generators/linspace.h new file mode 100644 index 0000000..9eb5086 --- /dev/null +++ b/include/utils/generators/linspace.h @@ -0,0 +1,37 @@ +#pragma once + +#include "utils/vector.h" + + +namespace utils{ + + template + void inplace_linspace(utils::Vector& a, T start, T stop, bool endpoint=true){ + + uint64_t N = a.size(); + T step; + + if (endpoint){ + step = (stop - start) / static_cast(N - 1); + }else{ + step = (stop - start) / static_cast(N); + } + + for (uint64_t i = 0; i < N; ++i){ + a[i] = start + (step*static_cast(i)); + } + } + + template + utils::Vector linspace(T start, T stop, uint64_t N, bool endpoint=true){ + + utils::Vector a(N); + + inplace_linspace(a, start, stop, endpoint); + + return a; + } + + + +} // end namespace utils \ No newline at end of file diff --git a/include/utils/utils.h b/include/utils/utils.h index cfd11c5..86904a4 100644 --- a/include/utils/utils.h +++ b/include/utils/utils.h @@ -3,3 +3,4 @@ #include "./utils/vector.h" #include "./utils/matrix.h" +#include "./utils/generators.h" \ No newline at end of file diff --git a/obj/main.d b/obj/main.d new file mode 100644 index 0000000..be73496 --- /dev/null +++ b/obj/main.d @@ -0,0 +1,56 @@ +obj/main.o: src/main.cpp include/./core/omp_config.h \ + include/./utils/utils.h include/./utils/vector.h \ + include/./utils/matrix.h include/./utils/generators.h \ + include/./utils/generators/linspace.h include/utils/vector.h \ + include/./numerics/numerics.h include/./numerics/initializers/eye.h \ + include/./numerics/matequal.h include/./numerics/abs.h \ + include/./numerics/transpose.h include/./numerics/inverse.h \ + include/./numerics/inverse/inverse_gauss_jordan.h \ + include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \ + include/./numerics/matmul.h include/./numerics/matvec.h \ + include/./numerics/min.h include/./numerics/max.h \ + include/./numerics/interpolation1d.h \ + include/./numerics/interpolation1d/interpolation1d_barycentric.h \ + include/./numerics/interpolation1d/interpolation1d_base.h \ + include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \ + include/./numerics/interpolation1d/interpolation1d_linear.h \ + include/./numerics/interpolation1d/interpolation1d_polynomial.h \ + include/./numerics/interpolation1d/interpolation1d_rational.h \ + include/./decomp/decomp.h include/./modules/mesh/mesh.h \ + include/modules/mesh/mesh1d.h include/modules/fluids/fluids.h \ + include/modules/fluids/diffusion1d.h include/core/global_config.h \ + include/utils/matrix.h +include/./core/omp_config.h: +include/./utils/utils.h: +include/./utils/vector.h: +include/./utils/matrix.h: +include/./utils/generators.h: +include/./utils/generators/linspace.h: +include/utils/vector.h: +include/./numerics/numerics.h: +include/./numerics/initializers/eye.h: +include/./numerics/matequal.h: +include/./numerics/abs.h: +include/./numerics/transpose.h: +include/./numerics/inverse.h: +include/./numerics/inverse/inverse_gauss_jordan.h: +include/./numerics/inverse/inverse_lu.h: +include/./decomp/lu.h: +include/./numerics/matmul.h: +include/./numerics/matvec.h: +include/./numerics/min.h: +include/./numerics/max.h: +include/./numerics/interpolation1d.h: +include/./numerics/interpolation1d/interpolation1d_barycentric.h: +include/./numerics/interpolation1d/interpolation1d_base.h: +include/./numerics/interpolation1d/interpolation1d_cubic_spline.h: +include/./numerics/interpolation1d/interpolation1d_linear.h: +include/./numerics/interpolation1d/interpolation1d_polynomial.h: +include/./numerics/interpolation1d/interpolation1d_rational.h: +include/./decomp/decomp.h: +include/./modules/mesh/mesh.h: +include/modules/mesh/mesh1d.h: +include/modules/fluids/fluids.h: +include/modules/fluids/diffusion1d.h: +include/core/global_config.h: +include/utils/matrix.h: diff --git a/obj/main.o b/obj/main.o new file mode 100644 index 0000000000000000000000000000000000000000..15bd2270574a5600f23fe359c812761fda8c12c1 GIT binary patch literal 21184 zcmc&+4R};nnZA=uet2KVk3^XWdYs~PQlv)##fpg)8QUWC{&@m*F5}O~KnS{{V zf=p7q-cHLtyRP+Fb*t;Hb@%D&y4FX1mYR?Zv{)-`rPg*;3NDaR`H>dFPnmtcbI+MM zxtXcO)yMNZcg{KQ`Mz_$bI$jD=jYyJYrt5X<8tL{BDl2MG?&aM@wPxZQ^S)r+B|J8 zr@aDs;wi@C5$O`dSK_%!S)jJu-Cep+xPf#JLlLjV$)N1Bvrn3)r3|ULs`!XL_z-@zBw=wd-ov zuODu!@J#FHbujZJGcT~Aqqlkn_5(5EW$#ddVIDb3u^x|>jV+BA`9|J2*R_wCyO?=) zIW3j;+l=mF3IkP7>%2@2kMHavMd}3@fi!UlP%o|ztP_zVT zuVH?wkTPy0yy#^vGnaYQ&ZQe@+X4P#ue12asV}@kG~l{9!pz^Ih$FJA$)sz_+E-lv znWOE_$s>9S4WY%B#7~pgutnCbd%T8KYFO7cSPdGB&;1!S!{Qwa3KOT-?M)n8w=Z$%o*_K1B`!1# zC(b`aN~tBztr<>SD9J@tHQvMp5L_s!D;ybmz`G1t?-I1Jt{`#xp1q0FcQNnCu2cV% zgTfLQ%5I)*Grhd7b?=yN-AUrC)^zh1EV;|@eCG}3`Ob^WINnA0%)uFwCSMv-45Gl@`Z6o8)o(S$rOq z*{T@3P4o1B6IsLTC@XB^_QC5CRH!ZB|GiCXEiNRRYSOLiAn`DY4A)UTnP67z@9dVA%FR6?%3DQl4ibbu=iy zXK_A}HmkXZ$Bw?*jarzTg{uprGfB!f!7JcEp?*_9V?$#w(q~!qQP02-ntZ=feIzB- zTi_B>*}}v|O4qK3n<7tAdrxhU`zlO^?N5H3Lme}Zrv+JDs8=UE+anm#(cSnNZ~s5$ zLvPP=iiZaB2BQDoM0&ga-+-1S6X_wa0)?!6AMY-;c7nJ)$z&3CZi{W@viLQSO?zfC zNodtnvg&=F{>MNFn&dvuGbcy4+nM)dcw1#Sf()nr3AtID&+ZHPansfdX7P4~l1hQN z<2MtEA7bVSUIYA4$tyWuhq^pknD^~M?VE+uNNzi%{R_N_?O#sY{$(rt%hhGfTEMI| zq_R11`kwwTrM$}dhqxHRu^6_|8P-}jqOG%p?j??RmKa4S}5z_95OM!77 z$MXtYT18R}k!-_#jI;QEuQu#{8DPD~Zi}lkc@huwfj8MdJ%lkG2rQKoDr zo2%HJGRnLE6j#A4H`K$+7Op6a&cYz3y@g@TEz{vK^_%ilhcSpTd%v>Xn3C-TjkN8? z#Q0}&7~}kf?MG7+p;zy|EPufV@^btIdTL$E%=)r%K4TQdJX-WVe8$773)o_3gz%Z2 z%qf-OcgS$+Kex|#D(7qP83zmU-Yw9^3WU$#Zf4*mPk#gxggwa^?G^d{q;^%KZobXT zmua&I7gA$ba3R(buTef``_X-#?Z1Y*Z7?rETU@&Nj$VGoz^#_ihGm&-&pd3ivGQMXPQetPJJqN@?CHM;6q+xLvLb*+lAJ>y-mz~ zk#?WY8s-5uG@8SP-pE@Tzjb=-4PW$HJziht8yR}Aq^OKlk3~!E{lmEZC-9!*XLjX! z`k&?v9K6^Nzxfx5vCVosK)a_BEU3`U{s&PC$jc$@MFD*Z5YKu?gGi>~D2-vB$0T5GSv(A#`$KA)a9UU{ z`y)GpIr=bl=1RvLU5)*_dOk?c&oyX*#q9Y>PmW%_$J2ivTA7|44jlnxPY$FX{RwiZ z`7tbPLdV!MFnbCfqukc-^4+ZZ70>q9x&FV=P=1Pb3)`O|jsMA*O}iD_ajxrCqye#m zRAzmqtcsaNSt*U1&p@h@W?5)IZoW&hE4c9PxC+{Kcm|eSNMdX9 zfVZfORPgq((fI&U`qwj=%8IE7dzi7&qvyv}0RC@Wl>C?^&*QJ^2jhujpr( z`6jlE`^mBKdy5Z%htV4M3DUjzt^28Qu2E1ed#!~#@}F5Acf&aMkW$)PL#|@j@zO@X z9Aef3*b7H|@=C!l>tM!pJ~ohyVxw@0>#>f5M(SWyo#?r}!E%??THb+|dZyP}HSU2y zPd|egNEcE1B1DOexjg;u`y@-KFZ1He@_RaD(D^2D!8>)gZoXt+ zVC*N=!3iVL;QFQR+Lbu8?pGt9#Z?15(*xdB;5>HEPx0(eoNn5iIJM?^vYH}(dx2_o zoVLL`%sF0euL! zCKq@Ic5@b*5PpAW#8<%~&=u;rHG}-Rpf(vEvK* zVYLH;5P_fm2U5Z$j;?vhu%6;2*1)qMFv^eeOUQxmBBKvo#2)&Hp8UDz^9!-v8yPxw z;~wY*-BZKZb`lLv;rBf2>0gS0L&-cU@ee}1p8g>U75g1|c^dFvAz z)XS(o5I$tscN6?fmpTKV<(>DxZreN8=U-z2Y*?@UZS-pUPDkF+*af)sGO7=S_hqVU z{P{IMb=CDM`et8sO@qg_&xoQnal@7BEF56s($>H>@`Ui@fhRlztHA?ThXKH#*jG6i z?M+Y!wr8LQ5f#wV@n*pPoc~Pif0J(*IrQR;F&9Sa+}52tFd)nfd>(1e20buL>`&J} z&f;;Zz$o9F(i5xU@ze}h1AWwr#YY`uK8cxgv06OM+m0`GF*x|%#mrj04VU3?^?y%p z-TXOsPB_58b(MNP*A8)LKS6hiP!HUV*6$jVK4OSFfU`Ys>_EAF8;Xt_@qf|KFxNm5 znq;JyjDX?;=nyOg-Tx#I>OPVu0a7yYvKPLk_h@Zuu=G9*C z?j%Z>m=ARs<{R+tdpLt#Z>FSDZramg-*G~=r>PcZqmW`}C9+%!dD))Drfo^3VKOU= zL3;+CpfUdJ2$`UNaBc0~asgeh$2W`(V@sRFtYAM;3p(0EU&Fg$p8a7V5^^)mVtWb| z+oAiR#Jk?IoW!r2>T|{t@7}cwTWL4uchE;ik3%(QcH?Y#7;WM$2!BISeu>4GfPki* z4gk=p#5r&4dF-yQCRBdqdDeU&p)*=Yjlp|WI zEgX%brG|mudHRoFSPbiFay7Ia7l%AtVe+G3e=+sh`l4NoBjs;}Ne(ST z7)U?CU_M6N)`|E|ZYLIZb+r8-qo1?PT)%B(__-M)u8~7SGmR)(biz#C_#^UuZ`*iLn&*(HS1*v&i#r4H)R|0|&Qce{ys>Mu)H^h{;g~V9q7l=sPr>mb%{L@!{b=VvQv> zQbCbTrLhVKOH?Uz(|!d=)S!?>SUW$q0;77`Ww~A~(zd(v9=$@_pO^Pmk@k9i-u5fB zM<(3_+?gsOs&04*orz_gy^UtZ6Q&IWR`8_;4`=Qx$D=O!f3z%1lr9-1TcIwKK z<*TW_LBg;o@Nl(=UH7llTzy{G^+g4RJJ2BU(7Rd+1$M--N|Bx}VoBq=e}%zoDELPz zxn9;o^ksr}8S16BwLr1=;hefF@)yN(S3Kg5<*l8JR+ml5d5PlXQw)$a37$GIaGu4* z-p6tl6_-AmyQtU~cL$1RJ(5>jT={T*ZE;nsU`g?c+@Ix4F0KMpTRbaJ>;vhdV(+42 zQQwoIzU$t9eP7F5Ujvt~N=(X)`KWj~^=I)F`TAqIkGkV|kL1S+a$8*J%kb{!fA+%9 zR}|%YV+9ZAJdz*Jd({0{ZtkjF>Zmd04vKD{vJF+x*@r?&7Wo#9?{9ta4p;6~u2hFH zCaDhBPRV%^1r2;8)Pj{$nzG<%A8kJ`WcVuj#K&AOnyjt%=I8WkOFg;0P5E9=PH&&L zIHz}$_llg}jozZ1-cK@bZcgut+|Q#+leGr)r*=ve>M9k6(StOv%dUWEd4-(Uu7%e} zHG?#3Ako~~DZWHIiE=f*<8#8$ar#nWV#w9#tA&Xnmw&~O38$|bCWc)8l|v?+zIvD# za`{&fnQ&=EEU;F;=QngyqOkWBO+%YttPo=F~B<7Q}SyQX!9!ek6_ruwhWf^W!z zKad6gTo(L`S@0iZ!C%ONAIXB}Kt(e3t27Jl&w?+>f;VQtHv*rbrTrnTzZ4E|Jni@R zItsiJ<(cHH%7T9a_zcZQRq{aVCxtI$q5s<~_&0!GsrX4+11Wq9_}8_xHCEcNrgb*A zcU81(+9VR)9j#i0QCCSnv*~9p{mi4EIrKB1e*E;afPSiYmgoJP++WQRPV48qeopAG z(t>v{Z;Z~jZ;kcV#%8@W5U8IQ=%s>rJ<;}#$oxC`$8IR7n-}P8ML}b9w!ghQ65P-n z2?aaaH-wwRTkr{7dz5o6YmCloiEh~x(pv&mt=+-R?X97%KuhM9X7iR>J67~|W+w1w z5oI=sQop7m9BQiw2Jh|b3vLR9Bi&ui9qrLA!QM*L*0q`Hgs_#tV02@+dvmZO)OBxk zV=xpBcZXLjgGj5pHskZL;L1>>r!zzXOz5)T?(*{aZ5=)Btr5Rpxsk_^%$vKYF4)u6 z-qs!N#CWvax4Ajoif;lj7@=S*F7Lyku9i@6cCb6#8VU#N+r#ZG8#_YLz)IOvLvu9T z-bcf}a(UzIV0yFzja5M||J)XQ+{RP>`N3d!B#P_)PVvd1-qINL`&%|Phl9~@b9*#` zK@Y41HV2c_-Xf{>Rl!x2!9e5u;EH99bAlbz)yD2>e7|R9;qF~OZo~W3y z?s@2VWTQO`>eTw@fSb=PpL$M!b|IlX)*G5zgDuSw(urV8Uvsdzqocb8uPK6|zLwCY zC|+sswjgO3r;YYb6wtKP`2&1rmKe7$sx^YQFJ9ajfbg}8jDRhilcC<$kX1T$dR$ea zVJ?nHV>%oHWSiI7mo=};Lp_~%7xVsSb)wPv=VCmSk4Xc|tD)CjJs7+`|AJt!3mVpp zFHNJNa5KI=)mOAO`h$T$=gMWRT4$)Uv%5D$am%JHG_F~?Rn;7cggQ5Lgx;^Ds4)_P z@&zMND0i@B<9)%l=JpQmm(b4L%b}d%(1N<|uD14jBT%EN)=*n>Pe(Ko=neSq75CB8(#>6>4P zFGHFF)lWWG;@b&CAe_ECmbg456MmI~&sXH6-J+ylso-?sB=I$hJbnf5%R(Q^LjRhA z(_U59e^$ZO`su|y1sWd}pDJ*Q)Nxp*;A+2CDD}(%k8HO|!RIP?t5T0T4jl@;Do;$o z)qefEf~)l3Q|edcrw?Q(P=8f9i?A=C;2ei4ffF0;+GM|Kve0*D!T(;t=PC7kO{u?3 z!4p~Nf1ZW@v_d~$spnkY+1`tY?9Oe_X+7KShBe)gDhP^QGz!@lv2-)gJ3PI?;YryH@S=ZiQae=K~6^>Ssk3 zJCgS(B#%npn1%k|6kOHkJqoVs^MNcleJVqN+Ew-ZS}|@Esr|ZM;55#v{@NtN) zsYmVC0)<|c^9}`9`}G+GSLvl)Q?V-le<<{-od2WXYCZguaMVG?s-1cj{a5Kfkp-u7 z0|go%)lS0#r~1qANPYft7W#kALcd3$SMBsg1y}9#Riz%)P7f*cs-30eMJm;3WTfje3_#*UWc^o@AC?-wi{IJnl9lb zeO$p+`a8w=P^9`jxlWNhs{TklCwf(XmMY_^>d&1DuGX_&;8YLYDM>k1dr|A5PcA8t zJoR`ay<7*W9y*;%T(u*0oPVIyqwe2cP;gZbf2H86-b(!^c~tr_LGQGeq`*l|)m{py z0tAXwdpV)Z*KH`1{gVEZ=s%|5laZ!CxY{q8QwmO8-YhzCc?0Ui$@fxl;v^ddC%)Q- zIj)XRvx2K~eku$8X$7xGT8_gV3SOb$ae>n~)8Fe!`fn-t^$LDK!PWKREd^KCska4A z{Z;+{6lK2Db!w`>Ngj3GEK}&!dgS{UqF477HCgE8xr69c`NImmy58Nd;HqE!hQO(w zsd(hLeOtj7D|jvuBG5Rfe#Rqkiav-(*7HFHC*6~Hzk;jt@`}KzeszB3ew&Kv7F^a- ztk@UHDZbH7e!zjt``R)GPCkf2g#-6-RMV;*ILSpp?%PQoA08>s?GAdHdkX6WPW+4T zNcvUz#1;T_1@O&fx zwQ;!?5VTRjX^dp}v`A1rVrV(?6$KX}bL3wXyg;NWQvI#8QNrJ$9|Wp@5>MGbw<>-~ zjv19aqTo{$`Y$SYk%D{0piq>4YlOXog45rK%dk$tX?|r0DLB!S&*h;avZXWHyg{qL zWvU%-l&NT_FRE3vHb2$?i5$>iMshSBl z%yg4Ya`-E7V?VC0jHFt{y)7-`5;X|fJ3?7&Rh!LnA@~=*zqb2@_Yvv-|1|}14)9+` zmdrEw*OFzN#$^-YFM1oEQKA1jl8hUS38wu2xuXB$7I3+D^J$-v@-0~g?N6(jtXB1G zf2PVWOZzWZa#iG>u8#KmLLFOF3~w>Q(Vk5cWw1RyJn0{Oqai1)33zN@3+z1~q)v~+iNhPtrc;LmtM75Ieiqx?OK$m8u##inpKUTj8XFG3rGZQ*9T ztE{*;+>MU>dBQL@10@-%su7r6`)S`)hDYu-WSotk34Xb@6aVzf;9n*9y&@+58JEG| znt`9r=TdGNXOlmgfuGLSl269j`1=LF+J8EiNj@28ze(`>P^JdizO3idNMz#QA^6oZ3)Lg}WSkWI&hdLf@XPPN zqzumSI}iR$e);{E6x%uerQpvb zKmIit53dSF8h^TnbIKnP?K}1FZn3ETG5e!`ut4%V`TGRFU$ieH*}j}d`u`4@`u|feK+edqiQiHVyj8y2yYTf)DU@=HL(hXsF+ z19$SjCgpd~mD0-%gp+o}$M*)Lo${Z~(*7sO*bwCRdUkH2|2v+g{lY4n!;5kXG=4Mj zINNu@PN`4AT_9zLcS+e0=-q@Goc+H_@XPyONn1mggb2b6?F@c>sg#|#Rq$Wwz&8ng zC%x1k886MkzavZgvi~xs`#I{9Q~p80FX^Os&ioqiXVTx&kJ}tlPig#b!Xx=5{&~^9 zbDW+P?F%<7o9D7pM$W_E3jP# -#include +#include "./modules/mesh/mesh.h" +#include "modules/fluids/fluids.h" -#include -//#include "./numerics/interpolation/interpolation_linear.h" +//#include +//#include +//#include + int main(int argc, char const *argv[]) { - utils::Md A; + utils::Vd a = utils::linspace(1, 10, 10, true); + a.print(); + mesh::Mesh1D mesh(a); + mesh.generate_vertices(0.5, 10.5); + double Gamma = 1.0; + -/* - int hw = omp_get_max_active_levels(); - if (hw <= 1) return 0; - - const uint64_t m=512, k=512, p=512; // ~134M MACs; adjust if needed - utils::Md A(m,k,1), B(k,p,1), C(k,p,1); - - omp_set_max_active_levels(1); - - auto t0 = std::chrono::high_resolution_clock::now(); - for (int i = 0; i < m*k*p; ++i){ - A==B - } - double t1 = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); + utils::Md A; + utils::Vd b, s(10,1); - omp_set_max_active_levels(2); - auto t0 = std::chrono::high_resolution_clock::now(); - for (int i = 0; i < m*k*p; ++i){ - A==B - } - double t1 = std::chrono::duration(std::chrono::high_resolution_clock::now() - t0).count(); - - omp_set_num_threads(prev); - - // Must not be notably slower with many threads - CHECK(tN <= t1 * 1.05, "rows_omp: multi-thread slower than single-thread"); + core::Configs& cfg = core::Configs::defaults(); + cfg.grid = core::GridKind::Uniform; + cfg.left = {core::FDKind::Forward, core::BCKind::Neumann, 0.0}; + cfg.right = {core::FDKind::Backward, core::BCKind::Neumann, 0.0}; + cfg.solver = core::SolverKind::LU; + fluids::Diffusion1D diffusion(cfg, mesh, Gamma); + diffusion.assemble(A, b, s); - utils::Md A(5,5, 1); - utils::Md B(5,5, 1); - utils::Md C(5,5, 2); - - bool result1 = (A==B); - bool result2 = (A==C); - - omp_set_max_active_levels(1): - - for (int i = 0; i < 100; ++i){ - (A==B) - } - - omp_set_max_active_levels(2): - - for (int i = 0; i < 100; ++i){ - (A==B) - } - - std::cout << result1 << std::endl; - std::cout << result2 << std::endl; - - - - -*/ - - - /* - - utils::Vector x(100, 0), y(100,0); - for (uint64_t i = 0; i < 100; ++i){ - x[i] = i+1; - y[i] = 1 + i*0.1; - } - //y[9] = 1.5; - - x.print(); - y.print(); - - double p = 5.5; - - - - numerics::interp_linear linear(x,y); - numerics::interp_polynomial polynomial(x,y, 3); - numerics::interp_cubic_spline cubic_spline(x,y); - numerics::interp_rational rational(x,y,2); - numerics::interp_barycentric barycentric(x,y, 2); - - std::cout << "=== interpolate: " << p << " ===" << std::endl; - - std::cout << linear.interp(p) << std::endl; - std::cout << linear.interp(p) << std::endl; - std::cout << polynomial.interp(p) << std::endl; // error = polynomial.dy - std::cout << cubic_spline.interp(p) << std::endl; - std::cout << rational.interp(p) << std::endl; - std::cout << barycentric.interp(p) << std::endl; - - p += 0.01; - - std::cout << "=== interpolate: " << p << " ===" << std::endl; - - std::cout << linear.interp(p) << std::endl; - std::cout << polynomial.interp(p) << std::endl; - std::cout << cubic_spline.interp(p) << std::endl; - std::cout << rational.interp(p) << std::endl; - std::cout << barycentric.interp(p) << std::endl; - - p += 0.01; - - std::cout << "=== interpolate: " << p << " ===" << std::endl; - - std::cout << linear.interp(p) << std::endl; - std::cout << polynomial.interp(p) << std::endl; - std::cout << cubic_spline.interp(p) << std::endl; - std::cout << rational.interp(p) << std::endl; - std::cout << barycentric.interp(p) << std::endl; - - p += 50.01; - - std::cout << "=== interpolate: " << p << " ===" << std::endl; - - std::cout << linear.interp(p) << std::endl; - std::cout << polynomial.interp(p) << std::endl; - std::cout << cubic_spline.interp(p) << std::endl; - std::cout << rational.interp(p) << std::endl; - std::cout << barycentric.interp(p) << std::endl; - -*/ + return 0; } \ No newline at end of file