From 1bb5780e58ec95af1f710abc3e2bbfca62ffb26c Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Fri, 11 Jun 2021 13:40:53 +0100 Subject: [PATCH 01/10] Add a method to update mesh tallies in openmc for changed moab mesh --- .../include/executioners/OpenMCExecutioner.h | 2 + openmc/src/executioners/OpenMCExecutioner.C | 47 +++++++++++++++++++ 2 files changed, 49 insertions(+) diff --git a/openmc/include/executioners/OpenMCExecutioner.h b/openmc/include/executioners/OpenMCExecutioner.h index b5f5bb78..c74a20a7 100644 --- a/openmc/include/executioners/OpenMCExecutioner.h +++ b/openmc/include/executioners/OpenMCExecutioner.h @@ -156,6 +156,8 @@ class OpenMCExecutioner : public Transient void updateMaterials(); void updateMaterialDensities(); + // Update mesh tallies + void updateMeshTallies(); // Set up OpenMC tallies void setupTallies(openmc::Filter* filter_ptr); diff --git a/openmc/src/executioners/OpenMCExecutioner.C b/openmc/src/executioners/OpenMCExecutioner.C index 91b94047..938fe080 100644 --- a/openmc/src/executioners/OpenMCExecutioner.C +++ b/openmc/src/executioners/OpenMCExecutioner.C @@ -475,6 +475,8 @@ OpenMCExecutioner::setupTallies(openmc::Filter* filter_ptr) } } + + void OpenMCExecutioner::setupTally(int32_t& tally_id, openmc::Filter* filter_ptr, @@ -865,6 +867,8 @@ OpenMCExecutioner::resetOpenMC() updateMaterials(); + updateMeshTallies(); + return true; } @@ -1073,6 +1077,49 @@ OpenMCExecutioner::updateMaterialDensities() } } +void +OpenMCExecutioner::updateMeshTallies() +{ + // Retrieve the current mesh id + int32_t mesh_id = openmc::model::meshes.back()->id_; + + // Update in place the mesh pointer + openmc::model::meshes.back().reset(new openmc::MOABMesh(moab().moabPtr)); + + // Set mesh ID to what is was before + openmc::model::meshes.back()->id_ = mesh_id; + + // Get a pointer to current mesh filter + openmc::Filter* filter_ptr = openmc::model::tally_filters.back().get(); + + // Upcast pointer type + openmc::MeshFilter* mesh_filter = dynamic_cast(filter_ptr); + + int nBinsBefore = mesh_filter->n_bins(); + + // Update the mesh in the mesh_filter + int32_t mesh_idx = openmc::model::meshes.size() -1; + mesh_filter->set_mesh(mesh_idx); + + int nBinsAfter = mesh_filter->n_bins(); + + // Update strides in tallies if number of bins has changed + if(nBinsAfter != nBinsBefore){ + for(auto tally_pair : tally_ids_to_scores){ + int32_t tally_id = tally_pair.first; + + // Look up tally index + auto tally_it = openmc::model::tally_map.find(tally_id); + if(tally_it == openmc::model::tally_map.end()){ + mooseError("Could not find tally ID in map."); + } + int32_t tally_index = tally_it->second; + + // Re-calcuate strides + openmc::model::tallies.at(tally_index)->set_strides(); + } + } +} bool OpenMCExecutioner::setupCells() From c877b3803108bdacfec2b8d6537dcc40ad83bc5c Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Mon, 14 Jun 2021 12:48:25 +0100 Subject: [PATCH 02/10] Only raise error about incompatible number of powers spanned by binning variable if we are doing log binning --- openmc/src/userobject/MoabUserObject.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/src/userobject/MoabUserObject.C b/openmc/src/userobject/MoabUserObject.C index ae7b95d5..14e72fde 100644 --- a/openmc/src/userobject/MoabUserObject.C +++ b/openmc/src/userobject/MoabUserObject.C @@ -115,7 +115,7 @@ MoabUserObject::MoabUserObject(const InputParameters & parameters) : powMin = int(floor(log10(var_min))); powMax = int(ceil(log10(var_max))); nPow = std::max(powMax-powMin, 1); - if(nPow > nVarBins){ + if(nPow > nVarBins && logscale){ mooseError("Please ensure number of powers for variable is less than the number of bins"); } nMinor = nVarBins/nPow; From f108f84df98bdcec54b70c081c51e8a38763b08d Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Mon, 14 Jun 2021 12:49:55 +0100 Subject: [PATCH 03/10] Add support for differently sourced openmc input xml files in unit tests. --- openmc/unit/include/BasicTest.h | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/openmc/unit/include/BasicTest.h b/openmc/unit/include/BasicTest.h index c8de0d0f..71b1db7d 100644 --- a/openmc/unit/include/BasicTest.h +++ b/openmc/unit/include/BasicTest.h @@ -150,9 +150,16 @@ class InputFileTest : public BasicTest { ASSERT_TRUE(fileExists(newName)); } - void fetchInput(const std::vector& inputFiles){ - for(const auto file : inputFiles){ - fetchInputFile(file); + void fetchInput(const std::vector& inputFilesSrc, + const std::vector& inputFilesDest){ + + ASSERT_EQ(inputFilesSrc.size(),inputFilesDest.size()); + + size_t nFiles = inputFilesSrc.size(); + for(size_t iFile=0; iFile openmcInputXMLFiles; + std::vector openmcInputXMLFilesSrc; std::vector openmcOutputFiles; std::string dagmcFilename; From 3f8411ddf405995a252606574da2f5e61e7755b5 Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Mon, 14 Jun 2021 13:37:55 +0100 Subject: [PATCH 04/10] Set up for new test using a deformed mesh --- unit/Makefile | 2 +- unit/inputs/aurora-thermal-expansion.i | 166 +++++++++++++++++++++++++ unit/inputs/dagmc_cube.h5m | Bin 0 -> 189488 bytes unit/inputs/openmc-thermal-expansion.i | 39 ++++++ unit/inputs/settings-deformed.xml | 23 ++++ unit/src/AuroraAppTest.C | 34 +++++ 6 files changed, 263 insertions(+), 1 deletion(-) create mode 100644 unit/inputs/aurora-thermal-expansion.i create mode 100644 unit/inputs/dagmc_cube.h5m create mode 100644 unit/inputs/openmc-thermal-expansion.i create mode 100644 unit/inputs/settings-deformed.xml diff --git a/unit/Makefile b/unit/Makefile index cb6d68ee..9b8810cd 100644 --- a/unit/Makefile +++ b/unit/Makefile @@ -37,7 +37,7 @@ RDG := no RICHARDS := no SOLID_MECHANICS := no STOCHASTIC_TOOLS := no -TENSOR_MECHANICS := no +TENSOR_MECHANICS := yes XFEM := no POROUS_FLOW := no LEVEL_SET := no diff --git a/unit/inputs/aurora-thermal-expansion.i b/unit/inputs/aurora-thermal-expansion.i new file mode 100644 index 00000000..5223a4b2 --- /dev/null +++ b/unit/inputs/aurora-thermal-expansion.i @@ -0,0 +1,166 @@ +[Mesh] + [cube] + # Length dimensions are cm + type = GeneratedMeshGenerator + dim = 3 + nx = 5 + ny = 5 + nz = 5 + xmax = 1 + ymax = 1 + zmax = 1 + elem_type=TET4 + [] + [cube_id] + type = SubdomainIDGenerator + input = cube + subdomain_id = 1 + [] + displacements = 'disp_x disp_y disp_z' +[] + + + +[Executioner] + type = Transient + solve_type = 'PJFNK' + petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart' + petsc_options_value = 'hypre boomeramg 101' + abort_on_solve_fail=True + num_steps=2 + dt = 0.1 + nl_abs_tol = 1e-10 +[] + +[Variables] + [disp_x] + [] + [disp_y] + [] + [disp_z] + [] +[] + +[AuxVariables] + [temp] + [] + [heating-local] + order = CONSTANT + family = MONOMIAL + [] +[] + +[Functions] + [temp-function] + type = ParsedFunction + value = T0*exp((x-x0)/(rdec)) + vars = 'rdec T0 x0' + vals = '2 1000 1' + [] +[] + +[Kernels] + [TensorMechanics] + displacements = 'disp_x disp_y disp_z' + [] +[] + +[AuxKernels] + [temp-kernel] + type = FunctionAux + variable = temp + function = temp-function + [] +[] + +[BCs] + [./x_bot] + type = DirichletBC + variable = disp_x + boundary = left + value = 0.0 + [../] + [./y_bot] + type = DirichletBC + variable = disp_y + boundary = bottom + value = 0.0 + [../] + [./z_bot] + type = DirichletBC + variable = disp_z + boundary = back + value = 0.0 + [../] +[] + +[Materials] + [copper] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity specific_heat density' + prop_values = '3.97 0.385 8.920' # W/cm*K, J/g-K, g/cm^3 + block = 1 + [] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 1e7 # 100 GPa = 10^7 N / (cm)^2 + poissons_ratio = 0.3 + [../] + [./strain] #We use small deformation mechanics + type = ComputeSmallStrain + displacements = 'disp_x disp_y disp_z' + eigenstrain_names = eigenstrain + [../] + [./stress] + type = ComputeLinearElasticStress + [../] + [./thermal_expansion_strain] + type = ComputeThermalExpansionEigenstrain + stress_free_temperature = 300 + thermal_expansion_coeff = 1.8e-5 # K^-1 + temperature = temp + eigenstrain_name = eigenstrain + [../] +[] + + +[MultiApps] + [openmc] + type = FullSolveMultiApp + app_type = OpenMCApp + execute_on = "timestep_begin" + input_files = "openmc-thermal-expansion.i" + positions = '0 0 0' + library_path = ../openmc/lib + [] +[] + +[Transfers] + [./to_openmc] + type = MoabMeshTransfer + direction = to_multiapp + multi_app = openmc + moabname = "moab" + apptype_from="AuroraApp" + [../] +[] + +[Postprocessors] + [volume_calc_deformed] + type = VolumePostprocessor + use_displaced_mesh = true + [] + [volume_calc_orig] + type = VolumePostprocessor + use_displaced_mesh = false + [] + #[total-heating] + # type = ElementIntegralVariablePostprocessor + # variable = heating-local + # use_displaced_mesh = true + #[] +[] + +[Outputs] + console=false +[] diff --git a/unit/inputs/dagmc_cube.h5m b/unit/inputs/dagmc_cube.h5m new file mode 100644 index 0000000000000000000000000000000000000000..2032062f3c152a8a3076d24306dd29ef9c6abc8c GIT binary patch literal 189488 zcmeF)37lPZdH?a7umvHAB7%q#lxipy1BxP66NEsJEe6ma)+9nCQ3wb@gA0p@EQ-aY z8j1_GRj}F?wN=oz;!=gGtx73Stx~ZHt+r9?Zp(k}{J!rwxliUY38H_AGkukl@AG-K z-*bNFo^xg<@uGPL9yI-aJKe9>|C%)p0l`_nXr$>VR-1m?}%e|~?#3$tpw^!n^G z!&rIx%9Z^;+mG{rdw87E+hjDC{kv`U?>P@`2}he6`Qt?U_l_-_|K$F?hbQ7b9#;3S zeR*L%*DP=4k`qtw|9R<2OHW?9V&&+c`+CJGCoJuY)0Z||pRYXg)KOcXFI#^4%2Q4| zb8O$cK8`x#$OHSI!S>yZP5RzzT_4fQN*~Q@^1mi9>^?Ssi{^1^uJ86aFm>pCY~I1m zf9!mm)$}rjCXXBS1bhxG&JF7NtajViv_H7NY`daq_m}Cpf{%RK(?-FRktel&igoTw zr?t4n?J)M=_BpWm*wfvttr~KlTo}D-db7_;%pzK>+l}5G8u9o zoi&$Gv$a>{8ZEMKu?Wu92=<9#`<=W*Ent@nqC_V2zswStrTcRdHc zk-63VJ16Vgal65V+P^1_ykGUn@u#dUAn7dD4f1aBCt$y9_cb{eBqr*N2uGpci^}2_&<3yhW3zs)PbKdwSHMyd@ zo4~Nof!nfw#Qt(V2R7@?Oo3(X+Y-vHbM*PNaQ}7CfWhPAyg38rb%f`}RR^qt1`Mu? zRk^Nw{?zMiAm~>p{5))4`9AE{y;-AOOzCaW+j``CKPN9)wP^VXz0uFzwx&7Zt=VTi zpdI_4EB#dGG&%C;RlWYVc6S{4315Hh{rR>TOzX`(;l!nVw~=?u{-O7bGnR~0==J8G za?+V6p0c9<&e@l}`DdQC-VZGQ#`hAtV3k_quk!eUHsXX3^(U zr?#}yl5dkKr#`IDd%ZcM|Mi~-=KE&-n1^lsnZ>=4dHnzU>aX+ablU%O|J>01^VUDN zbie1d_0$o4@Bv32*0+xiX5=vX=aHj*KI;jO-Rtp>-FvUTZ~tq|zVCkf?7QFIdu5xA z*6s6s_S}2V{!TN0$HR}@f8dI~^!|PC$*;)?Oio~O0+SP%oWSG+CMPgCfyoI>PGE8Z zlM|Snz~lse6%*)xKWOxwKQI1M@*TJ-Q=8$@dBOLdU$yVAJYdG`d~aY*@9|sp9)DK? zy6>;Fgx&X7roFwb@2W%BUT6LO%GSO1ow+H!El0m^b^6LBr>z|Q&Q*)rcH|op?Rye) zMh(7y;n*Adj^f%?hW~lk@8^Eu2`y$`LMQq@`i1*6N8c|QUU8Caj1#EeN1u^J&dl#= zS+k!9kNg~X#nR(fE`Razm1p*u{pXhld`;9z-i->zksl`cHD_ZS zdHA{e=d4rn0~XI&Rc|u#1Dbxa{?D)4xtcRN^x@A}SKm46&^sre^iADc&sQ@wnx8AS z?CbyEVb57>p4xKEOR5P!U!T$(?ceq9VA6AAm_Yq}osmUdI`?cw=Y#p#KgRLPem>ZDwo2n4?aw=MUAeK<{XTy5#%uL_ zP~WRI@cCe~iJlLx$=|2)`xI4g`18T6-Et0J;#~I6(!KS3u)6B^t{vPb5 z^V)e_mK{6M^Vy6;nq#|$5;S@EjhR4wKI^}K_rKO0)Qu@8eGE zU6{z#=e9q{e_ix#=QaDx_^ml|&~L|lobS7!*;mJJeiF;KDth(B&Ecx}-SYeG5Bz$U zHa+7F?YM9JddCd<&5WM+tY*I^eyiU+=yzT8ZE0*?{AN98&~IV%s^>QQwD?_k;Go~t z(O<~lHy{0Z)A1)QIsNnzd;jN5 z?O9_+-)-zXn~UF*-f@)nMt)Cv+tJ^Nt#R$VcK?3Hk1aTHlZW-2famUdzPsn{X`|nl zJ7MYbmz;6Z%8_r$nY1(3{}-Y9=X>U8!KqWW81+B#q*I=^o;%$!|p5dvX6}9-oCH+JY1FgMEyKGV(u}A?|;M*J`R87pNG?O zU0DC*as4LX^Ke>zn__16QTO>)KMzfMJwFe}cOUwBcy*pP*7x&p&QWchlb?s{cfL0^ zuVJ5u3$u@`_m}l~c;V>NU-!4UJRcr-$bv%-n0MIRqw|UQtA0M*G1rClPafBA0zMDx z`_J&tLz7<5&%^QEhkhPjm-nFc{XCq%u&vYk#QLpz(raUxK>a+Nk%cYGKC(`qhs#e` zH2P19ym-k;XY~I`9J_d5j_Y~!e7O497EqrL51xDE;foGDS-+g2s z-u?IK^Uv*#yj8ca4ga~t)o*P-x7l~#RNYs9Zeb3N{dx3%hB0-j0k>ri|2ftjnaAKi zp!n-MkMX~svL@f>Ti^TTIg8p{T>k@`VNz`L6R7t~_lwopM+V@2=~HbUa1W2}m)9?D z305VK_lNp>8ykAR%+I7pK1ciCd0gN7W%vE%g}E-$*5qM5C*Xcrf9^2+ereKn)&%#< z@!co9U;6tuj@|u1!nVNS_lHg2*7np}w_bY`zpnR(`ty=GOV;js^3Ur#fr*~O{Jyc{ zhVQE@-qZHgbz838(|u)Mtv-J3t|$Mzz7v>eU->;X#|__Ci$B=*)ouA3Nv*T{%D!4N z`Y&Rw?=&awvkL=*!$+R{y%8A*4MrLzUjTB{=Sy?&4czoVD8a} z9C7fXqmMkS-^PZ%Z+d&_zi;+GA7}KFjC}R~o`3&Wx`+Ke|52xL{-fXPU(nzkC$-~- z_SHVC``_2@|6li8)SHa_{!YKI-yZq@b4LEZ&%FDO*K2*BSe)y4By93)112!h=i{{g zFU^g7`J8QP^1r_+6R4k$GqUiRrwvZP=i|Ng^DCb>cbwAdOiMlgdq+p?fB4)!+0f6Q z?sIYK=;xw+y~)VmuIPXMY(Mh4+6 zPCe0kfct`Df8y=|ZKiWZ2Q&QpVT?2bm&54y}9g_?t0-dcX2%aLcN}^}AoWAIuxX z8?5ZtX)kM`3zNa~{K)@4-=YOa9+}@x`$g_s?YY@?(*K<7-?#efT>EpfJ|C~nb&<9v z59>Jr&mpT{+ybu2KI;BGv-^!agG}o!KYh`%B`Z!iX=(D-zb|5+b;tJm6waglv+V4a zV0}ODejNYF&%5$u_H-?=*9&)9chfcMq!@0f#Qe;)1oJA*rG`1g12$UJVl zkIZBI`_|(8oNayYTc^IZE!Xumu=q)_(NAF5ee0Tw2bX?t-?wD1|IS4Hd)D>u8P>mN zQ-6{olKF-*W8kU%i2yHT?c{+kI>v;IlW|MzSA zzbnxHyzskLYya**y^r0IK1bZg{c<0>uFS1EyUu?6)=D-`?K!3XyO_V|&wsY=%^Mw9 z|8EA*82NwyF5*Pz)%~GPd!qYo|9g`CukN^^zn6DM?k{#g%XeRWU+m`4*3sQ>2d{wP z_uJk#+I2UpCCGn&UFR{-P^<_@-&^6DX-8$s2 zeY5kc{j*=Y&bnG%^SN!E^|n4byFOe0x39BacKr_?ZyTFK&BISW^6+yUv|Cs0pWS|_ z{ogSA`J-RoTxY%2u5LTtxz2jo_3tPTUGogu^(%ipzwG>K|LoVUv#wV6hS{J0T>tCa z>#Voh_1XGcv(9?i_3tPTUGoguoiF*j^)^3t>#F^;U%Sq_tXKWEUv3=RcY}R&*Jbsm z2lG|`dYxH^{_1s!v)iA&UaP}yKE0m)2m7hlan}cb>!ckA^Qc$Xb?cDNiO%0#2J_-w ze}mk*ZX3(5&vxwW`fTw&XSzPw2lK1j`Rh~urn^4b2lI5-7rQ=P*X;Ui$2)VqbbYc< zyf?Xv)^1_U7sy(&DePv=C4orYoF|cbzKkYclQDPHvikl zuG7K#>wbQZ@89h&@%rUY#6?i2t=?`$j)1wLYDv z;r`YozdVEMtaYtRe~zltU0rz{c@FB~{`xc*ar$I;s`OT*t-KT2Z z?s}|!wtjD0r|#hS?)qexXRxl%HFNIj);`Ct+nrbW&69RicQ9|)C;LSH&9AE~t{uDc zJ(#Dv9@+Kjx?q>5areirPxguY^(lYtlYKBxxBqc{4j$L6bARNofA#e%o_!*Je)7pP zsAV-L`={3%KYx8r@2*GwGP+*a)on-j!8(Uqhr9=L{w=<%Th}#yUF*7SEWbY6v9s&b zbtgZ6_PQ>1{`!=E(C+iF*6pqf`PH@0*!9_t_II7tN7szjXHrevU2PeWknZYo6wZ^LKx&`&hls zf3RM6ztR8T@z!;o`}N}}A8j4%^x(19wT`aO!F=7i2JP1EdhVXv__aPfuX=vJame$Z z^|}vseYSk;U7zx|^>uv?+PgmO|K{K2XE!(fH$T^-^J!jpj(y*z`s=41>vKuhr@H## zKe+F^KG|!ZP1ny0=f(Bb<>6oZl%L(Y2JPKE_07+*##eje=IGevk)5#cHci4^t*lR`mT2Ij`Gm_*=gLqVs~u(&UNati+7ZVmWQ3j?Gtv# z)~EL=*T45(`_n$C&keVZ^;`3++m79RSlwp#`QM!f_U7kvMSa&vcYa#`KEFDD_VLd{ z*N;B+HNL;=wf1TKU0r$fKd8;s=W&;}))ha@j_ZHWug=Hkd{@^#8I1Gu8=tFnz4m9f zUU~Fmz3lSq(>`J6-?Y!kt}eU$H4nQyH9!A$?B09$H+%OySnZ9wpLF|>oqy}oeW>$i zHxK@`pQi1n?mc(Vzq?;myMAgOcK(j~;oth~-k15c#Ubu{+fqfo|>P3J9h7P{F}Xd-mmt-`rSTc=ijvdcAwuY z`r=>v(I-B>e|I0?Cr&>!&Q3eZ&%gO~?`!4=3(bw^Ygb4^uxdL?)|RX z<*#|z`PcmX)pO)8ewZDvdD!I<@5sL$uN`|%ZgK7p?vKsd>#@6EclTQ!j;*fygZBaX zx;(AE=WYJ_QOAAD=W6G#o^>=WzdY+cFI0bXl5a4N`*XM67GJMFb=kW<<+on>2kl*- zcqp_O4I+x#nlr;b6SyiaJmC%eqh4yFTTweGb~Y zJ_q-Gm!G}%*|g`7u21&5-a&iUr~I|gL3`Jyy!J)QYhTG*pC4Ucw~oEH%hQhbZ@rE> z|5nF3YkvKCFX)~t$L~|#n#Xm4cX?``@(v#DZ~MK=-~76K=Eg2RjkD8^wXS`|&acJa zIrjcj?eeQD4?BN$`T5%??EKrYdta+|`D-3_{xv`U=GR>({F}XgF4z6(JanJu^2*bW z_FL!AKIkWZ*N=Sd*zE&xb)D5d7_WY}zxmhv;_!ML#@XvpJUd;FxP4N2Td&Vm*SWdb zpVnEQq*@$=XB zM|17o*W_>Yy7#qe*N;BiarpbUe)!ABPUFMuwXQgJ+EJdGmz|&ZVRpRcVV6g|BY*pZ zoqs!a?`zdAf6c?rzvk!P{JQrw{>|RKuT^`~KBv3SLw5e{=)avmi@x~Re)QRnH;%m@ z4f=QQYt`O#S63c({_Ofu7iXvOVRrfXiDMts^?k7W9IN@A*I{<~)g9bV-RBv*Jp5_L z!9Kg6uT^`kI~d>9Wp`B9d3V$w-ST(8-_z{w*WKqFKgZ_RJ!i`|nBQ}zBY*E%-S^K1 z^LGA&_U?03y}_f;W1okef18Kpbm?Zk~hrx;_W(-8}94 z+9$hy#5e8w>CQWQ%`<54`jo%-}8$Isur zueCV)%JXgaoFeb|dGzPlDS;vLm(eRl6_gZ}G2H??^CzpE=xt;=5P4%Vsb?bg9h zzw*`fSJU{fK|F`~2^|ug~7{n0L1i+2v_{ zcJFKIvbST+-}S>^-{brH+^%`L?-$fOwa=QLzj}_=C4QJ4uX)(z5%0*qavJAf?eg=( z*=a|4_}eG!gML0IyE!*M`?K2*^3=NQ@~hj9-Tgq_X74^vt9>xO+lTD@TfOe{9cSm~ z*m&0u|F)j)Ivn(0_d0KJ-RDuwL(5b1&s`c^>P2pOb&>M}GTc{B?Bi zYlC@wPIh(WseQ7`uWs|}-lycZ5BXQS{520d|C*maUi0I_?DF%Ihn;qmhd=GePy8@D zUh}ZaBi@mJJ6=2X^R+?l{@Cr8=I?&ps`Pq3uW3i;-}7s?-{o=C?_fXObs%H)AGFu} z-MZwr-_)^QcKO(wpZ$zlm(SZ;m))G)x4b8I>rz)g`fNwvAMNt9>qMWf3-)$Y*Zs1~ z!+-pFs@wIcZtc_lV%HD5e%PB|_c^3b_L^tV-py0~#=AZT?e=f?{f?SnU$swtu@7qJ zzndqnA6kBP+EG9JY5QW(?}o9T<5YW#vroJJ*!j0Qo=-Y|_VMd>?`!(x=P18=!|Zs? z!!D0_NB-^Dy{`@WckgS}E`QC#&cEj8Z(r$$e~a(l*Q#Cqnunc#&CegNb@5?#`T5Di zPCLrOzxla-y8CdmyAId2w>*Ac)9pic^;(_oeNA2VcC7ike)tdeUH$Fvb?y48dD!(? z^YhoAGh=56#bdy6=bbul>m1*3-SO4f=QQYt^ownunc# z&CkCbyZ7Wl|8;-BA_ zy8Au9X78Rqs$IOJJT(8N{n!3;?4G~*%R|d&pRhYN|69j??pE#M9p$0#f&${ao$p*8J+KQ;)+qdp(M0r|S{FV)WmG zy6j{1*!o;EyV~i>X-9tbh&$rLH2bh)i+@l0?~ZM~`tR)gE2ka#)g$hR57X?!j{bc& z|2~<2+i&xBg2w*ccXs}`{~e0_&v=dbv&$=Po03-?JOBQ_IoSVFmsXd3dUG55cjjuI zCx!2uqde-c^JAY=@~XqG5Bd847GnR4eqfWM|7K$Uk57)q)ur*L#ZNrWuC9MS(VVwS z9IhW4e^C6y}@H@02|B&e7}#7H&WA#}A30c${5b^LHK} zmN;Bp8h?2F#N+I{6wc27E8)&ByST51^TXM#OF#Og@!gV#epEF3p@o|de|&!Y#N+Jh z>fd@Foj6=w8lM?I@i_bLg|qX&J=}cQ#oZas4`;V7{pgd%_e>soRy6y9!p(<2{;c?k z$Jy1@zxD2wI9y#C-#dQdarVa*&d&co!p(Mf0*x<8G^e-_P; z{lS?hEskAXpNEGe&kJ(I^+V&Q#7{h4?F-^JKgSo9xayDRf1b4`UY%+`Hu+pftD?nK zfBY5k6VJcekBguEe0GVe{`fiZ6R!^YxrMXy-zMkJ{%03AGn^mJz9{*Z4!i52k{e+SG)ab-@c{9Re$`g@e{93wVxe7`}XZ6uKMGD96#~uRQtK{ zJ156?mALAUzdL^7)nUJ~aCZK?=lt3K?BWgx=ZCYOpZxZleQ2NH_7(l!#QkZB$JzO@ ze<1#}I_&zef9$J2D|O`OkN;)-#N+Jhy+3~H&dQu{`=1~CVZ~1zySnzLefy!(mwx!; zAB~@QyxQ$g`}S{2T=mC45kK+jRQtQ)XWxFR#8rR%hWLqBr`oTKpZ)yV5?B54&&5x? zI_#e>oSpxHIe+#)ySV3u^TXNSoBZ~heMsY1$B(v;*lFB8{EL#Gzx?)z{X^r|mU#Kt z)n%8Tz9#Yf*~L9Q^Eo6(oLxNoGm4)$c6kqtpFA|~`jX$e^-0_3_CJkZm-W%|vC|)m zAN&0DhwGmo`!kE5IClQ}lAm@RxxQ%pGwF+#kDdN>{Me64U$}nwu^(6b#If_&m;Cf+ zQ`hxNUh`>d+@dSG#;Chw~H1u8*5j@0J{wCJvXMANvc6pE!2@-m`v|{L2%E^XJEY zN%0fMetYs+@7If7k$CInCvIgpKb-wPl8^tLMW2y){`|zfESw+C{@vu`zozJyC!Rk) zajy*LhqHe_`S^dY=<^cKpP#r_h4aJNy=U>C9)J4$X#VWt-Vx3ZXSW}=kDvFH*GA*^ z13&hwik~?4-NW_cIhnp9TtDpMJ`&CkXV;JCK+la=CJxsRKlV=+KXL2_hU@2mqCXm~ zAAaI)3g?Hj>&N}pefATH!}Y_D{j0@K9J|j+{rLQ(|2gscVHdY1oFC4vpF`s3`GKAv zyZbe}xUYnNIY*pbJp1j%PaM0v{~AAeXxx0{KPLX}%U_Si<>$wKNAVNK&fk6TxcGl7 zN1Q)D_U{xwaqRq^FZaQ3Cl2S&kNtbaPaHe{@5XOUjx^5ygmCw%hJMsGAC+@ItemJ{++;={i{7*_8t{;Bv#}q$t?Cv}IxheVS!xOI`c5x?$ z^TXNo2x?!Ig=}<>$x#(&8tMo&UF!-+l1R#Nqt; zv7cT1#If^tzT5}TN*vCgAN#9{pE!2@=f%%?r*Zy2%sSnteh}^RA9ugt$9_TT&=*Em zyL_(?=O>O`{-33fpXYc{;&A!-vA?ePiDT#QbI<4GC5glN^J9O1@e{|sMY#3457KW4 zw_bK}tHb%>?AwI%pI-E3(fs*|dq+4woPGOn{tqhpifI1)#JxM5AI|Q+#@~IAzB2Lr z*~MKI&JSmI-{HSk{OR{Z^Jf?L$#8x+yZ!K___<=h@;@J5+ zU+#n3CJyJ%kNv^LPaHe{?c?XX(>VVNvp)B!3yQ{nKX&>x(d=(1-2H_={^t0J$JzC9 zN&H@yBaN#o56=FkQb&G%;yzIP{ye(c<@<0rKY7^Y_qpnGn8xKd7y0#hS?b{Kr!@Yy z_=(5aFE5;(|6hgcpIzKvhx5bP^}%0$_h0-S$wT9R61#Yu{hfug^Z!V=e%Qr*Je(iS zt`Gk5Uzs@E^8k&%H-6%AcK30do&H$j^}{ai)8YJZc7FP~D*iOC?x(`pKUp*`zkP21 zejxt#18)EDWB+{o=^LV}UA`}c^ApD||7YXpIf2IIci!ciA=AKX0ovEPyW z^f#ibUB2&x^ApD||2O0JtsH4wesh-Jd7=L^@$$2a`(8LdoZX!HJ1_K_#Peqt_oHxr zIJ@(~-+7^bka+&=;(i{^4`=^rIRBp%jq`V2%-fvJ2RBcC?3?5{g}y8Gs$ITm@#iOw zU4GAxo*QXgesh-Jd7(WY%Fixtdivpqvzs%2=Y{rs$e&%@1CyU0&hC8hcV1}Ehy2;a z?U4NZaCXm){5?<7IDd0?emu|9xcRf2pXYWOH-C2X@ElI#=E-h-o~voxeAunab25!v z?~a+@Lvq|HM|OJWVyDGDw8U3Cy-WP~)9O4tT>fgO#qp=bJuG?Dt#(=*e_Gs(z?EJ*B^TU5DdD!*K&JSnjhnH`6x?Fmj=JzsE zrk9V!vA@o6PLjc-{vjjLDl-!F02gVXp{B|nYpM?8&N z4=(@K$!p#?jZZK2XF87+tfFo2NaDzFm{~Aw=0~+)e}!YC~^F88sEO;r*ZQVPvh$0@;^9vtp}&^9ZG!~ z*N=D_R}Yte$K)M0f86;NPviQ*XfvPviO#Pvh$0^8Zfq4x2yje2b@X{opjN9UU;V^ZI| za2mHBoW|9|Xh7E4lXBdzX#B~A)A(kE)3|!_(_1HwKThLYmHae5y>J>= zujc>V#90qc<4-C1X#?4bajq3-e zarJQV`=`En;WTbNIE|}^)3|!L^YOIQH=hHF#t)1gr|}0CPUGr{r*}*of1Jj5DEVpJ zyu{PEdbs>^lh=B18b7Gir*Zv=r*ZXg`43LsVe`kGZ}Bv)ADqV3!_}+TC%#$w!p#$> zasA*lt{yIaUh0?^PUF^t)3|y#jjM+{ABUvA`8>U7{Lt8O8sE8a8dpy|{fNZz$7y`m zlAp%SOFWIMhs%Fh@>&m0l{}IVMZ2q|OEuO~pgVVTrxO(;a#LWvg zPn^c}gVVTrxcDQpe)Ga<+F87emVcoC>mc7J5J+|DxAjE6Hot6;`rk< zzDLPVo+f)#;pgZarJN-R}X(u>OMEeML8}m8edX4jX$|? z8dpz#`l*TIkJI>5N`4ymz9XK-)x+g~Uh-NGPUFYN52tbc;54otF8>M1J8b^=()fv| zasA*lt{$%5^OLvEhsMoQJdNuIr*ZXg@h7IfdEqo}Jvfc4hts%vxbv~B*S^eWdC~X_ zV#jIxz`|)2C-1QN zCFWfwF8rKg_h10n8;54otPUGs~&c`X`{GVDh z{-W4%8b7pf8dpy|eMI8;<1~JF$xq|vC7#CB!{t9Md94Sh@zYCv8rP3_8dndOe`WI8 zpE!*>-#Cry2d8oMaQR=D`gJ}uZl2<4Tt7ICtA~p}BlXP-r*Z4SXf2PviO#Pvh$0 z^1m#3hs_^%zQxnHesCIB4_EKB)UWfQaq|>UU;%TwRHa2mHBoW|9|XfPV!n0PUGj6`ZTT| z@ieX;F8?c&ci8-K=UY6D>j$TC^>FphO#M0^8aGezG_D_<#?`~cpO^aPh10n8;54ot zPUGs~i&FPhIi8>6tBb}jD4fQZ6i(yn$xokl|AonGJvfbD zRO-{Xe#Fzbdbs=-C-1QNfz4E>r&r*e!pn^(%5kte?j3iuAX@M7+j4|hKPpq&3V7mdFqcAUmvR5*>RC!SuJIQ}?|pI-9QxOs`EarJQd zS0}Ia;52?&sZZnj5l`dl;qw1s^4gy`jXU2sjq3-earJQdU!VGQJ~VEg;%QtzIE|}^ zi@!Ye%?qb->%nPUJ)FkX!(WuTZ_Uws-c~gJ_SkV6UsX7bt0$g*S>pKPH2%_(pT^Bg zJdLY|%m0q#wH}U=KT6(V^T(ZU@ieXYEo%j-9OIJeEy_p{9Uo*G=6sBG_Ia_`jv^}kJI?M zB|nXumv|ai510Sl$!k40jbB;n)3|=b)3|!L{8uIKu=(Rxl=?KTANgrqJzTwaByXJ$ zjhm-<8rKg_)TeR%$WPf!Q#IC-rHr}1k_eHzz~cp6s^m;c)29X5a5`4&&( z`oU>jJzTwded6YYnjJzV@pvVQZzY212n8dndearJQL#$-^7m7 z_*)96arMO0f0#J_IE`Oc^3%9^iKlV(aQXi>d94Sh@sE}IG_D`l|8>c0f8sRm zeB(5(ADqV3!{z^Q>euP!6Z(cZ! zTMtg->fto59`1ZRB=ybbUy86Q^w_%4OhxO(F0uP2T_PUE+i{4{P};%Qtx zT>ft)@2(|3jsJV$G_D`{XCueD}g>Ts`^eI}^trr}6(N`Dxs| z#M8KXxcvW_yw-!$__s@a8rP3_8dndO|2xS$Z2q|OEuO~pgVVTrxO(;a#LWvgPn^c} zgVVTrxcL9d`ppZcaqGcpTs@q|)x({SCzkX7-$mp5#E#SW;|iy7^~BR_62~8>@$Z)W zG;UtvX8ZPaj~bspO~en+vCL^=kefCC++q8vk+0PviO#Pvh2u%m0()9X5abr=>oP>qmYXR}WY3 zXUSXVL*qA>`ZTT|`Dt7|T>Q_I*Sv5Vw;r6v)x&9AJ>2=2m-^=OzeVE*#g5bX4u#XW zdgAG?C5}H%F!nQ(EGul>9XAe2b@X z{opjN9xnfnQn${B#?4bajq3-earJQVQ?h>Z!fD)ka2i(+r*ZZ09kT8tb2Oi+>4(M- zj~`CsyB1F4>WQbnnK=G9jeDNMY23VU8dndOf0HS#0qenOeAD>hG_D_<#?`~+-z<3# zD)Xap=UY6D>j$TC^>F$BFLmpDXxu!-)3|o=b*(hrRv9Y37L_b8mk)e}#DCvp678uvVh)3|xzG_D>l|Fo>vdT<)wGJZIX>j$TC z^>F#`m%N9U`O&!ZEuO~pgVVTrxcv3{#LWvgPn^c}gVVTrxcIHIe)Ga<+F9o_cC|$**g8u_;K;WX?&l;X3l|NYCn@uM>r-1)|7Tt7ICtB1>9uTR{(aP!1zTt7ICtA~rjJzV|=m3ialk8fY<)3|=b)3|!L{0~mvIv*N0Pw_OaADqV3!^Q8A z`sRhxxb@&Pt{zU~>fz4EitJnS*)jdl_(}1@Y25Q2PUGs~bbUX=#p8FC`ZR7{;%Qtx zT>gh-z1D-%_)hV|Xeyo;54otE`Pl~ar45>6Q^2;?Bm2sH9#-}PetP_H8uxsM)3|y#UEj}e@wn$XoW{)y zr*ZXg`5&J3S`SX+yTlKtasA*lt{yJ`j52TB`M{lToW}Ko)3|!L{Pp_8%?meAoW}Ko z)3|!L_+7Jp^TKJ|dT<(752tbUaOdOY+1KXth_WB>m&OmLanE-+jjM;#_5BPNk9(fO zY23VU8dndO|B+d*_24wVTl{bu*AGtP>f!SLR+%^MeBjPEPUHH)Xfz4EdD++I^XT+LOt{yJ`W3pcB!D)QY_~A6JADqV3!{vW$nKy3!_^eW&#`Pnf z#?`~+-z$0Rd}!P}#nZTca2i(+7r%Gvn-@;w)`QcydN_@%hdUo{%Dy$9$E6<{e?$Cm z8uxsM)3|y#UEj}e@wn$XoW{3H-#Cq{hs*!?#90qcl{}al*ar4KW zZ=A;UgVVTrxcv3{#LWxeuB@NN^&_6f)x*UK&1c{A zL*tjl52ta@cQ}o!htu`_3>S}kp2KO}yl@&<510Q*S+DirG`?T_a2nSSPUGs~@;|xE z8#jO4`NnBnKRAu6hs$5DPu#rloyz)YTtDJzTs>U;?@|qV;7+j4|hI3l6`AFN2DJbzb1Y-jeEYsX3l|1*=pR#@$ z*N=D_R}UBetmHK>oW`vOr*ZXg8dnc@KCaKcHlKy*hsHl0Kb*#|D4fRC6HnLoGh96G zc@C#>^TKIdJzV}{vtH}LY5dvo!)aVUIE|}^%YR&%H*Ws8^NrKEesCIB50}4QpSXG9 z`<3<6xPHXbxO%ww=OnLr;WTbNIE|}^)3|!L^YM?p_GLc5mwssc^YO!J-1`+yjJzW09W!|{?U= zdVS*Nh0iYQr*Zv=r*ZXg@k^4|yl@(~9-PM2!)aVS{Kj7Ux;4j}b9`R%nPUJ)FkX!<~=YvTx1jr1V4M z{{9$F<6lW0oW|9|>625BKThK-;)m0?dEqp!9xneW$!k40jh|ZT)3|=b)3|!L{4Yvg z^Tui1`NnBnKRAu6hs$5DPu#q4^TcUfKRAu6hl@Wg>o+f)#;pgZarJN-R}XhSzLkAz zKBt%cfd6~^a2mh8a2i)nJiRjY_~SHwM*MIZH!qyV)x+g~aq?OZPU9~r^=Vu`;%Qtx zT>e$bJ8b^A^DUmn^@G#6dboP^`ozr(H&2|#^@G#6dbs#Avwri!Y212n8dndearJQL z)x+tRr5=Br#$O&koW{)yr*ZXg`CpN|)`QddS*1RW z>qk6|tB1>fcJdCJKkj^sr*Zw@G_D@5UcEkX^TN#&r*Zw@G_D>l{+z7eyl@(~9-PM2 z!)aVS-1+!H_Oe)jul3+G zetxM>fto59`1eQ=h?SE&GEwYL*qX#oW^e|oW|9YpRVs`xOn`;^o7&7_Z^(Z)x+h# zD0!_1r}2wReHzz~cp6s^m;W_o-njYWuPyayTtDJzTs>U=OOm(FhsKwe`ZTT|`Dt7| zT>R^j*Sv5Vw;r6v)x&9AJ>2=2lKSTJ`$gk-Wq;!|?)eUjJzV}blzHRM2kw01G_D_<#?`~+zbN(Vd}!P}#nZTc za2i(+7yri8H!qyVtp}%Z^>7+j4|hH`Pkr-wQ_=XQvEwxE`3|RX^>DhrpW)(h&vQ79 zn-@;w>f!SLLDpqGIE}x#te?j9Bc8_9!{vWVnK$ly;LbNrU=OH;qjhsMoQ zJdNuIr*ZXg@vBq+C1t)eZav~@Ts@q|)x({Sty0%~E-M<}GIpHCJ>TIpt{zU;_cL5P z?s*QUar44yTs>U=Kg_zU2dDAN%lc_tKjLXzJzV~`mU-jO2kw01G_D_<#?`~+uh%DT zUbuPUG_D_<#?`~czb)%GFPz4$2d8oMa2i(+cRsc)=l|_R@8C3kR_Tw%)e|rOI}&F-IE`OX^3%9}#M8KXxcq;Vy!IzfU= zdVS*Ng`ZW{PviO#Pvh$0;@_FP=7rO^_24wF9!}%x;oe7fNd4_|{NtkW?Fy%H&v)W! zTs@qAP4e=`Y25Q1PUGHpa2i(+m;X;v-+FKwe^*&Qjq67|jjM;t|L)`+Hh=ueQlG~4 zBR`F+hpTs0^49s#xOs}FasA*lt{yJ_J*jVAIE`BmPUGs~G_D@*eC(Y1=JTgT<2%KU z)41n5oW|9|>H2<#i^o0B;WU16>5scLy@38sf&bN3P*AGtP>f!3u>l5EJec^8|>!)%3h^KM&aPik9uX*7#Zap}StB2FL zdbsnkN9vo;wMFB*$BxtZHigr;dgAFnN*sTj#=Y<0G;UrvjjM;t|B=+U9-PKMTGmJ7 z`Vmj#>f!SLb@C3IKkj^sr*Zw@G_D@5UcEl?Ez=i%d09V=>qk6|tA~sKo8&bwoW`vO zr*ZXg8dneBK6UrbaaN9hTQt6B;WWNo;WVzE{PepM#~-J0&vQ79d*8unTs>U=kEOo# z;52?+Ss#t-M?8(Ihs*!*DW1mlgVVTrxcE<{ zzIov^Zap}StB2FLdbsoP#MC#RzbhKwCw831cPgC5)e}$G_cL5P?s*QU@pqN_G_Ia_ z`9GaF>%nRKGbKNb>qk6|tB1>feew>QKkj^sr*Zw@G_D@5UcEkX^TOX%)=%U55l`dl z;o@&dUh~3f+F9ocT?Yd{=R5@zu0ja-=%OGS5G|sfyD90Y25Q1PUGf< z)3|!L{GUyI>%nRKb7g%rt{?FIC-1QNFdGB(HhlG;TdOjjM;#xO(`k)IB`MLvy^fX#9}EX?)MZ zX5nIlKThMG=WrVLzJt@adbs@mlKR$z)A*Om`eh^l@38sf|61zP zxPIiParJQZzM8ysJ~VEg;%QtzIE|}^i~m~cn-@;w)`QcydN_@%hdUn&Qr~=TD;l34 zJ5J;K6i(yniKpxP87?08JcrZxb)`Oyt0!Lme@mS8;57dAlAp%)Bc8_9!{xs{d56s( zcfQ5bxPEXNR}WXOUZ1#m;n$V*)3|=b)3|!L_-`bydEqo}Jvfc4hts%vxbyL>)Hk1h zFB(54cAUodE1bsF6Hk9War|)__dJKwxOw3;t{yJ`H&fqwa2o$sSs#t-M?8(Ihs%FQ z@(!Cn?tF`g2e82RCf3B>b#`Pnf#?`~c|5x&w7f$2WgVVTr zIE|}^J0D9@-+cbNXnaxZIE^1zIE||(p1wJ8{BauhzJt@adEqp!9xnfy)VCg-#=l$E zN8|btPvh$0@_#RRhs_^%zQxnHesCIB4_B{VpZM(bg@3WEpT_kgp2pR~#eY9}%?qb- z>%nPUJ)FkX!eIM>Fz&N#0@e$DMESG_D_<#?`~stJf!PUbuPUG_D_< z#?`~cZ<_U+7f$2WgVVTrIE|}^J0E9e|C!HbWk29Ai62hmXB1B3>WQZ}Pd)xPjc*Y@ zoW{)yr*ZXg`KKkX_24wVWvNf&`Vmj#>f!R=FL{T}A9udR)3|b|qul3+GzC)=`g2+`Mq}#A#eVIE|}^i+@PgZ(cZ!TMtg->fto59`1eQb=kMC&2gvnL*o}0PU9yQ zPUGsyPuKS|Ts;1hlAp%C?}(>y^>Fzgn)O-_PUAbr52tbc;54otF8{;Iym9l#A71Ly zxPHXbxO%w!yCiR&4~?6rcpBFaPUGs~;%B73dEqo}Jvfc4hts%vxbyMG>?`xxHT}@| z>*I&hxaT{Z#?`~=`hJFs$34&CG;UrvjjM;t|A?&DdT<(lWc+X%*AGtP>f!S5R_2YH zKkj_vG_D_<#?`~+uh%DTUbuPUG_D_<#?`~c|5ny-UO0_g4^HFi;WVxu?tHAyzBZpn zmHmLfIes{ed%nYITs@qw?`ODp-18hxU=N2Gq84~?6rcpBFaPUGs~;(t5!w=DCeaqAIJ$M)7#vcH2<#i^sk1;55E%`o?KoJzW015@$U)jqhFZ z)3|=b)3|!L{EsX1#+?t``NnBnKRAu6hs$5DPu#rlZOi&;TtDJzTs>U;Q9^arJPzzMtXZanEx&jeFn0XfV3{{={vO*G zN0#-|xPHXbxO%wwxyfr@IE`BmPUGs~G_D@*d|Z=tn$JP$hsHk?Kb*!r-{CZ_9!}Ty zGh96Gc@C#>^TKIdJzV~SvtH}LX?$M%a2nSSPUGs~@*h&>jhjF2eB(5(ADqV3!{x8n zC;sM~Z+vE1KaJ~0JdLY|i+_6Znio#v)`QcydN_@%hdUpClXaTUq3MUlKN>%r#=T$R zG_D>_*Y`7AJnnr5r*ZSbX3f>%nRK@c7|0t{l|1rs1=R@P>DW1mlgVVTrxcFzLzIov^Zap}StB2FLdbsm(L-w`# zJS+Xs_-EpW)41n5oW|9|>H2<#i^o0B;WYll^o`TFdbs=x6K6d*jUQX`)3|=b)3|!L z{Le1)#?2phzHu7Y4^HFi;quq(6E`pXiDmsXt{?F#5+I*gqerWu2@xy7{^Bqp(>fv;KKf}f2p675HH!qyV)x+igy{y-Ia2kJZ{BRoA z4^HFi;qosk^Ty2|cfN5N*AGtP>f!R&>l1%p&Nu$lvVI!Zk9Zna4;Q~UdCd!_aqGcp zTs@q|)x({S8?#RHS(1Ke{EP9!Y25Q2PUGs~bbUX=#p9mma2hu+oW|9|<$qq*YdtuP zA0I!Q#`S~KxO%w!CzN^P=8rqyIF0KEr*ZXg`RnzGzdz?2KcK9i#`Pnf#?`~cFHK(a z!fD)ka2i(+r*ZXg=i`>F(|n$verWuk9@_bZ&n)x+ufeuj(3z3<>OZeBQztB1>f zV%BRtIE^ohA5P=?!D(DQT>j-{-njYW&Noiu`oU>jJzV~Jec~U=`Nrp!_0zb1#M8KX zxcC<&uX*7#Zap}StB2FLdbszIuVtP8n&S)84~>7ha2oe~C!WUD!|D2dhKt8N&*3!g zeFvv;^>F!5%6hE_r}2~Hhts%za2i(+mw!c>H*Ws;DWyJ*>qk6|tB1>fYVy|k(71Vu zr*Zw@G_D>l{za*8UO0_g4^HFi;WVxu?tFYB``Ua?OFuOJ_4wg5?)eUqk6|tB1>fMwvHm{OiZ*Sv5Vw;r6v)x&9AJ>2=YGyB?nUXp%j{9EzEY25Q2PUGs~ zbbUX=#p9mma2hu+oW|9|9uTT7QIp6qm%lc_tKjLXzJzV@b$!lIXjav^+ zjej?OIE{P1!f9MRoUZR@xOm+A4o>6dh10ltxcslodaVbi@$=$`)3|g4};$O`9#*Z)Sr*Zv=r*ZXg@#iP6dEqo}Jvfc4hts%v_>Z!_ z|IP7dIlemm(73-phST^DOFWIMC!W3__4wm7eqsD@8o#J;8dp!e{1+$AdT<(lP03H= z`Vmj#>f!RgHhIk(r}0ZleHzz~cp6s^m;ZIiTjxXL<|&^3f9>3R$li4s$ML5*opWZ> zvQZ2XhA@H=NF)}_V1y%(YOd7~Mq@1^8e0!ekRXKUuar(|mR2@9sX6CiW?E?nt=Y-) ztR2i!E3>kLht#a6eO?hC_apkFKZ5T-Zr}I&cU_3?&qI}cOD+|@_6;{PriNe zm^TlPd3ikMMPagBm!((0^um0@JZ(ls- z&BJ3}9*=o>y!%*G_y3V%eqr9{F~6WZ=H;o!=cLZ_Jmx=I{V{J}>M<{m*Z*VbYaSl+ z=T<)Eoku<9Qz|A}ILao*=Ke@J=E%Ttd(nL5w&m_NVzW8S{hV_qJw|EJQ|JUr%i zR6gdNM?L1{@%mqozJ5M=%)4(M^UlL#ULLQ1e?EEp;_Z{iyz}swm&dFBbmq4&9`okm zF)xqDygc50yr}N~GsXOpyw79)@bZ|KryhScb)M%j|GDaqdHYh2d3n74pHE-&@R+}_ z@-go`>M<{m*Z-pQ_4CPN-hK0!cOD+|@_7CG^U2#6Z=XEoorlM~JYN09ncu#6%$tYD zygVNB@_4_G9GUzh^7@5herb8k`+lb$^YVDyzn}5y`S(>m=Ka2-9`o{e{Vz#h^YED8 zS^1cE9`%@)$LoJ-`cCbi-&OgTcOLyQFOQdZS^D;VFn>nnW8Qi6$GkjV{pIOvUp(f` z!((0^k9m2#`#383_H#uszbx@tC(S9`o{e{jW^Ed3en4 zuK6(UJnAtokJtaI^qty2@4nSz-g$V;%j4x;lKkEe=Iv8G=ADPfygXk07n5&aJm$^A zV_qJQd3n71cuDf@=jvj9dEVzSzpy;!<*CQ#rOxv_=6#>zF>hZy=H>DFe<}Iq;W2+r z&4+pCQIC0fy#Cjw@6`Tz_pKiD&ckC~9xw08s!_C2m&d!0m(~4$y_i2b@AH^nQXcd2)Z>d%=XoCU7gm4F?I%$tYT|AzGS^T}i0ee;-i9v<`Zc>VkH$=esdqvprF^QgzXJYM~c>1$s+=FP)n zULKEmdHfN{ePv#c&Fi;{`D4mserb8k%hMlUnmW((nD>2-$Na9!$GklC`hPoh=HW5F zr}|^wdDLTG9CI?dSW&{K~w~V}4nA%*#`c`}Z?mJ@5M*kNKUIk9m3O_5VTY%)?{; z=IW1m=TVP&dA$BVOy8;f^X^+c=ADPfygXiBe?EEp;&;~kn0Fran3u<^|55td7ms=K z@R*m!V_qKbK3<)C`?;lV z@t8Lck9m1K=H>CnCU;$4*XH%UVt!3|%pX%8^YZk^H>J+=Jm!6$<1z2|9Uk-Yc>RBo zeDm;_|7FdGdFN4&d3n74_owgF{`m(gAM?(mKj!7}@*Ygz-Vf&OQ$6OLhsV4;Uj46< zZ(ls-&BJ3}9*=o>y!$vQ`S$bcVt#$z=P|#sJm%%8$Nl>mub%gPj>r6-%E!Ds_4@xN zb>`tQ|4{YEyz{8XygXk2-=^=>{(1MU9`nw_V_qIFuRou>eerv0e#|?Mdd$n?)jyoR z_Qhk~JUr&*@tBv#yN@>~-+q2q%x}p1Jmy!G$GklC`1aI!p2xiJb3Eqli^se?UjN@G z-#k3#|4{Q`-g(qxULLRiBk4P}f8Kqo$Gr3Kn3u=P>(3{@GUwuNuK6+VJnAtok5~W4 z^tCS@^XB0(FOSE(Jl=h5PQLy8shHoC_j$~(E{}P6>hV3P^E{7v-{*MD+ZT^{dA$CA zPQH0~%>Sk4!@TpT$GkjV|3}kzYX7|ZR*!k-;W00dm)D<9epSxJ-&XTu-g(qxULLRh zvGlbs9`okmF)xqDygc50ygm8$^LR1;w!F_{{)Fx?)nD_e*k9qszF)xqT z|B2+AhsXS1Yd*|7k9y3@RZNKjxiBJ?7={ z>Yqwq`{FTg9v<`Zc+AV=*CzL!dEJ`Vr;GV5N=XuQgKF4F;?>ju^ zhk55wk9m2#{xkc{{PW+`{`u$RIrW%#9v<`ZczN^ExA%j2`&5s4 z=ixCgk5|86^6iVqym@%c%i}RGk9QxZ)cw!TdF+cnInVK!Utb>c^3>!0{ft-7`##5G z-oALu%j5OmKl7T0$NT|#j>o+7@R*m!>%XA(&D%fkzIn_$5080yy#D?9MPagBm z!((0^ul~TyZ(ls-&BJ3}9*=o>y!&`pes1k&Vf{Str{*~x^Bc-zUY>gVNa{S#W8U{U z9`p9aV_qJw|3R78JUr$XIj>o*;uXxPMvo+7@R*m!>;J;qH}5`p_swJ8d3emr) zC%-B8!rLd0dFSCVFOOHhB=g%Bk9qU(n3u<6ULL>Z(7nHRy`f%D46n=U8}oY3{(GOl zdePoj^~-bqH)UV?otEbh-uw8>>BaeTSu@Ynzk~98+C9sA-)kQK*K^*VPt>daIZoM` z8TCB&AYPC+uzB=fSl)B|>$6vTwD+}n-yZuO=H{AP`CPw$|Nl9zRSp~*7#kQH7#kQH z7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH z7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH7#kQH z7#kQH7#kQH7#kQH7#kQH7#kQH_`lh}+~08V-`w9*@!#CvTk-c={C$=!8&7`4rmfqy zZJzyy+{eE;fBfvmQ#PNvarQ0y@cR9?fB)t9*{$2R?W3Q!)@^;q_AQ&yvtLi#KKs^< zo2DQ5U-$n{`LE00zge+r=|lIu{Kd!NGalZ3{$n@pTIuiP^xws?M-%(|IaB<1GS?TT literal 0 HcmV?d00001 diff --git a/unit/inputs/openmc-thermal-expansion.i b/unit/inputs/openmc-thermal-expansion.i new file mode 100644 index 00000000..9c4a2b3e --- /dev/null +++ b/unit/inputs/openmc-thermal-expansion.i @@ -0,0 +1,39 @@ +[Mesh] + type = GeneratedMesh + dim = 1 + nx = 1 +[] + +[Problem] + type = OpenMCProblem +[] + +[Executioner] + type = OpenMCExecutioner + variables = 'heating-local' + err_variables = '' + score_names = 'heating-local' + tally_ids = '2' + neutron_source = 1.0e18 + launch_threads=false + n_threads=1 + no_scaling = true +[] + +[UserObjects] + [moab] + type = MoabUserObject + bin_varname = "temp" + material_names = 'copper' + material_openmc_names = 'copper' + #output_skins = true + length_scale = 1.0 + graveyard_scale_inner = 1.05 + n_bins = 1 + var_max = 1200. + [] +[] + +[Outputs] + console=false +[] diff --git a/unit/inputs/settings-deformed.xml b/unit/inputs/settings-deformed.xml new file mode 100644 index 00000000..6aaa25a1 --- /dev/null +++ b/unit/inputs/settings-deformed.xml @@ -0,0 +1,23 @@ + + + fixed source + 1000 + 10 + 2 + + + 0.5 0.5 0.5 + + + + false + true + + + + + false + 0 + interpolation + true + diff --git a/unit/src/AuroraAppTest.C b/unit/src/AuroraAppTest.C index 5cb1189e..9012da7e 100644 --- a/unit/src/AuroraAppTest.C +++ b/unit/src/AuroraAppTest.C @@ -72,6 +72,35 @@ protected: }; +class DeformedMeshTest : public AuroraAppRunTest{ +protected: + DeformedMeshTest() : AuroraAppRunTest("aurora-thermal-expansion.i") { + + // Override source openmc input names + openmcInputXMLFilesSrc.at(0) = "settings-deformed.xml"; + + }; + + virtual void SetUp() override { + + // Base class setup. + AuroraAppRunTest::SetUp(); + ASSERT_FALSE(appIsNull); + + // Get the current dagmc file + fetchInputFile("dagmc_cube.h5m",dagmcFilename); + + // Should run without error + ASSERT_NO_THROW(app->run()); + + // Get the executioner + Executioner* executionerPtr = app->getExecutioner(); + ASSERT_NE(executionerPtr,nullptr); + } + +}; + + TEST_F(FullRunTest, UWUW) { ASSERT_FALSE(appIsNull); @@ -87,3 +116,8 @@ TEST_F(FullRunTest, Legacy) std::string dagFile = "dagmc_legacy.h5m"; checkFullRun(dagFile); } + + +TEST_F(DeformedMeshTest, SetUp) +{ +} From 4199fd395166683998f3011258775a47c0572a85 Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Mon, 14 Jun 2021 16:10:21 +0100 Subject: [PATCH 05/10] Ensure we use the displaced mesh if there is one in FunctionUserObject; add ability to set a tolerance to the underlying pointlocator --- include/userobject/FunctionUserObject.h | 7 +++++++ src/userobject/FunctionUserObject.C | 17 +++++++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/include/userobject/FunctionUserObject.h b/include/userobject/FunctionUserObject.h index 261be3cb..04f5858b 100644 --- a/include/userobject/FunctionUserObject.h +++ b/include/userobject/FunctionUserObject.h @@ -2,6 +2,7 @@ // MOOSE includes #include "GeneralUserObject.h" +#include "DisplacedProblem.h" // libMesh includes #include "libmesh/mesh_function.h" @@ -34,9 +35,15 @@ class FunctionUserObject : public GeneralUserObject private: + // Reference to the mesh + MooseMesh& mesh(); + // Name of variable we want to turn into a function std::string _var_name; + // Tolerance to pass to point locator + double tolerance; + // Pointer to the libMesh system containing our variable System* sysPtr; diff --git a/src/userobject/FunctionUserObject.C b/src/userobject/FunctionUserObject.C index d0985dfc..a53c647b 100644 --- a/src/userobject/FunctionUserObject.C +++ b/src/userobject/FunctionUserObject.C @@ -11,12 +11,15 @@ validParams() params.addRequiredParam( "variable", "Variable name to evaluate the mesh function on."); + params.addParam("tolerance", 1e-9,"Set point locator tolerance."); + return params; } FunctionUserObject::FunctionUserObject(const InputParameters & parameters) : GeneralUserObject(parameters), - _var_name(getParam("variable")) + _var_name(getParam("variable")), + tolerance(getParam("tolerance")) {} void @@ -24,7 +27,8 @@ FunctionUserObject::execute() { // Fetch a pointer to the point locator object - point_locator_ptr = _fe_problem.mesh().getPointLocator(); + point_locator_ptr = mesh().getPointLocator(); + point_locator_ptr->set_close_to_point_tol(tolerance); // Fetch a pointer to the system containing named variable sysPtr = &_fe_problem.getSystem(_var_name); @@ -69,3 +73,12 @@ FunctionUserObject::value (const Point &p) const throw std::logic_error("Point locator is null in FunctionUserObject."); } } + +MooseMesh& +FunctionUserObject::mesh() +{ + if(_fe_problem.haveDisplaced()){ + return _fe_problem.getDisplacedProblem()->mesh(); + } + return _fe_problem.mesh(); +} From 89ca08164393e9d95fa22e0b6f0fe78a65b0bf5e Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Mon, 14 Jun 2021 16:14:34 +0100 Subject: [PATCH 06/10] Add integration test that openmc tallies are succesfully updated when mesh deforms --- unit/Makefile | 2 +- unit/inputs/aurora-thermal-expansion.i | 22 ++++--- unit/inputs/openmc-thermal-expansion.i | 2 +- unit/src/AuroraAppTest.C | 82 ++++++++++++++++++++++++-- 4 files changed, 94 insertions(+), 14 deletions(-) diff --git a/unit/Makefile b/unit/Makefile index 9b8810cd..b91261c0 100644 --- a/unit/Makefile +++ b/unit/Makefile @@ -68,7 +68,7 @@ OPENMC_APP_LIB = -Wl,-rpath,$(OPENMC_APP_DIR)/lib -L$(OPENMC_APP_DIR)/lib -l$(OP include $(OPENMC_APP_DIR)/config.inc -ADDITIONAL_INCLUDES += $(OPENMC_APP_INC) $(OPENMC_INC) $(MOAB_INC) +ADDITIONAL_INCLUDES += $(OPENMC_APP_INC) $(OPENMC_INC) $(MOAB_INC) $(HDF5_INC) EXTERNAL_FLAGS += $(OPENMC_APP_LIB) $(OPENMC_LIB) $(MOAB_LIB) diff --git a/unit/inputs/aurora-thermal-expansion.i b/unit/inputs/aurora-thermal-expansion.i index 5223a4b2..b539ea79 100644 --- a/unit/inputs/aurora-thermal-expansion.i +++ b/unit/inputs/aurora-thermal-expansion.i @@ -20,7 +20,7 @@ [] - + [Executioner] type = Transient solve_type = 'PJFNK' @@ -38,9 +38,9 @@ [disp_y] [] [disp_z] - [] -[] - + [] +[] + [AuxVariables] [temp] [] @@ -49,7 +49,7 @@ family = MONOMIAL [] [] - + [Functions] [temp-function] type = ParsedFunction @@ -143,7 +143,15 @@ moabname = "moab" apptype_from="AuroraApp" [../] -[] +[] + +[UserObjects] + [uo-heating-function] + type = FunctionUserObject + variable = heating-local + execute_on = "initial" + [] +[] [Postprocessors] [volume_calc_deformed] @@ -159,7 +167,7 @@ # variable = heating-local # use_displaced_mesh = true #[] -[] +[] [Outputs] console=false diff --git a/unit/inputs/openmc-thermal-expansion.i b/unit/inputs/openmc-thermal-expansion.i index 9c4a2b3e..b392abe8 100644 --- a/unit/inputs/openmc-thermal-expansion.i +++ b/unit/inputs/openmc-thermal-expansion.i @@ -13,7 +13,7 @@ variables = 'heating-local' err_variables = '' score_names = 'heating-local' - tally_ids = '2' + tally_ids = '1' neutron_source = 1.0e18 launch_threads=false n_threads=1 diff --git a/unit/src/AuroraAppTest.C b/unit/src/AuroraAppTest.C index 9012da7e..13a8d31e 100644 --- a/unit/src/AuroraAppTest.C +++ b/unit/src/AuroraAppTest.C @@ -12,6 +12,13 @@ #include "AuroraAppTest.h" #include "Executioner.h" +#include "DisplacedProblem.h" +#include "FunctionUserObject.h" + +#include "openmc/position.h" +#include "openmc/mesh.h" +#include "openmc/tallies/tally.h" + TEST_F(AuroraAppBasicTest, registryTest) { @@ -93,11 +100,21 @@ protected: // Should run without error ASSERT_NO_THROW(app->run()); - // Get the executioner - Executioner* executionerPtr = app->getExecutioner(); - ASSERT_NE(executionerPtr,nullptr); + // Get the FE problem and mesh + FEProblemBase& problem = app->getExecutioner()->feProblem(); + if(problem.haveDisplaced()){ + meshPtr = &(problem.getDisplacedProblem()->mesh().getMesh()); + } + ASSERT_NE(meshPtr,nullptr); + + // Get the function user object + functionUOPtr = &(problem.getUserObject("uo-heating-function")); + } + MeshBase* meshPtr; + FunctionUserObject* functionUOPtr; + }; @@ -117,7 +134,62 @@ TEST_F(FullRunTest, Legacy) checkFullRun(dagFile); } - -TEST_F(DeformedMeshTest, SetUp) +// This test locates a point on the deformed mesh that is outside the bounds of +// the original mesh and checks that openmc can locate it, and get tally value +TEST_F(DeformedMeshTest, CheckExtremal) { + + // Create a point to having maximal x,y,z + Point extr(0.,0.,0.); + + // Loop over elems: find max point + auto itelem = meshPtr->elements_begin(); + auto endelem = meshPtr->elements_end(); + for( ; itelem!=endelem; ++itelem){ + // Get a reference to current elem + Elem& elem = **itelem; + // Loop over elem's nodes + size_t nNodes = elem.n_nodes(); + for(size_t iNode=0; iNode 0 + double xcoord = pointNow(0); + double ycoord = pointNow(1); + double zcoord = pointNow(2); + if(xcoord > extr(0) && ycoord > 1.0 && zcoord > 1.0){ + extr = pointNow; + } + } + } + + // Having found extr point subtract epsilon to ensure fully inside element + Point epsilon (5e-4,5e-4,5e-4); + extr -= epsilon; + + // Convert to Openmc position + openmc::Position pt(extr(0),extr(1),extr(2)); + + // Get OpenMC Mesh + openmc::Mesh* omcMeshPtr = openmc::model::meshes.back().get(); + + // Find bin of extr pt + int iBin = omcMeshPtr->get_bin(pt); + + // Upcast mesh ptr as unstructed + openmc::UnstructuredMesh* unstrMeshPtr = dynamic_cast(omcMeshPtr); + // Get the volume of the bin + double volume= unstrMeshPtr->volume(iBin); + + // Get the openmc tally result for this bin + ASSERT_FALSE(openmc::model::tallies.empty()); + openmc::Tally& tally = *(openmc::model::tallies.at(0)); + xt::xtensor & results = tally.results_; + int iScore = 0; + int nBatches = tally.n_realizations_; + double omc_val = results(iBin,iScore,1)/double(nBatches)/volume; + + // Now compare this result to that obtained from a function + double function_val = functionUOPtr->value(extr); + EXPECT_LT(fabs(function_val-omc_val),tol); + } From 2b9187c809194ae6ea8c21218e5318f7adb98dfa Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Tue, 15 Jun 2021 13:54:53 +0100 Subject: [PATCH 07/10] Resize containers in the event we refine mesh --- openmc/src/userobject/MoabUserObject.C | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/openmc/src/userobject/MoabUserObject.C b/openmc/src/userobject/MoabUserObject.C index 14e72fde..8b1cf53f 100644 --- a/openmc/src/userobject/MoabUserObject.C +++ b/openmc/src/userobject/MoabUserObject.C @@ -1056,7 +1056,8 @@ MoabUserObject::sortElemsByResults() for(const auto & elemSet : sortedElems){ elemCountCheck += elemSet.size(); } - if(elemCountCheck != mesh().n_elem()){ + + if(elemCountCheck != mesh().n_active_elem()){ mooseError("Disparity in number of sorted elements."); } @@ -1285,6 +1286,15 @@ MoabUserObject::resetContainers() // Update the serial solutions for(const auto& sol : serial_solutions){ System & sys = systems().get_system(sol.first); + + // Check if solution vector size has changed, e.g. due to mesh refinement + if(sys.n_dofs() != sol.second->size()){ + // clear + sol.second->init(0,false,SERIAL); + // resize + sol.second->init(sys.n_dofs(),false,SERIAL); + } + sys.solution->localize(*sol.second); } } From 0e53659c154b9dee906b70d51dfc0ec8156fb8fc Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Tue, 15 Jun 2021 14:32:01 +0100 Subject: [PATCH 08/10] Loop over active elements in mesh so it works for refined case --- openmc/src/userobject/MoabUserObject.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/src/userobject/MoabUserObject.C b/openmc/src/userobject/MoabUserObject.C index 8b1cf53f..0674e446 100644 --- a/openmc/src/userobject/MoabUserObject.C +++ b/openmc/src/userobject/MoabUserObject.C @@ -419,8 +419,8 @@ MoabUserObject::createElems(std::map& node_id_to moab::Range all_elems; // Iterate over elements in the mesh - auto itelem = mesh().elements_begin(); - auto endelem = mesh().elements_end(); + auto itelem = mesh().active_elements_begin(); + auto endelem = mesh().active_elements_end(); for( ; itelem!=endelem; ++itelem){ // Get a reference to current elem @@ -722,8 +722,8 @@ MoabUserObject::setSolution(unsigned int iSysNow, unsigned int iVarNow, std::ve bool procHasNonZeroResult=false; // When we set the solution, we only want to set dofs that belong to this process - auto itelem = mesh().local_elements_begin(); - auto endelem = mesh().local_elements_end(); + auto itelem = mesh().active_local_elements_begin(); + auto endelem = mesh().active_local_elements_end(); for( ; itelem!=endelem; ++itelem){ Elem& elem = **itelem; From 559d9ccee237e9bb4712c919ed99ee9f9241952f Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Tue, 15 Jun 2021 14:58:22 +0100 Subject: [PATCH 09/10] Add test fixture for case with refined mesh, having same base class fixture and setup as deformed test fixture. --- unit/inputs/aurora-refine.i | 147 ++++++++++++++++++++++++++++++++++++ unit/inputs/openmc-refine.i | 38 ++++++++++ unit/src/AuroraAppTest.C | 20 ++++- 3 files changed, 203 insertions(+), 2 deletions(-) create mode 100644 unit/inputs/aurora-refine.i create mode 100644 unit/inputs/openmc-refine.i diff --git a/unit/inputs/aurora-refine.i b/unit/inputs/aurora-refine.i new file mode 100644 index 00000000..53220506 --- /dev/null +++ b/unit/inputs/aurora-refine.i @@ -0,0 +1,147 @@ +[Mesh] + [cube] + # Length dimensions are cm + type = GeneratedMeshGenerator + dim = 3 + nx = 5 + ny = 5 + nz = 5 + xmax = 1 + ymax = 1 + zmax = 1 + elem_type=TET4 + [] + [cube_id] + type = SubdomainIDGenerator + input = cube + subdomain_id = 1 + [] +[] + +[Executioner] + type = Transient + solve_type = 'PJFNK' + petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart' + petsc_options_value = 'hypre boomeramg 101' + num_steps=2 + dt = 0.1 +[] + +[Variables] + [temp] + order = FIRST + family = LAGRANGE + initial_condition = 300 # Start at room temperature + [] +[] + +[AuxVariables] + [heating-local] + order = CONSTANT + family = MONOMIAL + [] +[] + +[Functions] + [heat-function] + type = ParsedFunction + value = (q1-q0)*(x-x0)+q0 + vars = 'q1 q0 x0 ' + vals = '1000 0 0' + [] +[] + +[Kernels] + [heat_conduction] + type = ADHeatConduction + variable = temp + [] + [heat_conduction_time_derivative] + type = ADHeatConductionTimeDerivative + variable = temp + [] + [heat-source] + type = ADMatHeatSource + material_property = volumetric_heat + variable = temp + [] +[] + +[BCs] + [temp-bc] + type = DirichletBC + variable = temp + boundary = left + value = 300 + [] +[] + +[Materials] + [copper] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity specific_heat density' + prop_values = '3.97 0.385 8.920' # W/cm*K, J/g-K, g/cm^3 + block = 1 + [] + [heating] + type = ADGenericFunctionMaterial + prop_names = 'volumetric_heat' + prop_values = 'heat-function' + block = '1' + constant_on = 'ELEMENT' + [] +[] + +[MultiApps] + [openmc] + type = FullSolveMultiApp + app_type = OpenMCApp + execute_on = "timestep_begin" + input_files = "openmc-refine.i" + positions = '0 0 0' + library_path = ../openmc/lib + [] +[] + +[Transfers] + [./to_openmc] + type = MoabMeshTransfer + direction = to_multiapp + multi_app = openmc + moabname = "moab" + apptype_from="AuroraApp" + [../] +[] + +[Adaptivity] + marker = error_frac + max_h_level = 1 + start_time = 0.1 + [Indicators] + [temperature_jump] + type = GradientJumpIndicator + variable = temp + scale_by_flux_faces = true + [] + [] + [Markers] + [error_frac] + type = ErrorFractionMarker + indicator = temperature_jump + refine = 1 + coarsen = 0 + [] + [] +[] + +[UserObjects] + [uo-heating-function] + type = FunctionUserObject + variable = heating-local + execute_on = "initial" + [] +[] + +[Outputs] + console=false +[] diff --git a/unit/inputs/openmc-refine.i b/unit/inputs/openmc-refine.i new file mode 100644 index 00000000..f0fdb613 --- /dev/null +++ b/unit/inputs/openmc-refine.i @@ -0,0 +1,38 @@ +[Mesh] + type = GeneratedMesh + dim = 1 + nx = 1 +[] + +[Problem] + type = OpenMCProblem +[] + +[Executioner] + type = OpenMCExecutioner + variables = 'heating-local' + err_variables = '' + score_names = 'heating-local' + tally_ids = '1' + neutron_source = 1.0e18 + launch_threads=false + n_threads=1 + no_scaling = true +[] + +[UserObjects] + [moab] + type = MoabUserObject + bin_varname = "temp" + material_names = 'copper' + material_openmc_names = 'copper' + length_scale = 1.0 + graveyard_scale_inner = 1.05 + n_bins = 1 + var_max = 500. + [] +[] + +[Outputs] + console = false +[] diff --git a/unit/src/AuroraAppTest.C b/unit/src/AuroraAppTest.C index 13a8d31e..24a574e6 100644 --- a/unit/src/AuroraAppTest.C +++ b/unit/src/AuroraAppTest.C @@ -79,9 +79,9 @@ protected: }; -class DeformedMeshTest : public AuroraAppRunTest{ +class CubeMeshTest : public AuroraAppRunTest{ protected: - DeformedMeshTest() : AuroraAppRunTest("aurora-thermal-expansion.i") { + CubeMeshTest(std::string inputfile) : AuroraAppRunTest(inputfile) { // Override source openmc input names openmcInputXMLFilesSrc.at(0) = "settings-deformed.xml"; @@ -117,6 +117,17 @@ protected: }; +class DeformedMeshTest : public CubeMeshTest{ +protected: + DeformedMeshTest() : CubeMeshTest("aurora-thermal-expansion.i") {}; +}; + +class RefinedMeshTest : public CubeMeshTest{ +protected: + RefinedMeshTest() : CubeMeshTest("aurora-refine.i") {}; +}; + + TEST_F(FullRunTest, UWUW) { @@ -193,3 +204,8 @@ TEST_F(DeformedMeshTest, CheckExtremal) EXPECT_LT(fabs(function_val-omc_val),tol); } + + +TEST_F(RefinedMeshTest, CheckBins) +{ +} From d983eaead99ae6e2d6de2e5bed69182a18dba753 Mon Sep 17 00:00:00 2001 From: helen-brooks Date: Tue, 15 Jun 2021 15:35:17 +0100 Subject: [PATCH 10/10] Completed integration test for uniformly refined mesh case --- unit/src/AuroraAppTest.C | 51 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/unit/src/AuroraAppTest.C b/unit/src/AuroraAppTest.C index 24a574e6..c42f1689 100644 --- a/unit/src/AuroraAppTest.C +++ b/unit/src/AuroraAppTest.C @@ -86,6 +86,8 @@ protected: // Override source openmc input names openmcInputXMLFilesSrc.at(0) = "settings-deformed.xml"; + hasDisplaced = false; + }; virtual void SetUp() override { @@ -104,6 +106,10 @@ protected: FEProblemBase& problem = app->getExecutioner()->feProblem(); if(problem.haveDisplaced()){ meshPtr = &(problem.getDisplacedProblem()->mesh().getMesh()); + hasDisplaced = true; + } + else{ + meshPtr = &(problem.mesh().getMesh()); } ASSERT_NE(meshPtr,nullptr); @@ -114,6 +120,7 @@ protected: MeshBase* meshPtr; FunctionUserObject* functionUOPtr; + bool hasDisplaced; }; @@ -150,6 +157,9 @@ TEST_F(FullRunTest, Legacy) TEST_F(DeformedMeshTest, CheckExtremal) { + // Check we indeed have a deformed mesh + ASSERT_TRUE(hasDisplaced); + // Create a point to having maximal x,y,z Point extr(0.,0.,0.); @@ -208,4 +218,45 @@ TEST_F(DeformedMeshTest, CheckExtremal) TEST_F(RefinedMeshTest, CheckBins) { + ASSERT_FALSE(hasDisplaced); + + // Total and active number of elements are different in refined case + ASSERT_NE(meshPtr->n_elem(),meshPtr->n_active_elem()); + + // Get the openmc tally results + ASSERT_FALSE(openmc::model::tallies.empty()); + openmc::Tally& tally = *(openmc::model::tallies.at(0)); + xt::xtensor & results = tally.results_; + int iScore = 0; + int nBatches = tally.n_realizations_; + + // Get OpenMC Unstructred Mesh + openmc::UnstructuredMesh* unstrMeshPtr + = dynamic_cast(openmc::model::meshes.back().get()); + ASSERT_NE(unstrMeshPtr,nullptr); + + // Check that everything has the correct number of bins + size_t nBins = meshPtr->n_active_elem(); + EXPECT_EQ(tally.n_filter_bins(),nBins); + EXPECT_EQ(results.shape()[0],nBins); + EXPECT_EQ(unstrMeshPtr->n_bins(),nBins); + + for(size_t iBin=0; iBinvolume(iBin); + openmc::Position pos = unstrMeshPtr->centroid(iBin); + + // Convert pos to a libmesh point + Point pt(pos[0],pos[1],pos[2]); + + // Get the result for this bin + double omc_val = results(iBin,iScore,1)/double(nBatches)/volume; + + // Now compare this result to that obtained from a function + double function_val = functionUOPtr->value(pt); + double reldiff = fabs(function_val-omc_val); + if(omc_val > 0.) reldiff /= omc_val; + EXPECT_LT(reldiff,tol); + } + }