From 70bdc34b4bf158e415921cdbd76058355138b7d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 8 Nov 2019 16:45:50 +0100 Subject: [PATCH] Remove obsolete code generation for periodic case. This is not needed anymore, because the code that is generated for the non-periodic case also works for a periodic domain. Moreover, the old special code was inefficient, not well tested, and it was not working if the spline degree was not the same along all directions. To verify that the periodic case is treated properly, a unit test is added: 2D Poisson's equation on a unit circle, with periodic BCs along the theta direction, homogeneous Dirichlet BC at r=1, and different spline degrees along r and theta. --- mesh/circle.h5 | Bin 0 -> 14400 bytes psydac/api/ast/fem.py | 112 +++--------------- .../api/tests/test_api_2d_scalar_mapping.py | 32 ++++- 3 files changed, 45 insertions(+), 99 deletions(-) create mode 100644 mesh/circle.h5 diff --git a/mesh/circle.h5 b/mesh/circle.h5 new file mode 100644 index 0000000000000000000000000000000000000000..daebcc8ac460a5cba36ceb8b2f64177dcad5a9ab GIT binary patch literal 14400 zcmeHNXIK^I);%;O1RDxs#}0N76Xo)`kBXXD5er6CL_G>oihv?k>@_NuSfa6a#oo~^ z*bt+l(!@p;g$TxqCGj3+uQ8(6B;WmUpM3XoemJww-tXFLt#{6x**wnr%B5Wivx;V# z!o$>5W2`ByF;D{e=M4^R^P;g@G5{Cqnwf~V7|Ml#9<4BEjnB*I@e7al>w37gZL85& zs{Pb6mT8(8%aw|we*WqF|HFZ=+O~2nRN5PIaXV^aaG$9z4AuGj>VgA98-@A~FZ?Xn zKfr&u|IoLwuU9~TkKa(wCRm`}w>F61KN)C>$@c+;1{zbPL|td3j|;EM>i_AfKbs&n zG&EF7OmXKM84QqxO4aXOEZ)Ie(DzI~jb_>J$-m)M_by#s3^awR<*+gOW0JN=&nzYv zs&512F`(Iv&(-hsXp=w2`)|qqZ|^;|~DYBW`G{vEHlf;U(uKKk!*SLb;XPP>CMuA?zKK_1Q%CoVsrL}%RA+th!VQ8vp*Q!%@{g1Jxhg;9W zlaf$i&dT$R}ZH;=S|m!9KoDHV$S1>e(xJkhd<#$*iG`pO>%B$)>BXZm5@; zjg9(X_8IK-_PUAw^We|`eO32>;XZ!4&hqd!8|XjMZ?IQjXd#N%XrCY_o4QTv>T7kQ zgY}mU?F+93y>ELjUtcdzd;JFig`L%Y|7zWb%N#zsOsV~7yQY6--G}WQ-(NP!-_K7s zDA;F|PjIM{&Dil~iivd{>Y5E!Y*Fe1^)?3U0`(^73*P!QVx(7~&P;#L>jHiJ1D*7V zckO=n`n?T%8%>v3Aw2c+09WCR$lnTx;FXrN^u4 zpJ}C^nRel}eG8Mvnrj88q_%`@%_|2IN@=BUY3*sZ^E1vbEdyK&?W8gWbq~9g1%5eg z%edxeH~mr$eJr)F=LP7N{Ah)K<+VpTZ!wvzb#|!1#cx|kt zRsHyM{^Yx9|6kpKBc7;D!6+bDE+2 zjZL#gZLOR?oYO`Z$^WN=_vd>}**|80F<#T#@{jA(@3*(rAJ_f2WdH2Gstuv)`QN{< z*5STVCu9NsyWLkG`q5k=YG5Rj<@&zO^MBv30k{=YtG#=#edt&3chjG)KRNKp zfq&sZ>%#^OAIGhsIloLy?YnFohaR^!vaw%F299~(RXtK6QqB8Gp5`etF=H8osQ zys55W4(&>rKe7ImAE;ZWO#6gpi>XWZxwldbmQvr%LC1Q3x10nooW`|UmR4|f6$Luw z6I`nRho&C?ey(A?eI$5Wsk*(xR;Qi^DZ*VhvPSA^x^Oh&*`n!(sa)>o z%cm@uMqMrEB~ROUgbXESSHzp8o_|8Uy~1xRt9WCiAFKULgbz;vAD!^!>fq}mdbk$! zK(Fj?*eK_r$eY6tb$vBi&niCj@&xpsE_%5I_*E9Yd;t8Rm+L@Z7tzQ0IieqW`3LwD zD|)%H;86|WjfTH7(9anDPZGU682eKfe0$VsUC`L%3>7@8INaJGzqsZH{--Hfrj!-S{{Keg(Bgrw`H7DXB zQ-XqbfF9Vsnd`AKn#!KHJ5+hWHtOEaxuAByY5HQh<2K)da5`0{`W+X~Gt^LGjz(-M z_3VcFc){a*;8l};`~~{$5I#H%d=i8&R|Q|_C zY={00q(2MZ0ft_-27l<~4$v1R`nVGGpGN#9M*MLVy&MPs9*JJA0)L^Gx50m3(aZ7J zpX%T{+VIr8t-cqj|6_lbZYEnfyyGyF0?SKyKKnS2LcWdU=~cTG+6vvdD*|5{&ED7ETa&4hIQ=w!m)eefQP;S{? za%o)ssllqBD91Tt#&X-`lx+R#%GunDlq50VKwMqwxi#wNEBsDyci?48KUVwo5k7nY ze4YzmR{4$;J=_y|7KvWj->_M#UTy(D3b4MW=;eOEj}d&%DcYZxfL|x!$MwO#jqv9; z(05(*aTN5YAg%;|)cU!Cr~1cf@HbfeW!3*n;PV*!V+X!n^VcpclY5o!Ja6A}zvXVW zO&DwB;(iVL`%PpQL(WWU(^`iKxFQKI#y1-=vA z)^ff^Z6xN`h>uA)9{8Uzci0Wbw>Bg}3bE>Bg;%=dMp8gUx|Ci5&(x_e1{Ir%ImI z=G}eOcUb~eYiYUYjMFLVb9Zm{#cJ2+(Y@uhOkyw4f=6o}&RutzN=><=X&8QmW+-^) z$aJT5r3<6ussDoZ-&bpXoWjFKKJPm=k?OWyvv&09^)zEql80_`5?M>kGZFWadY*^+ zX@bXQz?&ugxHS5K4?hQ=mBN>8z_-8XVN=mdGes}HzbjR}{8IGtDy*+4>)BKE@-KqV z^U=Sg^yk^&S4{Zv5b$3h{JA;wK`$FZKlHLI{8@cf zz5y3%*v;vFgIYfsI_knd4sljfQ523y;&5pU|!$NWVzI(=cbq&lrJ%>@ffLRQ`Dyk9=8DA4(Z3Q(QmKt zVN>vFBz(Cp_(Cs#2|f2juk5dOmGdys%QfIfQ>?Eh>vp_iwE zKlF0R#fn}oCwe&^`pY98S=1lC4gZ>pUWQ&WC{X-W`A~82vBUn<1K;s220b}9;|>Ks z3b9DKbd=ZHXB(ERn)>$rO3iyno}!L~x$g|TMfK|>-MQU7mP~Gaf6lsG3Vm1n?Cp?8 zSE=H>eT|QYCs5Y#$i`EiB$437(U9X4y3fkFK|l5KoMrj*DVpHZ=SlM|w<$iW%ft;& zwvpRK{e)3Gm7*nPHI9*bR_lLM_{Fg|@G3|@zK4Fs!iU4b=aleel`r(LH}tF%z4Cn4 ztaF|qey|(-h{pPwqL)7dKJ;=G^v{(3ydL}p3qQUJ{^LY1=Ru#V=;h7O-yQKO_%jCe zD+P~+1J4@%W}_eU(qhrejQwc@zA@LXcpBbLqm2(|M5KG2;BH@cFmiB7SKqH9&muod z^0edfrs$kqsnob)^Xa*VE|Pol%r9HcyG7;04W~JTCeTc4-JF6E$#glbUi(IMZctYR z?<#$NW{UeMK#b@c2-;mJmo@un>%Q}8X=SGlK*HX0>*Gs1#C1y3Yka||@ zy9geuc$1_b$D^Nz@Zl=pvq|`J6YzyzUI9I>q8Gov8=TLHUJipFTd=;CtY=T)TMItV zK>th9pWVQ3k?3U`@P8?K`6BdLie4TA{pAstf;br z`>Cn4rXtmU@Ft1!)>&mYb-O`FdoAyO-TgKRUOZ(Ed32^`L>djrTG76e^Lc8zb@k}k z-)7M1(-RtbTkN6KxHD_^ch00k60;iXq@LCKmkPglz6-nu(vKt2FIo8T7VrrazFZl6 zFN+=?1U-gcieA~@lN+5cA>RXj1Y^CmtmhQquN8cLjQ&TYKMP(wrwczG0sfAnm#;zJ z9ns61it1+%_@map5WSQIyj=Ku0{v9~H-nEK_NOZNCbzPvPZoD+#1TpdWg9GyL%vy}S|r{RTWY_zS%} z3H~<&pXNpPhcia+UC^(=JzDwHBWm(@r@8-vAmbko-oy7dG%^+W*^(y{+djoNFS<)M zV|PZ*tQb#?bi3NP|CT`q8XgWASmG9C+;HyNb5<&SZFB7G(thdWq2MJ`Ise+}hC{Q+ zG&s?u`rJ#jvx-GXL0C2|v+7l4!nXbNZM?&Zr;F~<9EsT%@ozy&J#R&QU%_K*;N6yf ztoCaxeE2x{ys4}3s(7>Hi3DC^k~_=g0a)6u`3=;eXn zXDR%+r0C`C!k>3QAM|o4^bZ%k{8;o-1nT?1KMUa9g1?>64|?eY_-L>{wZS)Mm)o8(;FecD$FfKd zP^HS^J5+wH<*$J+(kQx5pquZr3^G^nQpomk%YeW(*|fdNm`7z+#8Iw$$5D;v-lsVn zRc_X(gXGZbzLnLX92z4rJ0jjF^{m!^E_iGQyt2}dPorP5@Zomgb4mDeZSaL2-T^(( zi|=oQbFS#+AK*s})>oGGd=vPs1)nRSf0Xp+lHgZI`0;M=&k?=60{Wnr2SERqqLc}BWb7_QkL-WP057V=x)r-DJ`I!bv%+C>b zmU>p}4=Vg_aXa9BE&W*SH(mJfFz{I*eAy0spNbxihMq{#E5E;CyPcmP9|1oWVSP1O z&ohBPRPY&k)&6_}{9J?|yMuoX(aSTT4|=%+^uG|j+#deaMEyR&Q~hJ;r555Zhrxeu z@CnEMGyvaUt1OB-JSLCsJ_@Q|F6TTCO(rqCiedAc-6y-OFPfcqO{XdK&-ISDo=H0|r}|jiW|NPCcZa@tHDYeJ*Ev*> zx%5i8*4L>2Psf8-t$RS8WtJPQE_swre$%V8`>%P_QDQ!U_=(i>0Mri=Jnjd)iPDeN ze&EB6z^9M!WjpYNUcL!E!Y?21&wHF-h+aMmKl);QOr%5{N06qs{hc-;7g6bH%yuL@-hF$V*XYCm#;=6 z^I$&7{2Pn;_W{*Y=BpGkROZ1P%)jYWS((qWF#l#_{!PXFn~VAP8ckH@yFARlM@i`I<7%USpYm^Eg+T zcaJjk@)n6P|3Bm?sVB95mCWblqs+g#Yz2SsEAwwQtNt%k=HGNaE_$hn><|AFg6v Expecting a kernel') + if isinstance(discrete_space, (tuple, list)): + space = discrete_space[0] + else: + space = discrete_space + + if not isinstance(space, (SplineSpace, TensorFemSpace, ProductFemSpace)): + raise NotImplementedError('Only Spline, Tensor and Product spaces are available') + # ... + obj = SplBasic.__new__(cls, kernel.tag, name=name, prefix='assembly', mapping=mapping, is_rational_mapping=is_rational_mapping) - - # ... get periodicity of the space - dim = kernel.weak_form.ldim - periodic = [False for i in range(0, dim)] - if not( discrete_space is None ): - if isinstance(discrete_space, (tuple, list)): - space = discrete_space[0] - - else: - space = discrete_space - - if isinstance(space, SplineSpace): - periodic = [space.periodic] - - elif isinstance(space, TensorFemSpace): - periodic = space.periodic - - elif isinstance(space, ProductFemSpace): - periodic = space.spaces[0].periodic - - else: - raise NotImplementedError('Only Spline, Tensor and Product spaces are available') - # ... - obj._kernel = kernel obj._discrete_space = discrete_space - obj._periodic = periodic obj._comm = comm obj._discrete_boundary = kernel.discrete_boundary obj._backend = backend @@ -1124,10 +1102,6 @@ def kernel(self): def discrete_space(self): return self._discrete_space - @property - def periodic(self): - return self._periodic - @property def comm(self): return self._comm @@ -1157,7 +1131,6 @@ def _initialize(self): form = self.weak_form fields = kernel.fields fields_coeffs = kernel.fields_coeffs - fields_tmp_coeffs = kernel.fields_tmp_coeffs vector_fields = kernel.vector_fields vector_fields_coeffs = kernel.vector_fields_coeffs zero_terms = kernel.zero_terms @@ -1326,56 +1299,11 @@ def _initialize(self): # ... kernel call mats = tuple(element_matrices.values()) - if not( self.comm is None ) and any(self.periodic) : - # ... - stmts = [] - for i,il,p,span,n,per in zip( indices_i, - indices_il, - test_degrees, - indices_span, - npts, - self.periodic ): - - if not per: - stmts += [Assign(i, span-p+il)] - - else: - stmts += [Assign(i, Mod(span-p+il, n))] - -# _indices_i = [i for i,s,p in zip(indices_i, starts, test_degrees)] - _indices_i = [i-s+p for i,s,p in zip(indices_i, starts, test_degrees)] - for x,tmp_x in zip(fields_coeffs, fields_tmp_coeffs): -# stmts += [Print([_indices_i, ' ', indices_i, starts])] - stmts += [Assign(tmp_x[indices_il], x[_indices_i])] - - ranges = [Range(0, test_degrees[i]+1) for i in range(dim)] - for x,rx in list(zip(indices_il, ranges))[::-1]: - stmts = [For(x, rx, stmts)] - - body += stmts - # ... - - # ... - f_coeffs = tuple(fields_tmp_coeffs) - # ... - - # ... TODO add tmp for vector fields and mapping - gslices = [Slice(sp-s-d+p,sp+p+1-s) for sp,d,p,s in zip(indices_span[::ln], - test_degrees[::ln], - test_pads[::ln], - starts)] - vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs]) - m_coeffs = tuple([f[gslices] for f in kernel.mapping_coeffs]) - # ... - - else: - gslices = [Slice(sp-s-d+p,sp+p+1-s) for sp,d,p,s in zip(indices_span[::ln], - test_degrees[::ln], - test_pads[::ln], - starts)] - f_coeffs = tuple([f[gslices] for f in fields_coeffs]) - vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs]) - m_coeffs = tuple([f[gslices] for f in kernel.mapping_coeffs]) + gslices = [Slice(sp-s-d+p, sp+p+1-s) for sp,d,p,s in + zip(indices_span[::ln], test_degrees[::ln], test_pads[::ln], starts)] + f_coeffs = tuple([f[gslices] for f in fields_coeffs]) + vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs]) + m_coeffs = tuple([f[gslices] for f in kernel.mapping_coeffs]) if is_bilinear: if not unique_scalar_space: @@ -1543,12 +1471,6 @@ def _initialize(self): fields_shape = tuple(FunctionCall('len',[p[0,Slice(None,None)]]) for p in points) for F_value in self.kernel.vector_fields_val: prelude += [Assign(F_value, Zeros(fields_shape))] - # ... - if not( self.comm is None ) and any(self.periodic) : - for v in self.kernel.fields_tmp_coeffs: - stmt = Assign(v, Zeros(orders)) - prelude += [stmt] - # ... # ... if self.debug: diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 69ba30c5d..2706f7959 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -30,13 +30,20 @@ # ... #============================================================================== -def run_poisson_2d_dir(filename, solution, f, comm=None): +def run_poisson_2d_dir(filename, solution, f, dir_boundary=None, comm=None): # ... abstract model domain = Domain.from_file(filename) V = ScalarFunctionSpace('V', domain) + if dir_boundary is None: + B_dirichlet = domain.boundary + elif len(dir_boundary) == 1: + B_dirichlet = domain.get_boundary(**dir_boundary[0]) + else: + B_dirichlet = Union(*[domain.get_boundary(**kw) for kw in dir_boundary]) + x,y = domain.coordinates F = element_of(V, name='F') @@ -56,7 +63,7 @@ def run_poisson_2d_dir(filename, solution, f, comm=None): l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - bc = EssentialBC(u, 0, domain.boundary) + bc = EssentialBC(u, 0, B_dirichlet) equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) # ... @@ -353,7 +360,6 @@ def test_api_poisson_2d_dir_collela(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -# TODO h1 norm not working => regression due to the last changes in sympde def test_api_poisson_2d_dir_quart_circle(): filename = os.path.join(mesh_dir, 'quart_circle.h5') @@ -370,7 +376,25 @@ def test_api_poisson_2d_dir_quart_circle(): expected_h1_error = 0.009473407914765117 assert( abs(l2_error - expected_l2_error) < 1.e-7) -# assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_circle(): + filename = os.path.join(mesh_dir, 'circle.h5') + + from sympy.abc import x,y + + solution = (1 - (x**2 + y**2)) * cos(2*pi*x) * cos(2*pi*y) + f = -laplace(solution) + + dir_boundary = [{'axis': 0, 'ext': 1}] + l2_error, h1_error = run_poisson_2d_dir(filename, solution, f, dir_boundary) + + expected_l2_error = 0.0015245737751297718 + expected_h1_error = 0.06653900724243668 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== def test_api_poisson_2d_dirneu_identity_1():