Time-harmonic electromagnetics with exact controllability and discrete exterior calculus

This is a sel -a chi ed e sion o an o iginal a icle. This e sion may di e om he o iginal in pagina ion and ypog aphic de ails. Au ho (s): Ti le: Yea : Ve sion: Copy igh : Righ s: Righ s u l: Please ci e he o iginal e sion: CC BY 4.0 h ps://c ea i ecommons.o g/licenses/by/4.0/ Time-ha monic elec omagne ics wi h exac con ollabili y and disc e e ex e io calculus © The Au ho s 2023 Published e sion Mönkölä, Sanna; Räbinä, Jukka; Rossi, Tuomo Mönkölä, S., Räbinä, J., & Rossi, T. (2023). Time-ha monic elec omagne ics wi h exac con ollabili y and disc e e ex e io calculus. Comp es Rendus. Mecanique, Online i s . h ps://doi.o g/10.5802/c meca.234 2023 Comp es Rendus Mécanique Sanna Mönkölä, Jukka Räbinä and Tuomo Rossi Time-ha monic elec omagne ics wi h exac con ollabili y and disc e e ex e io calculus Published online: 13 Decembe 2023 h ps://doi.o g/10.5802/c meca.234 Pa o Special Issue: The scien i ic legacy o Roland Glowinski Gues edi o s: G egoi e Allai e (CMAP, Ecole Poly echnique, Ins i u Poly echnique de Pa is, Palaiseau, F ance), Jean-Michel Co on (Labo a oi e Jacques-Louis Lions, So bonne Uni e si é) and Vi e e Gi aul (Labo a oi e Jacques-Louis Lions, So bonne Uni e si é) This a icle is licensed unde he C ea i e Commons A ibu ion 4.0In e na ional License. h p://c ea i ecommons.o g/licenses/by/4.0/ Les Comp es Rendus. Mécanique son memb es du Cen e Me senne pou l’édi ion scien i ique ou e e www.cen e-me senne.o g e-ISSN : 1873-7234 Comp es Rendus Mécanique Published online: 13 Decembe 2023 h ps://doi.o g/10.5802/c meca.234 The scien i ic legacy o Roland Glowinski / L’hé i age scien i ique de Roland Glowinski Time-ha monic elec omagne ics wi h exac con ollabili y and disc e e ex e io calculus Élec omagné ique ha monique empo elle a ec con ôlabili é exac e e calcul ex é ieu disc e Sanna Mönkölä ,a, Jukka Räbinä aand Tuomo Rossi ,∗,a aFacul y o In o ma ion Technology, Uni e si y o Jy äskylä, P.O. Box 35, FI-40014 Uni e si y o Jy äskylä, Finland E-mail: uomo.j[email p o ec ed] (T. Rossi) Abs ac . In his pape , we apply he exac con ollabili y concep o ime-ha monic elec omagne ic sca e - ing. The p oblem is p esen ed in e ms o he di e en ial o ms, and he disc e e ex e io calculus is u ilized o spa ial disc e iza ion. Acco dingly, he physical p ope ies o he p oblem a e main ained. Despi e we con- side ime-ha monic p oblems, we concen a e on ansien wa e equa ions ea ed by he exac con ollabil- i y app oach. Essen ially, we use a con olled a ia ion o he asymp o ic app oach wi h pe iodic cons ain s, in which he ime-dependen equa ion is simula ed in ime, un il he ime-ha monic solu ion is eached. Résumé. Dans ce a icle, nous appliquons le concep de con ôlabili é exac e à la dispe sion élec omagné- ique empo elle. Le p oblème es p ésen é en e mes de o mes di é en ielles e le calcul ex é ieu disc e es u ilisé pou la disc é isa ion spa iale. En conséquence, les p op ié és physiques du p oblème son main- enues. Bien que nous considé ions des p oblèmes ha moniques empo els, nous nous concen ons su les équa ions d’ondes ansi oi es ai ées pa l’app oche de con ôlabili é exac e. Essen iellemen , nous u ili- sons une a ia ion con ôlée de l’app oche asymp o ique a ec des con ain es pé iodiques, dans laquelle l’équa ion dépendan du emps es simulée dans le emps, jusqu’à ce que la solu ion ha monique empo elle soi a ein e. Keywo ds. Maxwell equa ions, Elec omagne ic sca e ing, Di e en ial o ms, Disc e e ex e io calculus, Exac con ollabili y. Mo s-clés. Équa ions de Maxwell, Di usion élec omagné ique, Fo mes di é en ielles, Calcul ex é ieu dis- c e , Con ôlabili é exac e. Funding. Academy o Finland (G an ag eemen nos. 259925 and 260076). Published online: 13 Decembe 2023 ∗Co esponding au ho . ISSN (elec onic) : 1873-7234 h ps://comp es- endus.academie-sciences. /mecanique/ 2Sanna Mönkölä e al. 1. In oduc ion Wa e equa ions a e impo an o modeling acous ic, elas ic, elec omagne ic, and quan um mechanical sys ems in di e en ields o science, enginee ing, and echnology. The e has ecen ly been g owing in e es in di e en ial o m -based app oaches o conside ing wa e equa ions (see, e.g., [1–3]). The e a e se e al easons why he amewo k is a good al e na i e o he adi ional ec o ield p esen a ion. The me ic- ee na u e o he di e en ial o ms allows he cons uc ion o di e en ial ope a o s ha a e independen o he coo dina e sys em [4]. Fu he , he disc e e spaces and exac di e en ial ope a o s can be cons uc ed o mimic hei con inuous coun e pa s (see, e.g., [5, 6]). This p ope y also implies ha ce ain p ope ies o he sys em, such as, ene gy, a e conse ed, p o ided ha he scheme is s able. In his s udy, we conside ime-ha monic elec omagne ic wa es. We p esen he Maxwell equa ions as a i s -o de sys em in e ms o he di e en ial o ms (see, e.g., [7]) o elec omag- ne ic sca e ing p oblems. Tha is, ins ead o he ec o p esen a ion E=(E1,E2,E3) o he elec- ic ield, we p esen a 1- o m e E=E1dx1+E2dx2+E3dx3, whe e dis he ex e io de i a i e. Fu - he , he elec ic lux densi y is conside ed, ins ead o he ec o D=(D1,D2,D3), as a 2- o m e D=D1dx2∧dx3+D2dx3∧dx1+D3dx1∧dx2, whe e ∧is a gene aliza ion o he c oss p oduc , known as he wedge p oduc o ex e io p oduc . We also p esen a 1- o m e H o he magne ic ield and a 2- o m e B o he magne ic lux densi y. Fo he spa ial disc e iza ion, he disc e e coun e pa s o he a iables p esen ed as di e en- ial o ms a e appoin ed a geome ic objec s such as poin s, lines, su aces, and olumes. We apply, in pa icula , he disc e e ex e io calculus (DEC) ollowing he g oundwo k p esen ed by Ma sden and his g oup [8] and pionee ed o elec omagne ics simula ions by Bossa i and Ke - unen [9]. This cons uc ion p o ides, e.g., he conse a ion o ene gy, he elimina ion o non- physical modes, and exac di e en ial ope a o s a he disc e e le el [10]. The o he ela ed ap- p oaches include he ini e in eg a ion echnique (FIT), he disc e e geome ic app oach, and mixed ini e elemen echniques wi h Nédélec and Ra ia –Thomas elemen s (see [11] and he e e ences he ein). F om a ma hema ical poin o iew, i is con enien o use a ou -dimensional ( h ee spa ial dimensions and one empo al dimension) space– ime p esen a ion. The amewo k also allows he applica ion o he gene alized heo ies p oposed in [12,13]. S ill, he ou -dimensional sys em leads o la ge nume ical sys ems ha also equi e sophis ica ed compu a ional g ids. To demon- s a e how he con en ional h ee-dimensional p oblem is inhe i ed om he gene alized wa e model, we i s p esen he ou -dimensional backg ound a he modeling phase and go on o p oceed in ime wi h h ee-dimensional spa ial disc e iza ion. Ea lie we ha e p esen ed a nu- me ical scheme o ansien elec omagne ics wi h he DEC and high-quali y g ids [14, 15]. In his pape , we sys ema ically apply and es hese echniques o simula ing ime-ha monic elec- omagne ic sca e ing. We ollow he exac con ollabili y concep (see, e.g., [16]) and sol e he ime-ha monic p ob- lem by using he model p esen ed in he ime domain. This concep has a his o y da ing back o he ea ly 1990s, when Glowinski e al. (see [17–19]) i s conside ed a nume ical me hod based on i . I s i s applica ion in h ee-dimensional elec omagne ics was, o ou knowledge, p esen ed in he INRIA epo [20]. The d awback o he app oach was ha i used H1-con o ming ini e elemen s, ende ing i applicable only o limi ed classes o p oblems due o spu ious modes ha can a ise in he nume ical solu ion. Addi ional compu a ional cos s a e due o he second-o de o mula ion o he Maxwell’s equa ions ha was used, leading o slowe con e gence and he need o p econdi ioning o he minimiza ion p ocess, since he ene gy space is no o L2- ype. To imp o e he me hodology, Glowinski and Rossi [21] conside ed he wa e equa ion as a i s - o de sys em, and hey de i ed he associa ed exac con ollabili y scheme. The main bene i o Sanna Mönkölä e al. 3 hei app oach is ha he sys em’s ene gy space, whe e he minimiza ion is ealized, is o L2- ype, and he minimiza ion con e ges wi hou p econdi ioning. The d awback, howe e , is he need o use mixed me hods (e.g. Ra ia –Thomas ini e elemen s [22]) in disc e izing he p oblem. Recen ly, a i s -o de o m based on a discon inuous Gale kin me hod was de eloped o he h ee-dimensional Maxwell case [23]. Ins ead o using mixed ini e elemen s, ano he app oach, based on DEC, s a ed o p og ess wi h he PhD hesis o Räbinä [24]. The idea o his hesis a ose om he heo e ical a icle [25]. Räbinä also conside ed exac con ollabili y echniques o he compu a ion o ime-ha monic wa es. Se e al imp o emen s o he ime in eg a ion o he Maxwell’s sys em, such as adap i e local ime-s epping, ha monic Hodge s a app oxima ion, and spa ial ilings inspi ed by c ys al- log aphy, we e hen in oduced in [14]. In [15], he DEC app oach was u he gene alized o a class o wa e p oblems co e ing acous ics, elas odynamics, elec omagne ism, and e en quan- um mechanics ( he Weyl equa ion). Fu he applica ions in quan um mechanics we e hen con- side ed in [26–28]. The sys ema iza ion o scien i ic so wa e de elopmen co e ing undamen- al conse a ion laws was ske ched in [29], and a gene ic ca ego y- heo e ic amewo k o space– ime linea wa e phenomena was inally p oposed by [30] in 2022. Such o malisms allow, e.g., mo ing and de o ming compu a ional domains [29]. Ano he esea ch ack conside s highe - o de disc e e ex e io calculus in ol ing highe o de Whi ney o ms and no el Hodge s a o - mula ions [31–33]. The exac con ollabili y app oach was ex ensi ely s udied in he la e 2000s and ea ly 2010s by Mönkölä e al. In pa icula , highe -o de spec al elemen disc e iza ions and luid-s uc u e in e ac ion p oblems we e conside ed in [34–36]. The app oach was ecen ly ex ended o h ee- dimensional isco-elas ic equa ions [37]. Essen ially, he app oach is a con olled a ia ion o he asymp o ic app oach wi h pe iodic cons ain s, in which he ime-dependen equa ion is simula ed in ime un il he ime-ha monic solu ion is eached. The ime disc e iza ion is ealized wi h a s agge ed leap og- ype scheme wi h asynch onous ime s eps (see, [14,38,39]). In p ac ice, he esidual o he exac con ollabili y algo i hm de ines a each i e a ion how a he solu ion is om a pe iodic one. I also gi es an an impulse o he sys em o accele a e he con e gence a e. Roughly speaking, he numbe o compu a ional ope a ions g ows linea ly wi h he numbe o deg ees o eedom in ol ed in he spa ial disc e iza ion. The es o his a icle is o ganized as ollows. In Sec ion 2, we p esen he ou -dimensional di e en ial o m o mula ion. The model is simpli ied in o sepa a e ime and space dimensions in Sec ion 3, whe e we also ca y ou disc e iza ion. To elimina e disc e iza ion e o s, we ad- jus he scheme o ime-ha monic p oblems wi h spa ial and empo al co ec ions. This ype o s a egy has been used o imp o e he accu acy o ime disc e iza ion in elec omagne ics, e.g., by Ma and Chen [40] and Appelö e al. [41]. The spa ial co ec ion e m is de i ed wi hin he disc e e ex e io calculus ha we use o space disc e iza ion, and he empo al co ec ion e m a ises om adap ing he leap og ime disc e iza ion o ime-ha monic equa ions. Ne e heless, bo h co ec ion e ms can be encapsula ed in o disc e e Hodge ope a o s. In Sec ion 4, we ocus on a ime-ha monic p oblem and conside how i can be sol ed e icien ly in he ime domain wi h he exac con ollabili y me hod. The p oblem is e o mula ed as a leas -squa es op imiza ion p ob- lem ha is sol ed by he conjuga e g adien algo i hm. In Sec ion 5, we apply he app oach o h ee-dimensional elec omagne ic sca e ing p oblems and es he pe o mance o he me hod wi h se e al ypes o compu a ional meshes. The concluding ema ks a e p esen ed in Sec ion 6. 2. Model We conside a linea combina ion o he di e en ial o ms in he ou -dimensional Euclidean space (x0,x1,x2,x3)Tas ˜ F=˜ u1+˜ u2dx0+˜ u3dx1+˜ u4dx2+˜ u5dx3+˜ u6(dx0∧dx1)+˜ u7(dx0∧dx2)+ 4Sanna Mönkölä e al. ˜ u8(dx0∧dx3)+˜ u11(dx2∧dx3)+˜ u10(dx3∧dx1)+˜ u9(dx1∧dx2)+˜ u14(dx0∧dx2∧dx3)+˜ u13(dx0∧ dx3∧dx1)+˜ u12(dx0∧dx1∧dx2)+˜ u15(dx1∧dx2∧dx3)+˜ u16(dx0∧dx1∧dx2∧dx3), whe e he coe icien s ˜ ui,i=1,...,16 a e scala - alued. The basis 1- o ms a e he coo dina e di e en ials dx0,dx1,dx2, and dx3, and highe o de basis k- o ms, k=2,3,4, a e cons uc ed as ex e io p oduc s o kbasis 1- o ms, such ha , dxi∧dxi=0 and dxi∧dxj=−dxj∧dxi. This o mula ion allows us o p esen a se o linea models wi h he i s -o de di e en ial ope a o in a sho o m, (d+δ)˜ F=˜ J, (1) whe e d=P3 i=0(dxi∂i∧) is he ex e io de i a i e, δis he code i a i e o d, and ˜ Jis a sou ce e m p esen ed as a linea combina ion o he di e en ial o ms. Acco dingly, Equa ion (1) is equi alen o he sys em                −∂0∇· −∂0∇· −∂0−∇× ∇ −∂0−∇× ∇ ∇ ∇× −∂0 ∇ ∇× −∂0 ∇· −∂0 ∇· −∂0                               ˜ u15 ˜ u16 ˜ u3 ˜ u6 ˜ u11 ˜ u14 ˜ u1 ˜ u2                =                −˜ b16 ˜ b15 −˜ b6 ˜ b3 −˜ b14 ˜ b11 −˜ b2 ˜ b1                , (2) whe e ˜ u3:=  ˜ u3 ˜ u4 ˜ u5 ,˜ u6:=  ˜ u6 ˜ u7 ˜ u8 ,˜ u11 :=  ˜ u11 ˜ u10 ˜ u9 , and ˜ u14 :=  ˜ u14 ˜ u13 ˜ u12 , and ˜ b3:=  ˜ b3 ˜ b4 ˜ b5 ,˜ b6:=  ˜ b6 ˜ b7 ˜ b8 ,˜ b11 :=  ˜ b11 ˜ b10 ˜ b9 , and ˜ b14 :=  ˜ b14 ˜ b13 ˜ b12 . We in e p e ∂0as a ime de i a i e and ∇,∇·, and ∇× as h ee-dimensional g adien , di- e gence, and cu l ope a o s, espec i ely. By he choice ˜ u1=˜ u2=˜ u15 =˜ u16 =0, ows 1, 2, 5, and 6 o Equa ion (2) gi e he comple e se o he Maxwell equa ions wi h he magne ic ield s eng h ˜ u3=(H1,H2,H3)T, he elec ic ield s eng h ˜ u6=(E1,E2,E3)T, he magne ic lux densi y ˜ u11 =−(B1,B2,B3)T, and he elec ic lux densi y ˜ u14 =(D1,D2,D3)T, espec i ely. The magne ic cu en densi y is −˜ b14, he elec ic cu en densi y is −˜ b11, and magne ic and elec ic cha ges a e modeled by −˜ u15 and −˜ u16. A Maxwell-like sys em is also ob ained om ows 3, 4, 7, and 8 o Equa ion (2). We e u n o Equa ion (1) and conside ˜ F=e H+dx0∧e E−e B+dx0∧e D, (3) ˜ J= −e JE−dx0∧e JH−e JρH−dx0∧e JρE−e JB−dx0∧e JD−e JρD−dx0∧e JρB, (4) whe e espec i ely, elec ic and magne ic ield s eng hs a e ep esen ed by spa ial 1- o ms e E= E1dx1+E2dx2+E3dx3and e H=H1dx1+H2dx2+H3dx3, elec ic and magne ic lux densi ies and elec ic and magne ic cu en densi ies by spa ial 2- o ms e D=D1dx2∧dx3+D2dx3∧dx1+ D3dx1∧dx2,e B=B1dx2∧dx3+B2dx3∧dx1+B3dx1∧dx2,e JE=JE1(dx2∧dx3)+JE2(dx3∧ dx1)+JE3(dx1∧dx2), and e JH=JH1(dx2∧dx3)+JH2(dx3∧dx1)+JH3(dx1∧dx2), and elec- ic and magne ic cha ges, ρEand ρH, by spa ial 3- o ms e JρE=ρE(dx1∧dx2∧dx3) and Sanna Mönkölä e al. 5 e JρH=ρH(dx1∧dx2∧dx3). Fu he , in Equa ion (3), he e a e also sou ce o ms o cu en den- si ies and cha ges co esponding o he dual componen s o he elec ic and magne ic ones, e- e ed o wi h subsc ip s Band D, espec i ely. F om he choices made abo e, we ge om Equa- ion (1) o he di e en ial o m o mula ions        de H−∂0e D=−e JE, dx0∧(de E−de B)=−dx0∧e JH, −dx0∧de B=−dx0∧e JρH, −dx0∧de D=−dx0∧e JρE, (5)        −⋆(dx0∧de H)=−⋆(dx0∧e JρB), −⋆(dx0∧de E)=−⋆(dx0∧e JρD), ⋆(dx0∧(de B−de E)) =−⋆(dx0∧e JD), ⋆(de D−∂0e H)=−⋆e JB, (6) whe e ⋆is he Hodge s a ope a o . The sys ems (5) and (6) a e dual o each o he , and in wha ollows, we p oceed wi h Equa ion (5) in sepa a e space and ime dimensions. Fo u he in o ma ion on di e en ial o ms in elec omagne ics, he eade is e e ed o, e.g., [42–44] and he e e ences he ein. 3. Disc e iza ion The quali y o he disc e e model, based on he s uc u e o he compu a ional g id, has a signi - ican impac on he e iciency o he me hod. As p esen ed in [14], we gene a e pa ly-s uc u ed non-uni o m polygonal g ids imi a ing he close packing in c ys al la ices and disc e ize he op- e a o s in he spa ial domain wi h he DEC. We disc e ize he h ee-dimensional spa ial domain by a pai o p imal and dual g ids. The elemen s o a g id a e called k-cells, whe e k=0,1,2,3 is he dimension o he cell. Basically, a k-cell is an o ien ed, con ex, polygon s yle objec , and i is de ined in a ecu si e manne by a lis o k−1-cells. A 0-cell (node) is a e ex. A 1-cell (edge) is de ined as a line segmen be ween wo 0-cells and o ien ed om he i s o he second node o he lis . A 2-cell ( ace) is a con ex polygon su ounded by a ini e se o edges. The bounda y edges a e p esen ed wi h coun e - clockwise o ien a ion, also p o iding, he o ien a ion o he ace elemen . A 3-cell (body) is a con ex polyhed on su ounded by a ini e se o aces. The g id cons uc ion has a c ucial ole in he accu acy o he me hod. In p ac ice, we gene a e a pa ly s uc u ed g id ha consis s o s uc u ed a eas sepa a ed by uns uc u ed ones. The app oach allows us o model bounda ies as accu a ely as wi h ully uns uc u ed g ids, and i is also possible o modi y elemen sizes inside he domain. The g id gene a ion is based on he Delaunay iangula ion (see [45]). Wi h he help o a Vo onoi diag am, i is possible o c ea e sophis ica ed polyhed al g ids o ill he s uc u ed pa s o he g id. By main aining he numbe o uns uc u ed elemen s a a ela i ely low le el, he g id gene a ion p ocess is conside ably sped up. By cons uc ing he dual elemen s ha a e o hogonal o he co esponding p imal elemen s, we ge diagonal disc e e Hodge ope a o s, p o iding a signi ican sa ing in compu ing ime. 3.1. Disc e e di e en ial o ms The disc e e di e en ial k- o ms (k-cochains) a e a iables associa ed by he de Rham map wi h he k-cells in he g id (see, e.g., [24, 46]) ha , in he h ee-dimensional cons uc ion, a e 1-cells (edges) and 2-cells ( aces). We deno e he i h edge o he g id by ei,i=1,...,ne, and he j h ace by j,j=1,...,m , and he co esponding dual elemen s by ∗ iand e∗ j, espec i ely. Acco dingly, 6Sanna Mönkölä e al. he componen s o he ec o - alued disc e e di e en ial 1- o ms Eand H, p esen ing he disc e ized elec ic and magne ic ields and associa ed wi h he p imal and dual edges, a e Ei=Zeie E,Hj=Ze∗ je H,i=1,...,ne,j=1,...,m . (7) Respec i ely, Dand Bp esen he disc e ized elec ic and magne ic luxes, such ha , Di=Z ∗ ie D,Bj=Z je B,i=1,...,ne,j=1,...,m , (8) a e he disc e e di e en ial 2- o ms associa ed wi h he dual and p imal aces o he g id ele- men s. Also, he ec o s o sou ce o ces, JEand JH, a e conside ed as disc e e di e en ial 2- o ms JEi=Z ∗ ie JE,JHj=Z je JH,i=1,...,ne,j=1,...,m , (9) associa ed wi h he dual and p imal aces. The ela ion be ween he disc e e 1- o ms Eand Hand disc e e 2- o ms Dand Bis p esen ed by he cons i u i e equa ions D=⋆ϵE,B=⋆µH, (10) whe e he disc e e Hodge s a ope a o s ⋆ϵand ⋆µco e he ma e ial p ope ies, pe mi i i y ϵ and pe meabili y µ, espec i ely, and he me ic p ope ies o he space o coo dina e sys em. In his pape , he disc e e Hodge s a ope a o s a e diagonal ma ices based on he o hogonali y o he p imal and dual elemen s and de ined as ⋆α=| | |e|| ⋄e|Z ⋄eαne·ned , (11) whe e neis he uni o ien a ion ec o o edge e, and ⋄eis a con ex hull including bo h e and . The i h diagonal componen o ⋆ϵis composed by se ing α=ϵ,e=ei, and = ∗ iin Equa ion (11). The j h diagonal componen o ⋆µis composed by se ing α=µ,e=e∗ j, and = jin Equa ion (11). The key componen o he disc e e ex e io calculus is he disc e e ex e io de i a i e, o incidence ma ix, d ep esen ing he neighbo ing ela ions and ela i e o ien a ions o he p imal edges and aces. We use he no a ion d1 o emphasize ha in his case, he disc e e ex e io de i a i e ope a es disc e e 1- o ms. The ma ix en y, d1ji , is non-ze o i and only i he edge eiis included in he bounda y o he ace j. Fu he , he non-ze o en ies ha e he alue ±1, indica ing he ela i e o ien a ion de ined by coun e -clockwise ci cula ion. Wi h he no a ions p esen ed abo e, he elec omagne ic sys em a e space disc e iza ion in a h ee-dimensional domain Ωis ⋆ϵ ∂E ∂ −dT 1H=−JE, in Ω×[0,τ], (12) ⋆µ ∂H ∂ +d1E=−JH, in Ω×[0,τ], (13) wi h he ini ial condi ions, E(0) and H(0), se a he ini ial ime 0. The inal s a e is conside ed a he inal ime τ. By deno ing M=µ⋆ϵ0 0⋆µ¶,K=µ0−dT 1 d10¶, (14) we can w i e Equa ions (12)–(13) wi h he ini ial condi ions as (M∂0+K)µE H¶=−µJE JH¶, in Ω×[0,τ], (15) µE(0) H(0)¶=µE0 H0¶, in Ω, (16) Sanna Mönkölä e al. 7 he solu ion o which, a ime , is µE( ) H( )¶=e− M−1KµE0 H0¶−Z 0e−( −s)M−1KM−1µJE(s) JH(s)¶ds, (17) whe e e− M−1K=I+³0⋆−1 ϵdT 1 −⋆−1 µd1 0´+³0⋆−1 ϵdT 1 −⋆−1 µd1 0´2 2! +³0⋆−1 ϵdT 1 −⋆−1 µd1 0´3 3! +³0⋆−1 ϵdT 1 −⋆−1 µd1 0´4 4! +... =Ãcosµ³⋆−1 ϵdT 1⋆−1 µd1´1 2 ¶⋆−1 ϵdT 1³⋆−1 µd1⋆−1 ϵdT 1´−1 2sinµ³⋆−1 µd1⋆−1 ϵdT 1´1 2 ¶ −⋆−1 µd1³⋆−1 ϵdT 1⋆−1 µd1´−1 2sinµ³⋆−1 ϵdT 1⋆−1 µd1´1 2 ¶cosµ³⋆−1 µd1⋆−1 ϵdT 1´1 2 ¶!. (18) The de i a i e o he o al ene gy d d E(E,H)=ET⋆ϵ∂0E+HT⋆µ∂0H=−(ETJE+HTJH) (19) is ob ained by mul iplying Equa ion (15) om he le by (ETHT). The igh -hand side o Equa- ion (19) includes he sou ce unc ions se on he abso bing bounda y o laye , such as he Sil e – Mülle bounda y condi ion [47] o pe ec ly ma ched laye (PML) [48]. By assuming hem o be o he o ms JE=⋆σEEand JH=⋆σHH, whe e σEand σHp esen he elec ic and magne ic con- duc i i ies, espec i ely, we ge (d/d )E(E,H)= −(ET⋆σEE+HT⋆σHH)≤0, implying ha he ene gy is conse ed o dissipa ed as he wa es a e abso bed om he domain. The o al ene gy o he elec omagne ic sys em is E(E,H)=1 2(ET⋆ϵE+HT⋆µH). (20) 3.2. Spa ial ha monic co ec ion To imp o e he accu acy o he simula ions in ol ing ime-ha monic wa es wi h he ime- dependence o he o m eiω , whe e ωis he angula equency and i is he imagina y uni , we o mula e spa ially ha monic co ec ions o he disc e e Hodge ope a o s. Since he elec omag- ne ic wa e speed is c=1/pµε, we can p esen ime as dis ance mul iplied by pµε. F om his basis, he concep o spa ial co ec ion is based on ob aining he componen s o Eand Hin Equa- ion (7) by in eg a ing o e ˆ Eeiω and ˆ Heiω , wi h complex- alued ˆ Eand ˆ H, ins ead o in eg a ing o e e Eand e H. Respec i ely, in Equa ions (8) we in eg a e o e ϵiˆ Eeiω and µjˆ Heiω , ins ead o e D and e B, o ob ain he componen s o Dand B. Acco dingly, by minimizing he squa ed no m o e o s based on he cons i u i e equa ions, we ind spa ially ha monic disc e e Hodge s a ope - a o s ⋆spa ϵand ⋆spa µ o sa is y he cons i u i e equa ions. These a e diagonal ma ices wi h he diagonal componen s ⋆spa ϵi=⋆ϵiκi,i=1,...,ne,⋆spa µj=⋆µjκ∗ j,j=1,...,m , (21) whe e he spa ial co ec ion e ms, κiand κ∗ j, a e κi=  1−κ ∗ 5+κ2 ∗ 56 1−κ ∗ 10 −κe 120 +κ2 ∗ 280 +κ ∗κe 1680 +κ2 e 22,400   ,κ∗ j=  1−κ 5+κ2 56 1−κ 10 −κe∗ 120 +κ2 280 +κ κe∗ 1680 +κ2 e∗ 22,400   , wi h κe=ω2ϵµ|ei|2,κ ∗=ω2ϵµ 2 ∗ i,κe∗=ω2ϵµ|e∗ j|2, and κ =ω2ϵµ 2 jand he squa ed adii o he aces ∗ iand j, 2 ∗ iand 2 j, espec i ely (see [14,24]). Respec i ely, ⋆spa σEi=⋆σEiκi,i=1,...,ne,⋆spa σHj=⋆σHjκ∗ j,j=1,...,m . (22) The la ge he g id elemen s, he mo e ad an ages a e gained by eplacing ⋆ϵby ⋆spa ϵ,⋆µby ⋆spa µ,⋆σEby ⋆spa σE, and ⋆σHby ⋆spa σHin ime-ha monic simula ions. 14 Sanna Mönkölä e al. Figu e 2. Rela i e e o o he Muelle ma ix, in eg a ed o e all sca e ing di ec ions. bounda y o he PML a he adius 2.7. To epo he accu acy o he me hod, he nea - ield so- lu ion is ans e ed o a - ield sca e ing esul s by applying a nea - ield o a - ield ans o ma- ion [51,52]. The a - ield sca e ing da a a e applied o p oduce he Muelle ma ix [53] cha ac- e izing how a ma e ial in e ac s wi h elec omagne ic wa es and p o iding in o ma ion abou sca e ing in ensi ies and pola iza ion in all sca e ing di ec ions. The accu acy is e alua ed by compa ing he Muelle ma ix solu ion wi h he analy ical Mie sca e ing solu ion [54]. Figu e 2 illus a es he ela i e e o o he Muelle ma ices compa ed o he exac Mie solu ion o each g id ype wi h i e di e en mesh esolu ions (numbe s o elemen s pe wa eleng h). We see ha he ha monic co ec ion imp o es he accu acy o he me hod wi h all he g id ypes. The mos accu a e esul s, ob ained wi h he C15- and Z-g ids and ha monic co ec ion wi h espec o space and ime, a e almos an o de o magni ude mo e accu a e han he leas accu a e esul s ob ained wi h he cubical g id and con en ional Yee’s disc e e Hodge ope a o . 5.2. Non-con ex obs acle In he second es , we use he S an o d bunny [55] as a non-con ex obs acle (see Figu e 3). Inside he bunny, we se ϵ=4 and µ=1, while he media su ounding he bunny is modeled wi h ϵ=1 and µ=1. The bunny, wi h heigh 6.173 ( om ea ip o he g ound), leng h 5.265, and dep h 4.782, is cen e ed in a ec angula p ism o heigh 9.453, leng h 9.508, and dep h 8.107. The ou e mos laye , o hickness 1.5, o he ec angula p ism ea u es he PML condi ion o Equa ion (57). In he in e io o he obs acle, we cons uc ed o each g id ype a se o g ids wi h a ying esolu ion. Ou side he obs acle, we used a egula g id wi h an edge leng h o one en h o he wa eleng h. As an inciden wa e, we used a ci cula ly pola ized plane wa e p opaga ing in he di ec ion o he posi i e x1-axis. The simula ion esul is illus a ed in Figu e 3. The ed, g een, and blue componen s in he illus a ion p esen he x1,x2, and x3componen s o he elec ic ield, espec i ely. The e can be seen a li le o e i e sca e ing wa es in he ho izon al di ec ion inside he bunny, which is in good ag eemen wi h he wa eleng h and he bunny’s wid h. To assess he accu acy, he e e ence solu ion is compu ed wi h he BCC g id wi h 40 elemen s pe wa eleng h, and he esul s o he elec ic and magne ic wa e componen s compa ed wi h he e e ence solu ion a e shown in Figu e 4. Again, we obse e a signi ican di e ence be ween he e o s compu ed wi h he con en ional Yee’s disc e e Hodge ope a o and wi h he ha monic Hodge ope a o . When he ha monic Hodge ope a o is used, he cubical g id is he mos inaccu a e, bu he e a e only sligh di e ences in accu acy be ween he o he g id ypes. In he las examples, we used he BCC ype bunny mesh wi h 50,190 nodes, 269,621 edges, and 403,215 aces o compa e he con ol me hod wi h he asymp o ic app oach and o conside Sanna Mönkölä e al. 15 Figu e 3. The S an o d bunny objec is illus a ed on he le -hand side. The x1–x3plane c oss-sec ion o he in e io sca e ed ield is illus a ed on he igh -hand side. Figu e 4. Rela i e e o o he in e io sca e ed ields. he numbe o i e a ions equi ed by he con ol me hod wi h inc easing angula equency. In hese es s, we unca ed he bunny by applying he Sil e –Mülle bounda y condi ion wi hou addi ional abso bing laye s. We ound ha despi e each i e a ion o he con ol me hod in ol es sol ing bo h he s a e and adjoin s a e equa ions, i emains mo e e icien han he asymp o ic app oach (see Figu e 5). This is because he numbe o i e a ions equi ed by he con ol me hod is less han hal o he numbe o ime pe iods needed o he asymp o ic app oach o achie e he accu acy le el o he s opping c i e ion o he con ol me hod. Figu e 6 shows he con e gence o he con ol me hod a ou di e en angula equencies. The numbe o i e a ions inc eases wi h he angula equency. 6. Conclusion We p esen ed a ou -dimensional o igin o he elec omagne ic wa e p oblem and de i ed om i a di e en ial o m o mula ion whe e he spa ial and empo al a iables a e sepa a ed. We comple ed he disc e iza ion in space by he disc e e ex e io calculus. Since ou main ocus was on ime-ha monic p oblems, o imp o e he accu acy, we applied an exac ime-s epping scheme and a ha monic co ec ion on he Hodge ope a o ma ices used a he disc e e le el. In p inciple, we used a ansien wa e sol e wi h space and ime disc e iza ion adap ed o ime-pe iodic 16 Sanna Mönkölä e al. Figu e 5. Compa ison o CPU usage be ween he exac con ollabili y me hod and he asymp o ic app oach a a ious angula equencies. Figu e 6. Numbe o i e a ions o he exac con ollabili y me hod a a ious angula equencies. p oblems. We u he accele a ed he con e gence o he ime-ha monic solu ion by using he exac con ollabili y app oach ealized by he conjuga e g adien me hod. The ene gy no m is a weigh ed L2-no m, and we minimize he disc e e quad a ic unc ional spanned by a diagonal mass ma ix. Thus, he conjuga e g adien algo i hm ope a es in an L2- ype Hilbe space, and no p econdi ioning is needed. The nume ical esul s demons a e he accu acy imp o emen s o he ha monic co ec ions. They also show he capabili y o he exac con ollabili y me hod wi h he chosen disc e iza ion s a egies in non-con ex domains. Decla a ion o in e es s The au ho s do no wo k o , ad ise, own sha es in, o ecei e unds om any o ganiza ion ha could bene i om his a icle, and ha e decla ed no a ilia ions o he han hei esea ch o ganiza ions. Sanna Mönkölä e al. 17 Dedica ion The manusc ip was w i en h ough con ibu ions o all au ho s. All au ho s ha e gi en app o al o he inal e sion o he manusc ip . Acknowledgmen s This p ojec has been pa ially unded by he Academy o Finland, g an s 259925 and 260076. The au ho s app ecia e wo anonymous e iewe s o hei insigh ul and cons uc i e ema ks. Re e ences [1] S. C. Chen, W. C. Chew, “Elec omagne ic heo y wi h disc e e ex e io calculus”, P og . Elec omagn. Res. 159 (2017), p. 59-78. [2] L. da Sil a, C. Ba is a, I. González, A. Macêdo, W. de Oli ei a, S. Melo, “A disc e e ex e io calculus app oach o quan um anspo and quan um chaos on su ace”, J. Compu . Theo . Nanosci. 16 (2019), no. 9, p. 3670-3682. [3] P. D. Boom, O. Kosmas, L. Ma ge s, A. P. Ji ko , “A geome ic o mula ion o linea elas ici y based on disc e e ex e io calculus”, In . J. Solids S uc . 236 (2022), a icle no. 111345. [4] J. B. Pe o , C. J. Zusi, “Di e en ial o ms o scien is s and enginee s”, J. Compu . Phys. 257, Pa B (2014), p. 1373- 1393. [5] D. N. A nold, P. B. Boche , R. B. Lehoucq, R. A. Nicolaides, M. Shashko (eds.), Compa ible Spa ial Disc e iza ions, The IMA Volumes in Ma hema ics and i s Applica ions, ol. 142, Sp inge , New Yo k, USA, 2006. [6] S. H. Ch is iansen, H. Z. Mun he-Kaas, B. Ow en, “Topics in s uc u e-p ese ing disc e iza ion”, Ac a Nume . 20 (2011), p. 1-119. [7] H. Ca an, Di e en ial Fo ms, Ke shaw Publishing Company, London, 1971. [8] M. Desb un, A. N. Hi ani, M. Leok, J. E. Ma sden, “Disc e e ex e io calculus”, 2005, p ep in , h ps://a xi .o g/abs/ ma h/0508341 2. [9] A. Bossa i , L. Ke unen, “Yee-like schemes on a e ahed al mesh, wi h diagonal lumping”, In . J. Nume . Model. 12 (1999), no. 1–2, p. 129-142. [10] M. Desb un, E. Kanso, Y. Tong, “Disc e e di e en ial o ms o compu a ional modeling”, Disc e e Di e . Geom. Obe wol ach Semin. 38 (2008), p. 287-324. [11] S. H. Ch is iansen, F. Rape i, “On high o de ini e elemen spaces o di e en ial o ms”, Ma h. Compu . 85 (2016), no. 298, p. 517-548. [12] R. Pica d, S. T os o , M. Wau ick, “Well-posedness ia mono onici y—an o e iew”, in Ope a o Semig oups Mee Complex Analysis, Ha monic Analysis and Ma hema ical Physics, Sp inge , Cham, Swi ze land, 2015, p. 397-452. [13] R. Pica d, S. T os o , M. Wau ick, “On a connec ion be ween he Maxwell sys em, he ex ended Maxwell sys em, he Di ac ope a o and g a i o-elec omagne ism”, Ma h. Me hods Appl. Sci. 40 (2017), no. 2, p. 415-434. [14] J. Räbinä, S. Mönkölä, T. Rossi, “E icien ime in eg a ion o Maxwell’s equa ions by gene alized ini e-di e ences”, SIAM J. Sci. Compu . 37 (2015), no. 6, p. B834-B854. [15] J. Räbinä, L. Ke unen, S. Mönkölä, T. Rossi, “Gene alized wa e p opaga ion p oblems and disc e e ex e io calculus”, ESAIM: Ma h. Model. Nume . Anal. 52 (2018), no. 3, p. 1195-1218. [16] M. O. B is eau, R. Glowinski, J. Pé iaux, “Con ollabili y me hods o he compu a ion o ime-pe iodic solu ions; applica ion o sca e ing”, J. Compu . Phys. 147 (1998), no. 2, p. 265-292. [17] R. Glowinski, “Ensu ing well-posedness by analogy; S okes p oblem and bounda y con ol o he wa e equa ion”, J. Compu . Phys. 103 (1992), no. 2, p. 189-221. [18] M. O. B is eau, R. Glowinski, J. Pé iaux, “Nume ical simula ion o high equency sca e ing wa es using exac con- ollabili y me hods”, in Nonlinea Hype bolic P oblems: Theo e ical, Applied, and Compu a ional Aspec s: P oceed- ings o he Fou h In e na ional Con e ence on Hype bolic P oblems, Tao mina, I aly, Ap il 3–8, 1992, Vieweg+Teubne Ve lag, Wiesbaden, Ge many, 1993, p. 86-108. [19] M. O. B is eau, R. Glowinski, J. Pé iaux, “Sca e ing wa es using exac con ollabili y me hods”, in 31s Ae ospace Sciences Mee ing, Ame ican Ins i u e o Ae onau ics and As onau ics, Washing on, USA, 1993. [20] M. O. B is eau, R. Glowinski, J. Pé iaux, T. Rossi, “3D ha monic Maxwell solu ions on ec o and pa allel compu e s using con ollabili y and ini e elemen me hods”, Tech. Repo RR-3607, INRIA, 1999. [21] R. Glowinski, T. Rossi, “A mixed o mula ion and exac con ollabili y app oach o he compu a ion o he pe i- odic solu ions o he scala wa e equa ion. (I) Con ollabili y p oblem o mula ion and ela ed i e a i e solu ion”, C. R. Acad. Sci. Pa is 343 (2006), no. 7, p. 493-498. 18 Sanna Mönkölä e al. [22] S. Kähkönen, R. Glowinski, T. Rossi, R. Mäkinen, “Solu ion o ime-pe iodic wa e equa ion using mixed ini e- elemen s and con ollabili y echniques”, J. Compu . Acous . 19 (2011), no. 4, p. 335-352. [23] T. Chaumon -F ele , M. J. G o e, S. Lan e i, J. H. Tang, “A con ollabili y me hod o Maxwell’s equa ions”, SIAM J. Sci. Compu . 44 (2022), no. 6, p. A3700-A3727. [24] J. Räbinä, “On a nume ical solu ion o he Maxwell equa ions by disc e e ex e io calculus”, Phd hesis, Uni e si y o Jy äskylä, 2014, h p://u n. i/URN:ISBN:978-951-39-5951-7. [25] D. Pauly, T. Rossi, “Theo e ical conside a ions on he compu a ion o gene alized ime-pe iodic wa es”, Ad . Ma h. Sci. Appl. 21 (2011), no. 1, p. 105-131. [26] J. Räbinä, P. Kuopanpo i, M. Ki ioja, M. Mö önen, T. Rossi, “Th ee-dimensional spli ing dynamics o gian o ices in Bose–Eins ein condensa es”, Phys. Re . A 98 (2018), a icle no. 023624. [27] M. Ki ioja, S. Mönkölä, T. Rossi, “GPU-accele a ed ime in eg a ion o G oss–Pi ae skii equa ion wi h disc e e ex e io calculus”, Compu . Phys. Commun. 278 (2022), a icle no. 108427. [28] M. Ki ioja, R. Zamo a-Zamo a, A. Blino a, S. Mönkölä, T. Rossi, M. Mö önen, “E olu ion and decay o an Alice ing in a spino Bose–Eins ein condensa e”, Phys. Re . Res. 5(2023), no. 2, a icle no. 023104. [29] T. Rossi, J. Räbinä, S. Mönkölä, S. Kiiskinen, J. Lohi, L. Ke unen, “Sys ema isa ion o sys ems sol ing physics bounda y alue p oblems”, in Nume ical Ma hema ics and Ad anced Applica ions ENUMATH 2019: Eu opean Con e ence, Egmond aan Zee, The Ne he lands, Sep embe 30–Oc obe 4, Sp inge , Cham, Swi ze land, 2020, p. 35- 51. [30] L. Ke unen, T. Rossi, “Sys ema ic de i a ion o pa ial di e en ial equa ions o second o de bounda y alue p oblems”, In . J. Nume . Model.: Elec onic Ne wo ks, De ices and Fields 36 (2023), no. 3, a icle no. e3078. [31] J. Lohi, L. Ke unen, “Whi ney o ms and hei ex ensions”, J. Compu . Appl. Ma h. 393 (2021), a icle no. 113520. [32] L. Ke unen, J. Lohi, J. Räbinä, S. Mönkölä, T. Rossi, “Gene alized ini e di e ence schemes wi h highe o de Whi ney o ms”, ESAIM: Ma h. Model. Nume . Anal. 55 (2021), no. 4, p. 1439-1460. [33] J. Lohi, “Sys ema ic implemen a ion o highe o de Whi ney o ms in me hods based on disc e e ex e io calculus”, Nume . Algo i hms 91 (2022), no. 3, p. 1261-1285. [34] E. Heikkola, S. Mönkölä, A. Pennanen, T. Rossi, “Con ollabili y me hod o he Helmhol z equa ion wi h highe -o de disc e iza ions”, J. Compu . Phys. 225 (2007), no. 2, p. 1553-1576. [35] S. Mönkölä, E. Heikkola, A. Pennanen, T. Rossi, “Time-ha monic elas ici y wi h con ollabili y and highe o de disc e iza ion me hods”, J. Compu . Phys. 227 (2008), no. 11, p. 5513-5534. [36] S. Mönkölä, “An op imiza ion-based app oach o sol ing a ime-ha monic mul iphysical wa e p oblem wi h highe - o de schemes”, J. Compu . Phys. 242 (2013), p. 439-459. [37] J. H. Tang, R. B ossie , L. Mé i ie , “Fully scalable sol e o equency-domain isco-elas ic wa e equa ions in 3D he e ogeneous media: a con ollabili y app oach”, J. Compu . Phys. 468 (2022), a icle no. 111514. [38] A. Lew, J. E. Ma sden, M. O iz, M. Wes , “Asynch onous a ia ional in eg a o s”, A ch. Ra ional Mech. Anal. 167 (2003), no. 2, p. 85-146. [39] A. S e n, Y. Tong, M. Desb un, J. E. Ma sden, “Geome ic compu a ional elec odynamics wi h a ia ional in eg a o s and disc e e di e en ial o ms”, in Geome y, Mechanics, and Dynamics: The Legacy o Je y Ma sden, Sp inge , New Yo k, USA, 2015, p. 437-475. [40] C. Ma, Z. Chen, “S abili y and nume ical dispe sion analysis o CE-FDTD me hod”, IEEE T ans. An ennas P opag. 53 (2005), no. 1, p. 332-338. [41] Z. Peng, D. Appelö, “EM-Wa eHol z: a lexible equency-domain me hod buil om ime-domain sol e s”, IEEE T ans. An ennas P opag. 70 (2022), no. 7, p. 5659-5671. [42] I. Lindell, Di e en ial Fo ms in Elec omagne ics, IEEE P ess Se ies on Elec omagne ic Wa e Theo y, Wiley, New Je sey, USA, 2004. [43] K. F. Wa nick, P. H. Russe , “Di e en ial o ms and elec omagne ic ield heo y”, P og . Elec omagne . Res. 148 (2014), p. 83-112. [44] C. on Wes enholz, Di e en ial Fo ms in Ma hema ical Physics, S udies in Ma hema ics and i s Applica ions, No h Holland, Ams e dam, Ne he lands, 1978. [45] A. N. Hi ani, K. Kalyana aman, E. B. Vande Zee, “Delaunay Hodge s a ”, Compu . Aided Des. 45 (2013), no. 2, p. 540- 544, Solid and Physical Modeling 2012. [46] S. Mönkölä, J. Rä y, “Disc e e ex e io calculus o pho onic c ys al wa eguides”, In . J. Nume . Me hods Eng. 124 (2023), no. 5, p. 1035-1054. [47] B. Hanouze , M. Sesques, “Abso bing bounda y condi ions o Maxwell’s equa ions”, in Nonlinea Hype bolic P ob- lems: Theo e ical, Applied, and Compu a ional Aspec s, Vieweg+Teubne Ve lag, Wiesbaden, Ge many, 1993, p. 315- 322. [48] J.-P. Be enge , “A pe ec ly ma ched laye o he abso p ion o elec omagne ic wa es”, J. Compu . Phys. 114 (1994), no. 2, p. 185-200. [49] J. Räbinä, S. Mönkölä, T. Rossi, A. Pen ilä, K. Muinonen, “Compa ison o disc e e ex e io calculus and disc e e- dipole app oxima ion o elec omagne ic sca e ing”, J. Quan . Spec osc. Rad. T ans . 146 (2014), p. 417-423. Sanna Mönkölä e al. 19 [50] P. Mullen, P. Mema i, F. de Goes, M. Desb un, “HOT: Hodge-op imized iangula ions”, ACM T ans. G aph. 30 (2011), no. 4, a icle no. 103, p. 1–12. [51] K. Umashanka , A. Ta lo e, “A no el me hod o analyse elec omagne ic sca e ing o complex objec ”, IEEE T ans. Elec omagn. Compa . 24 (1982), no. 4, p. 397-405. [52] A. Ta lo e, K. Umashanka , “Rada c oss sec ion o gene al h ee-dimensional sca e e s”, IEEE T ans. Elec omagn. Compa . 25 (1983), no. 4, p. 433-440. [53] C. F. Boh en, D. R. Hu man, Abso p ion and Sca e ing o Ligh by Small Pa icles, Wiley & Sons, New Yo k, 1983, 53-56 pages. [54] C. Mä zle , “MATLAB Func ions o Mie sca e ing and abso p ion. Ve sion 2”, Tech. Repo 2002–11, Ins i u e o Applied Physics, Uni e si y o Be n, 2002. [55] G. Tu k, M. Le oy, “Zippe ed polygon meshes om ange images”, in P oceedings o he 21s Annual Con e ence on Compu e G aphics and In e ac i e Techniques, SIGGRAPH ’94, ACM, New Yo k, NY, USA, 1994, p. 311-318.