From c1c0ffdad3975f94c336d53cc216a06a139f68ed Mon Sep 17 00:00:00 2001 From: Daniel Schwen Date: Fri, 3 Jan 2025 14:14:51 -0700 Subject: [PATCH 1/2] Add time transformation (#29639) --- .../include/userobjects/SolutionUserObject.h | 18 ++++-- .../src/postprocessors/ParsedPostprocessor.C | 2 +- .../src/userobjects/SolutionUserObject.C | 55 +++++++++++++----- .../gold/solution_aux_exodus_interp_out2.e | Bin 0 -> 75364 bytes test/tests/auxkernels/solution_aux/tests | 9 +++ 5 files changed, 62 insertions(+), 22 deletions(-) create mode 100644 test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e diff --git a/framework/include/userobjects/SolutionUserObject.h b/framework/include/userobjects/SolutionUserObject.h index caa608425e36..c1349f0a440a 100644 --- a/framework/include/userobjects/SolutionUserObject.h +++ b/framework/include/userobjects/SolutionUserObject.h @@ -11,6 +11,7 @@ // MOOSE includes #include "GeneralUserObject.h" +#include "FunctionParserUtils.h" // Forward declarations namespace libMesh @@ -22,12 +23,13 @@ class MeshFunction; template class NumericVector; } +class Function; /** * User object that reads an existing solution from an input file and * uses it in the current simulation. */ -class SolutionUserObject : public GeneralUserObject +class SolutionUserObject : public GeneralUserObject, public FunctionParserUtils { public: static InputParameters validParams(); @@ -321,15 +323,13 @@ class SolutionUserObject : public GeneralUserObject /** * Updates the times for interpolating ExodusII data - * @param time The new time */ - void updateExodusTimeInterpolation(Real time); + void updateExodusTimeInterpolation(); /** * Updates the time indices to interpolate between for ExodusII data - * @param time The new time */ - bool updateExodusBracketingTimeIndices(Real time); + bool updateExodusBracketingTimeIndices(); /** * A wrapper method for calling the various MeshFunctions used for reading the data @@ -453,7 +453,7 @@ class SolutionUserObject : public GeneralUserObject /// Pointer to second serial solution, used for interpolation std::unique_ptr> _serialized_solution2; - /// Interpolation time + /// Time in the current simulation at which the solution interpolation was last updated Real _interpolation_time; /// Interpolation weight factor @@ -507,6 +507,12 @@ class SolutionUserObject : public GeneralUserObject /// Map from block names to block IDs. Read from the ExodusII file std::map _block_id_to_name; + /// function parser object for the resudual and on-diagonal Jacobian + SymFunctionPtr _time_transformation; + + /// Time at which to sample the solution + Real _solution_time; + private: static Threads::spin_mutex _solution_user_object_mutex; }; diff --git a/framework/src/postprocessors/ParsedPostprocessor.C b/framework/src/postprocessors/ParsedPostprocessor.C index 21b795e22d56..5a221eddd749 100644 --- a/framework/src/postprocessors/ParsedPostprocessor.C +++ b/framework/src/postprocessors/ParsedPostprocessor.C @@ -31,7 +31,7 @@ ParsedPostprocessor::validParams() {}, "Vector of values for the constants in constant_names (can be an FParser expression)"); params.addParam( - "use_t", false, "Make time (t) variables available in the function expression."); + "use_t", false, "Make time (t) variable available in the function expression."); params.addClassDescription("Computes a parsed expression with post-processors"); return params; diff --git a/framework/src/userobjects/SolutionUserObject.C b/framework/src/userobjects/SolutionUserObject.C index 2d8d9d47f203..c4dafb809664 100644 --- a/framework/src/userobjects/SolutionUserObject.C +++ b/framework/src/userobjects/SolutionUserObject.C @@ -16,6 +16,7 @@ #include "MooseUtils.h" #include "MooseVariableFE.h" #include "RotationMatrix.h" +#include "Function.h" // libMesh includes #include "libmesh/equation_systems.h" @@ -36,6 +37,7 @@ SolutionUserObject::validParams() { // Get the input parameters from the parent class InputParameters params = GeneralUserObject::validParams(); + params += FunctionParserUtils::validParams(); // Add required parameters params.addRequiredParam( @@ -86,6 +88,13 @@ SolutionUserObject::validParams() 0.0, "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector."); + params.addCustomTypeParam( + "time_transformation", + "t", + "FunctionExpression", + "Expression to transform from current simulation time to time at " + "which to sample the solution."); + // following lines build the default_transformation_order MultiMooseEnum default_transformation_order( "rotation0 translation scale rotation1 scale_multiplier", "translation scale"); @@ -110,6 +119,7 @@ Threads::spin_mutex SolutionUserObject::_solution_user_object_mutex; SolutionUserObject::SolutionUserObject(const InputParameters & parameters) : GeneralUserObject(parameters), + FunctionParserUtils(parameters), _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")), _mesh_file(getParam("mesh")), _es_file(getParam("es")), @@ -137,12 +147,12 @@ SolutionUserObject::SolutionUserObject(const InputParameters & parameters) _initialized(false) { // form rotation matrices with the specified angles - Real halfPi = std::acos(0.0); + Real halfPi = libMesh::pi / 2.0; Real a; Real b; - a = std::cos(halfPi * -_rotation0_angle / 90); - b = std::sin(halfPi * -_rotation0_angle / 90); + a = std::cos(halfPi * -_rotation0_angle / 90.0); + b = std::sin(halfPi * -_rotation0_angle / 90.0); // the following is an anticlockwise rotation about z RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1); // form the rotation matrix that will take rotation0_vector to the z axis @@ -151,8 +161,8 @@ SolutionUserObject::SolutionUserObject(const InputParameters & parameters) // back _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z); - a = std::cos(halfPi * -_rotation1_angle / 90); - b = std::sin(halfPi * -_rotation1_angle / 90); + a = std::cos(halfPi * -_rotation1_angle / 90.0); + b = std::sin(halfPi * -_rotation1_angle / 90.0); // the following is an anticlockwise rotation about z RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1); // form the rotation matrix that will take rotation1_vector to the z axis @@ -164,6 +174,18 @@ SolutionUserObject::SolutionUserObject(const InputParameters & parameters) if (isParamValid("timestep") && getParam("timestep") == "-1") mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply " "remove this parameter altogether for interpolation"); + + // setup parsed expression for the time transformation + _time_transformation = std::make_shared(); + setParserFeatureFlags(_time_transformation); + + // parse function + const auto & expression = getParam("time_transformation"); + if (_time_transformation->Parse(expression, "t") >= 0) + mooseError("Invalid parsed function\n", expression, "\n", _time_transformation->ErrorMsg()); + + // only parameter is time + _func_params.resize(1); } SolutionUserObject::~SolutionUserObject() {} @@ -344,7 +366,7 @@ SolutionUserObject::readExodusII() _es2->init(); // Update the times for interpolation (initially start at 0) - updateExodusBracketingTimeIndices(0.0); + updateExodusBracketingTimeIndices(); // Copy the solutions from the first system for (const auto & var_name : _nodal_variables) @@ -452,7 +474,7 @@ SolutionUserObject::timestepSetup() { // Update time interpolation for ExodusII solution if (_file_type == 1 && _interpolate_times) - updateExodusTimeInterpolation(_t); + updateExodusTimeInterpolation(); } void @@ -572,11 +594,11 @@ SolutionUserObject::getSolutionFileType() const } void -SolutionUserObject::updateExodusTimeInterpolation(Real time) +SolutionUserObject::updateExodusTimeInterpolation() { - if (time != _interpolation_time) + if (_t != _interpolation_time) { - if (updateExodusBracketingTimeIndices(time)) + if (updateExodusBracketingTimeIndices()) { for (const auto & var_name : _nodal_variables) @@ -603,12 +625,12 @@ SolutionUserObject::updateExodusTimeInterpolation(Real time) _es2->update(); _system2->solution->localize(*_serialized_solution2); } - _interpolation_time = time; + _interpolation_time = _t; } } bool -SolutionUserObject::updateExodusBracketingTimeIndices(Real time) +SolutionUserObject::updateExodusBracketingTimeIndices() { if (_file_type != 1) mooseError( @@ -619,7 +641,10 @@ SolutionUserObject::updateExodusBracketingTimeIndices(Real time) int num_exo_times = _exodus_times->size(); - if (time < (*_exodus_times)[0]) + _func_params[0] = _t; + Real solution_time = evaluate(_time_transformation); + + if (solution_time < (*_exodus_times)[0]) { _exodus_index1 = 0; _exodus_index2 = 0; @@ -629,12 +654,12 @@ SolutionUserObject::updateExodusBracketingTimeIndices(Real time) { for (int i = 0; i < num_exo_times - 1; ++i) { - if (time <= (*_exodus_times)[i + 1]) + if (solution_time <= (*_exodus_times)[i + 1]) { _exodus_index1 = i; _exodus_index2 = i + 1; _interpolation_factor = - (time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); + (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); break; } else if (i == num_exo_times - 2) diff --git a/test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e b/test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e new file mode 100644 index 0000000000000000000000000000000000000000..3967749bb9555b7570f0c930de0d3f65538bd6dd GIT binary patch literal 75364 zcmeHwTdy3)b?(;Hq$Jz2E!*-XPDgSq$r7)kE|k5!hoZJ*+Po+vCCiqb@$8wNy<41{ z=0a{NHlTlSUeBKhke?vr$&Y$45Cg{#gkc1EI0%fu2@;&(;Cx?IeLdYZ-7{0Y$37G} zO|IS3U8}xYwN`ai^{PvseeLyoH#RmNa`ZWz4|U>Rz1Qf*8yolFa^vtg>BjZZI3DoI z9-rz>y7gAl<<)(l!@TbGTk*(-{}QJ|ycc(3GSxs2+3WgF=Nhgo_@@=GK8nXKF6A-z zM@b9m=I(#SVI)1moG>^w{Qa!MsSn~I(lh~waQ_kBw>Tr-eAK&l(%pMl@q31C<14tI zxktAAz3Yu3?@?Bh&!pGx5c0*`ok9T^8vB%XM7OzdAOZ)peq}Y31j#dx{dv$ zOZ&j4bqKdR-ox!QeuUqT0Q|*MFmhu&p&5jYpBoR1lW~U$;$VKC81*}oF)~qaO!n*X ze!n#t)e%1)4(k2M_~cQ?`!IOF271VppBsC92SS-l!NQrxgCOPs$3OiMe-6hxQ$cV| zUg}B*Wh3D@^cw$we=zJ1XkYPC=lp%YEyqUg`_7Tu##Q|9(m`G)-=XvT z4gPF;t;n8e*EFU_u>HQQp4}m|M9(k8duLa%*n|+j;eFIOO@Ojk%+3))o zdvU0Pv*o>M?#$tRg!!<%A1va1^t<4ojc5v%EtVnKQ~7wP{(92syxr(bY~djf)8~QT z`uW&5aZUYav1F?yu?WrMiv|n2&0@(8fmmNA)cv)UT}QfqiWKK7H_{2h(>K-y;7Z!#^?=Cy?Dy-X6od5yzd~6a9IhHm-UK{ zEbjlr%|i}%5$|uvYxSGTOFiizZp&UgZRB;l!u_-uanADo0@shQK(R5b+eU(V{*2?7 zEf$n_M&a@FE&TcP^+pG!)Z)>&ILJUF=&$d$saQ^-=d-3gCI8zp%l$P%EdhW6Hu?;Re`;mRVKUcWYReFoAeZEEf zi^5lhw{^=Red5oq1>5?hpa0 zwxQ*J(|q}I=-N-$OUFmP?eUbaW%V7JU+K%mlYVCOa!|f?KBU^*dF6UmIg!3g?3NQ^ zeX-WSS_5khtTnLKz*+-$iw5Y!qyLWnM*2JHzo*}w{_;a0`kCqPrEiu#PWngb^OWE4 zbNKyv&=)}T8-5Y=CD4~akApbg_Eiwai@px}3(z+}e+h~}e+7C1^d#t;puYxv3-lD| zY0xvEZ-br%JqJ1hItn@lIu4@SlrGTcL8m}3fL;W>1bP{C8uSY2J0Ol%Y=XWA+5*); zXF#ul&VpV8ZG&D1y#e|@=p5(=pdW(HgD!wBf-ZsH1pN)@GUy8ED(EfH+n^tT-U0m> z^b^p#pr3+%2Krl29rSZh1GEEbf?6QD@Y|qW&>koOy$8Am>VUeS9;go*fZhiULGy$5tJ=sr-9%^4=UlZWL&_Gh_Zn9705hRTD=h{}b^ zips~Mh|4sItFpp!MqHH{h7p%#hdc};F3S;l7)D%{De^FkxGZ1fVHnFK%NltY#;;Ia?oDx@MkzvGr5=0(`5%-%Q@-U3JEUV;U7;#x{$-^+F{S=5i3}f0X&*X_f z#C-WcdmKa_hB0kE8{}aa)8=zRo(Dh= zf~Z&O-hXx&mbEvZtCv80o?dWaS$i|=H0Tu&c^H3?uG25P29zo66n{Bkm7C~5y-bU z!-)G1h&;Z%8AjY6gUI9Cn_*1*Cm{0p_GX&IeHZjokZ*5>5%*^x^7!^<7;*m=M4lQb zYj1`z?K+4&zP*_yaeoeKfP8y1jJP|XpM!jRGmN-R5P5uiGmL4sK;-f59phdbvxUL>}MXOp~~M&;aDyn_$G11rB<@Yn2O!_x z3?uH}fym?An_o&y~L9R(c&9S5BNodi7(It6+G^djgb(959HpjSZO0eu&=3HlyL z&+S3;PkkMP4k|xwbrACNFMG?*vajq)d(xJ)qijgq(QdSv>_t1#M$BIfYJ>Rf^7&3c zd`?wQbn%-%F!;g-=HNv#j zW25HYjn=U@hyC4QqdPj*?e|CV5lsd?!dae2PM&`8`P0Y7@n~#=>->47KNv?xlE|X^p2UCqvK(tH){8XT}*`S_cq7JPDYm|ob?L8 z-_>^-qj*!Il*~(wjvR@H!~Rg`BYjM|9b%3#Zgi5Jaevs{JJO1;18^g0#!+)Hp(@|D zGj^kMm=$FenA%bH`JJyxBV4dhw-3)bB+nkDokss3Tw2N3@$W! zjop|Xm%<~$1KS%9`dYMp4qEY}a#QWq>!Mr+fW+ zD{eO?oiRIR3l0J%CbQpJN+A2tFhCcpqVN-!X<2U09o(FyHNY}xc>f) zIN9ApqvwvzD6GnXQ6{5tzl((?{pPh4MsT!ps16+<=GCH2nF8roxDGkv8WmXGAPQ!B zMmc=z>ZJvy%>+kL36LPe0nL2ua@24@&R>@}1?SD^jWcguxsu7${GThY9OfriA@tod zll>pY!(QAOt-2zsAWcz}D0=tkv0iU=($TJSt5{RPTBV#ZRdA!#GQGZQMp09ZK`cs4 zz`0dy_@lR=iVd$E064OpUz;V$oH@d)uPB=PtZYVC+&V6(tH`Y0J<38gx6!bm6Fqxu zM#sG<^Dl{D>eg8(d#8#~OncH8-Ymqd0@&P(o7ZeJ1z0SL7>$QXlka~xsu-n}U};5h zBTB$<9mYM=;VjH?B?)8Gi;e4eXP7-Uqo>ZFySRPka&-0Fh3zX>x8JO9Uwn=Jo#KNC z?ID^`EQU#n1g&BeE{Ch{#$&8GtdAz`cCw!lw92D5YBsQbvZQENvGORkB_t!)^k)?- z>*aQelKHjCck2L87o$!_F;?)72AxJTZq>0gZ!ZTJUJeOpmx1=t^93eXld7$*lI5}z zW8~b7!d!c7q>dh2T^qIW$>=$}?Vh^Zv+wZ#spo(OPkOO!ZRW+pAI^>TM5MZ{fW!-x5@*u~*)=DSxMcJ~pG* z&t1N9)mK&NHFd})01NH;bM~l`_8fBHwcF#`fw zY_6g`o{N$RECxiC%K#o#vfLOB8)fAhfZ=V*%KL}#j1Kxe;|#ai%>J#?(8I|P7xj}2 zZj~&rygS#vp3S58yNZO^1>e8{OEcQ9PdUon>*l9vb$BYcASeYLLTzSu7^$q=> z+2`1g zY$c7|`ltb)b28f!ROx-`wagCMGKtJlCCe*sQ{sJ$%3@y<_(Jw&i&Vx+5*om~aSsHT zbO)WJi9yh!qe??xl35E6uLxSNaommAC0ks47;6AuQd}$_dtafW73jg16hR|mm2$Xp zC&tQs4yI19CYb#j-q1%8I>3cqDyPG+=`@T~Qlzc4n-R)D(n3$guUCq%??xwv9aG0E zXgpj=Q@@~@=c8w5_+%2cjeh1v);wsUCgVP~=)#UU?rk~$9vS{D&+Ru(aT;gv(9*^CDU5gY2& z;d6$I)9%Wbzz=Q4`XT1huPTr*2RKNY4Q%dL3CpO#mBU8sdZUNwOS4CKyUc>kX;;m9 zp?4%4(2k}wCRu-dM$iEDxh(U2elo_^ht_)od*`xV?a{*sm)XD(n@=_{DXea07Qot_ zMt7&xsAG!H2$t3~s?svZ0Q(N$19rxJ6~+oWw9^vMFRo8vrb{-PBPR{^(^TWJ8+dJ&6a?Axh+G z0SoFdW0aukE7^t4=)kb=7Na^8iRmiizFO6=DmIEb^@75ZwQ^_^(80U>j485!FCSmN zYT%PnT=S=hR)@k)X-{#q%8ruO(CO zV)^O;R_?F{4+xRv|J9!CTEJ#N{IzKHP*sg*}|gG$W?=3~zrYUQ&W zY82C`)v+K84z&V39uN>&o>>ZUAPHyR_bNze9V=x$rttdEz2g=1hE2H%x3U8)t8qe} z5|(rfE!ZrVRS+i2K|NR?HTI2*t8kKJxP%4}4M5r_wvU|thC3w9U`pBQ6h0}TkDhsj zuU{Y4-)p4qFLV^Kj<2GXQ=b49W)if&lQep0@-dHfFCJEtkbN8U2yOs>_Q53q0QN=np(iAmZhofFR>b5cQoT@Cc^72intsJaBy!5?d zGrG9_qpO!LDocyb$TCJ+$NWgmG@IF9;pGq?<*<(8;y!G6^##Cr&Z-JVF%Mp;Z>b7k z+-TjjQ>L*lZwz;5cGl$q0nIy5&J}h*Ko|x?9$8;|1*6#aQ1khfYo5*rcoX&1l(Z zjHMgRcr@y}1(4yDLl$tR0brIxhcd0Z#=34AWAl_4mF~p5^lyh}nH9D1j5IJ{G3&># z28dRJTrY^oOi)LL(Qb z6S?T9VwClS#{4b}Rdi@)NS4h0ssgANks>3rp)V_Ph4%@vU_zrECMb+uaFc+mUzK_H zO=x7_;m^5{fWNs(p{J4+b5VL!YaBNuWvPH5yVORH++smukj>91u?mHiVd9a+=xcrjv}|kL7j;=7XWUM7EcAJx7r#C3shMk|uH9 zQ4lWSEFcPF+wM(O$;vCA`a3ccr^_}m3eumg4~-ADh0T$S3w9X8YmaYO(;eGrbIXWf z+~&TB%NoWxbyJiXPv^2gs%D+b)Hv?bvX}i4Zsk&Uu7$;dXbcj}R`rs|KVu1HFTns; zSGvkd$KT;uW(8|J#Xgk14hQ6BG<@i=hz9jT{xLBUhPEkLzyhyicztL9Cp#_7q2|0a zU&ez1&d68+(CnFpO%NyD9XL*X7k+rF*XRuP8ueixTYqHz*RpxIggfzg4j#&3sFmj~ zDX-M(O_C0FpXy+C2TCK_d+ckM=|O}S>?RgD4EqUI7pUGson@y$Q>WosZgx1AiVq*= z(6M34KP(j_Ok2IQnX|24+6Y3*hK%LcHtR@{u<`SDNuN^j0J zCwOP3(Pm|{zeBRibp8ahojD9VD_DAG6Z8Cc;!%HsiC0ILct1l|q=8L@SNSZf2aJL} z#E_#XA`zgJ|{$>N-!cd7;vdn(W(VLfZ>f)29PtoG14%7 zkPTpHr@h(l6z)0dV8d%0`fV_WB_~$qq7vZz9lS%QB{(0fUsFrSaHBLw^EQI2VC7c7 z-EEXC8iknQ63*Rtlz`bSM+?hx?k4)!`~B9WOoIem{vCzLTn$$^)pBZ-bgDbnIAa?C}d3q*R_B#aOkDdpJ9c#Nl=4VdK0l_DpoG`p{04p z#0mSFP$rUhAAW3pEx>FFl!9b`gd61s#}!QJE`$y**cg*CL!p?HP}JQH?>SIF&BHcp zJM%STc$VGEBkEW#vdF5>{0;~Flk6+P%4SA3dm~(j3pXkWeP5Po-3NkTY?;VdC%o7Q zFfpI{n#uS90cSroQ9+09s}&fL^0P~Et4lb*w4Ldyg;Z2*`FQT)xvS^SoX@?MmhSz? zw^4o_E=GrNfUU{!38)my=(9Q?SH+>rnKSW$vTOl8`KZIHS$2C)u(RpDedw|r`T{N; z&H@@`ve+1#cr|m`<s<`m8ja2V*^IU?zJ2cUrHdD~ zFD|qDJ9j?WF|Pnc{ozW*#;xX5Uma^+;L5+AV3Cik9J2euhTWS51v&-v>c{ZP!I#im z5El-}fx=<`M!v0^8V9w_vWr?x(ZptH8l&A1-ITfJUDr2f=(Gv!_aCDo-UH(p6=j!08jBUzH_H zjV&a1t_gS7kwKd9B{y1yFWQg^06GA@=DtLOXvIA!;)jKM5j6o?QoVl zSS(_=Q3Pnc=Q;l?Au9mpyK#-tv>#RtaMa{(v`HJgqw)pKW_0>DTtYg6rVB@BrcqY_ zEN@S2mQNgg>GVr4p614nOKiu9X{-ED2WGg<`q50j!z%|{J(s>qRRt`vM^y>w`@+^p znF%LT8&@fZqtuz^{3BU70){?%nXYgi%3h(n*-k$`&VC=ZtKOeB=I*ToE*aN~;j@uk zM~@NLNq6>h2EYujz6>BnmrFcg;RV|O77Q%A{A2IF4u(VKSAv;Pmt}5zl5i>Kj2gp@ zVk9i-dRZN!_qBYj?R~CqM!*X7p?yEeH+uD4tr%n%E(}}gqVg!3`o7AXxh!0IXh%mv9kPgQ(b17)71d#-5I`@@UF!PV?FZ3ep|AZw654&W^1k+iNLc8}4t@4mHte1lu4>?VI5;OcEe`&nY6ZK@?PUsL z6&=o<-)FwEJ+bRZmz&)eMRRA*g-@KAb&cnI?XYs74|WzX*CbTY;ap#P6~MW!@ha_d z>N=c`Fl5vis6>{v0>G&nt7N$VFyqM%M5Khyt&3G`F}|YY)!X4~d{V%fN1BOMcpFUt zi>EYJc~a(jvIAAP!Ca+)+?t*@vPYG2=m0T29&2wgXw40DVk}iR>oX4kn3kTbfCn*3 zzF8J(l{Tfg`|l7dWf|`iH}I8{-0Q?t!l0#Q8;m#C*S;A&@4d~Wqxk)a`Iv-LTkcdvhjS-gR;g{Q^S@xqlU>CqTzQJ7lUrJ0-ig=2siMP`y2e>H=6bTLJbI?Z ztM9eiFBF!;&{ps6`@sEMqj-49Oi#W|KXh*f7ck!1Tg}jdb?Pqu!2MJ&c06EujZfV} znoM_xC2u_zkbnL*qbu%PH&exBR5tUIu`z)fP@1Yj7y8qDVY=}86 zSjmR(LrBJDbtmgE9tM~3sYCX(3qDqR06;7TEBN^BfdG@HGCDn10W@<1rV#mG;k(Ch zEmwa~$Jb4`<=SY{#)ls>m{mGU+E(g#9c$O}14z|gE2e<9${tg8`#Zgd~Eq? zwYFU*;7Z5!(vKwAZ#H4Co6JN%SlNWTJs^z%=ULojnYm3jjWijZ{sG|e2><4RWEIMo zlX?5CR;j)SXy%Cw`j8VVTTykVUZHTnyRS?9wmwh0svc+7yTt8AHz}Kbh#LF+nKv$lxXgZ%~8)?zOu~P8LMMPW2IR-BfE;vVeV4Y&FDfO-)m9Xdllb$ zn`>@LD5IF&h|RcWrj=o_xw^f4o{IsmPlrLOWVtaMHcF<|R_o5q-P*H?gw1{zUmwOS z@1k`+Rmx$U;G4))X33UVGktG*u+DgIq-lgbu10#Y{rS8^z5c%qsmu(-oOy zcX-KICp^oyea~U*maROtK6d9{dDp*o*E)+&*S~S$AGx#Q>n#3%@c2A_8@jQvB|e?S z=clK1f9r6Sj`;lag`@j^eBs=3c)}6iB78r-aDM0Tgd@I1_~Q5Z6n~LCAGvguU&WJL zKECi2F1Z!Y=U4deT{?>A^Os$V-_M`q(0$3T`--RZeSF1J_#*if{?9HSil=mhUv@2i z=|}vML-!@W?kk?s_wf}^;fv%`_`kS(D4x=B_%?L?cddc72JU1HZ2iIJZtIpit3J`S z>Ia?u`>ID&Zz;ac;{RyF9BKc-{iV)A@hfn5q2Y^V1iO?)&kD^9P3~ z9P!=q?h8lv{rJN9Hy2Mh;#-6-exFbAm9D~ldHnn;p5*fJg{N@It$04a!vE;fOZWZ!DxT!>@fAc zy-XGNO*V(_X`tFtm+ERR-#s9Gj7oX0K&xYcCU3@x=Pw6W?o&Vd$Oyw7!pT2N( z-;eKbY8Gx|OE}_FK6GC=y6?poPFl|jNAe0s_l4u%_pcS-=lA9D^Xti%_UDA}@Mh#z zdXm%eD_^zyt-=&<%jYk-PWhF-@E7Du@s$2ddV%t#@;~VE;o3!|BmA=KR1V3n{4U6s z;-%@^Ts~wMFF!7s^=YkvwFd5B4XEBo+cCXos9yRnoB1uB{rh?^ zkl#`9br%1x9gg^Pc6?ItmiTlQpVC))?p(X)LlZTXUwnT0!qI&%zHm}~2uJBDAG$9b zci)dMoZ5pMGioD02V3$R>it9Yzs}zM+Pzjg>3h23>n#4?czk-lbbQkBmiTlQ zpV|eNp3<#7;Bb|W=`d)f|{MsRl zBYpVkEyACP|A@^|?XwP7?VHk(UZp>$CqLbpYZrf!Jia|<CNXayH+|bec2{p z0qOX*-aGaFC)`{1N6M#hRLAGoHnzm4v&MauzTRteR(HGM?Q}`nJ6#gZb55<#w!Y{iPzvLIc=owCiBD(8m!|L1tLd!wk$WwI_%yG%=I;CQ&BPau>cIy;G=M4J0{ry-OYwDn ztd_@7It%dAbbUU>Um(v+eicu0J#56|#q;?U{-qCX{)FrEmtBkB&)V#Z&q|zVH;jNIr#s&17Svrg%!n;oH#l-?awT z8n}}+P*cBN??-ODRqws&ydBNM(R+xl-8f!4Z_UL|#}|eBd#lBi?~m@<@k__)d@qkb zl|$o(I?Lba$u*O{mtM_}FC6K^%b%|gA3ybbD!$IDU;XqJ;HUEWe3COXzltZhB#-cg zr*O%wcs{?vzhW(ymLs3P>{{t4ec?+EcYj9Dil_8_eBmj4k$ejO+J`11)ArEiOZZCH wzg9YuU*VEN_q}|j=_q|4U-1;aNIr#s!^y39lF#AW(DmQ72G$z5lQr=F0AdIpKL7v# literal 0 HcmV?d00001 diff --git a/test/tests/auxkernels/solution_aux/tests b/test/tests/auxkernels/solution_aux/tests index 3e64f0ec14d1..27a873f1b54e 100644 --- a/test/tests/auxkernels/solution_aux/tests +++ b/test/tests/auxkernels/solution_aux/tests @@ -110,6 +110,15 @@ "temporal interpolation." [] + [exodus_time_transformation] + type = 'Exodiff' + input = 'solution_aux_exodus_interp.i' + exodiff = 'solution_aux_exodus_interp_out2.e' + cli_args = 'UserObjects/soln/time_transformation=t/2 Outputs/file_base=solution_aux_exodus_interp_out2' + requirement = "The SolutionAux object shall be capable of setting an auxiliary variable with " + "temporal interpolation and a time transformation function." + [] + [exodus_interp_restart] requirement = "The system shall be capable of initializing an auxiliary variable from an " "existing solution that" From 4f4d28a1c01f2ea0f2c636b0a48808b27506a1da Mon Sep 17 00:00:00 2001 From: Daniel Schwen Date: Tue, 7 Jan 2025 09:40:02 -0700 Subject: [PATCH 2/2] Refactor SolutionUserObject (#29639) --- framework/include/auxkernels/SolutionAux.h | 4 +- .../include/auxkernels/SolutionScalarAux.h | 4 +- .../Axisymmetric2D3DSolutionFunction.h | 4 +- .../include/functions/SolutionFunction.h | 4 +- framework/include/ics/ScalarSolutionIC.h | 4 +- framework/include/ics/SolutionIC.h | 4 +- .../include/userobjects/SolutionUserObject.h | 494 +------ .../userobjects/SolutionUserObjectBase.h | 514 +++++++ framework/src/auxkernels/SolutionAux.C | 4 +- framework/src/auxkernels/SolutionScalarAux.C | 4 +- .../Axisymmetric2D3DSolutionFunction.C | 4 +- framework/src/functions/SolutionFunction.C | 4 +- framework/src/ics/ScalarSolutionIC.C | 4 +- framework/src/ics/SolutionIC.C | 4 +- .../src/userobjects/SolutionUserObject.C | 1220 +--------------- .../src/userobjects/SolutionUserObjectBase.C | 1265 +++++++++++++++++ .../userobjects/AdjointSolutionUserObject.h | 10 +- .../userobjects/AdjointSolutionUserObject.C | 33 +- .../src/userobjects/SolutionRasterizer.C | 2 + .../postprocessors/TestDiscontinuousValuePP.h | 4 +- .../postprocessors/TestDiscontinuousValuePP.C | 4 +- 21 files changed, 1839 insertions(+), 1755 deletions(-) create mode 100644 framework/include/userobjects/SolutionUserObjectBase.h create mode 100644 framework/src/userobjects/SolutionUserObjectBase.C diff --git a/framework/include/auxkernels/SolutionAux.h b/framework/include/auxkernels/SolutionAux.h index 1ab91effa858..147e6ed13bc5 100644 --- a/framework/include/auxkernels/SolutionAux.h +++ b/framework/include/auxkernels/SolutionAux.h @@ -11,7 +11,7 @@ #include "AuxKernel.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * AuxKernel for reading a solution from file. @@ -41,7 +41,7 @@ class SolutionAux : public AuxKernel virtual Real computeValue() override; /// Reference to the SolutionUserObject storing the solution - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name of interest std::string _var_name; diff --git a/framework/include/auxkernels/SolutionScalarAux.h b/framework/include/auxkernels/SolutionScalarAux.h index 91b66380c3ed..80e4b33ec975 100644 --- a/framework/include/auxkernels/SolutionScalarAux.h +++ b/framework/include/auxkernels/SolutionScalarAux.h @@ -11,7 +11,7 @@ #include "AuxScalarKernel.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * AuxScalarKernel for reading a solution from file. @@ -32,7 +32,7 @@ class SolutionScalarAux : public AuxScalarKernel virtual Real computeValue() override; /// Reference to the SolutionUserObject storing the solution - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name of interest std::string _var_name; diff --git a/framework/include/functions/Axisymmetric2D3DSolutionFunction.h b/framework/include/functions/Axisymmetric2D3DSolutionFunction.h index 1cec3d6984ab..3dece4349f07 100644 --- a/framework/include/functions/Axisymmetric2D3DSolutionFunction.h +++ b/framework/include/functions/Axisymmetric2D3DSolutionFunction.h @@ -11,7 +11,7 @@ #include "Function.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Function for reading a 2D axisymmetric solution from file and mapping it to a @@ -44,7 +44,7 @@ class Axisymmetric2D3DSolutionFunction : public Function protected: /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; /// Factor to scale the solution by (default = 1) const Real _scale_factor; diff --git a/framework/include/functions/SolutionFunction.h b/framework/include/functions/SolutionFunction.h index 05ab6802777e..e648f25567b4 100644 --- a/framework/include/functions/SolutionFunction.h +++ b/framework/include/functions/SolutionFunction.h @@ -11,7 +11,7 @@ #include "Function.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** Function for reading a solution from file * Creates a function that extracts values from a solution read from a file, @@ -55,7 +55,7 @@ class SolutionFunction : public Function protected: /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; /// The local SolutionUserObject index for the variable extracted from the file unsigned int _solution_object_var_index; diff --git a/framework/include/ics/ScalarSolutionIC.h b/framework/include/ics/ScalarSolutionIC.h index 927317959c1a..1f79ddb49485 100644 --- a/framework/include/ics/ScalarSolutionIC.h +++ b/framework/include/ics/ScalarSolutionIC.h @@ -11,7 +11,7 @@ #include "ScalarInitialCondition.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Class for reading an initial condition from a solution user object @@ -25,7 +25,7 @@ class ScalarSolutionIC : public ScalarInitialCondition protected: /// SolutionUserObject containing the solution of interest - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name extracted from the SolutionUserObject const VariableName & _solution_object_var_name; diff --git a/framework/include/ics/SolutionIC.h b/framework/include/ics/SolutionIC.h index 7bd64b9ae4c7..bf5f7372834c 100644 --- a/framework/include/ics/SolutionIC.h +++ b/framework/include/ics/SolutionIC.h @@ -11,7 +11,7 @@ #include "InitialCondition.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Class for reading an initial condition from a solution user object @@ -27,7 +27,7 @@ class SolutionIC : public InitialCondition protected: /// SolutionUserObject containing the solution of interest - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name extracted from the SolutionUserObject const VariableName & _solution_object_var_name; diff --git a/framework/include/userobjects/SolutionUserObject.h b/framework/include/userobjects/SolutionUserObject.h index c1349f0a440a..799e68dc312e 100644 --- a/framework/include/userobjects/SolutionUserObject.h +++ b/framework/include/userobjects/SolutionUserObject.h @@ -10,509 +10,25 @@ #pragma once // MOOSE includes -#include "GeneralUserObject.h" +#include "SolutionUserObjectBase.h" #include "FunctionParserUtils.h" -// Forward declarations -namespace libMesh -{ -class ExodusII_IO; -class EquationSystems; -class System; -class MeshFunction; -template -class NumericVector; -} -class Function; - /** * User object that reads an existing solution from an input file and * uses it in the current simulation. */ -class SolutionUserObject : public GeneralUserObject, public FunctionParserUtils +class SolutionUserObject : public SolutionUserObjectBase, public FunctionParserUtils { public: static InputParameters validParams(); SolutionUserObject(const InputParameters & parameters); - virtual ~SolutionUserObject(); // empty dtor required for unique_ptr with forward declarations - - /** - * When reading ExodusII files, this will update the interpolation times - */ - virtual void timestepSetup() override; - - /** - * Returns the local index for a given variable name - * @param var_name The name of the variable for which the index is located - * @return The local index of the variable - */ - unsigned int getLocalVarIndex(const std::string & var_name) const; - - /** - * Returns a value at a specific location and variable checking for multiple values and weighting - * these values to - * obtain a single unique value (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValueWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type = weightingType(), - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValue(Real t, - const Point & p, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable for cases where the solution is - * multivalued at element faces - * Use pointValue for continuous shape functions or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map - discontinuousPointValue(Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable for cases where the solution is - * multivalued at element faces - * Use pointValue for continuous shape functions or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map - discontinuousPointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable checking for multiple values and - * weighting these values to - * obtain a single unique value (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradientWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type = weightingType(), - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradient(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradient(Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable for cases where the gradient is - * multivalued (e.g. at element faces) - * Use pointValueGradient for continuous gradients or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map discontinuousPointValueGradient( - Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable for cases where the gradient is - * multivalued (e.g. at element faces) - * Use pointValueGradient for continuous gradients or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map discontinuousPointValueGradient( - Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Return a value directly from a Node - * @param node A pointer to the node at which a value is desired - * @param var_name The variable from which to extract a value - * @return The desired value for the given node and variable name - */ - Real directValue(const Node * node, const std::string & var_name) const; - - /** - * Return a value from the centroid of an element - * @param elem A pointer to the element at which a value is desired - * @param var_name The variable from which to extract a value - * @return The desired value for the given element and variable name - */ - Real directValue(const Elem * elem, const std::string & var_name) const; - - /** - * Returns a value of a global variable - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param var_name The variable from which to extract a value - * @return The desired value for the given variable - */ - Real scalarValue(Real t, const std::string & var_name) const; - - // Required pure virtual function (not used) - virtual void initialize() override; - - // Required pure virtual function (not used) - virtual void finalize() override; - - // Required pure virtual function (not used) - virtual void execute() override; - - /// Initialize the System and Mesh objects for the solution being read - virtual void initialSetup() override; - - const std::vector & variableNames() const; - - bool isVariableNodal(const std::string & var_name) const; - - static MooseEnum weightingType() - { - return MooseEnum("found_first=1 average=2 smallest_element_id=4 largest_element_id=8", - "found_first"); - } - - /** - * Return the spatial dimension of the mesh file - * @return The spatial dimension of the mesh file - */ - unsigned int getMeshFileDimension() const { return _mesh->spatial_dimension(); } - - /** - * Return the name of the mesh file this object read the solution from - */ - const std::string getMeshFileName() const { return _mesh_file; } - - /** - * Get the map from block name to block ID. Only works for ExodusII files. - * - * @return Map from block name to block ID - */ - const std::map & getBlockNamesToIds() const - { - return _block_name_to_id; - } - - /** - * Get the map from block id to block name. Only works for ExodusII files. - * - * @return Map from block ID to block name - */ - const std::map & getBlockIdsToNames() const - { - return _block_id_to_name; - } - - /** - * Get the type of file that was read - * @return Returns a MooseEnum that specifies the type of file read - */ - MooseEnum getSolutionFileType() const; -protected: /** - * Method for reading XDA mesh and equation systems file(s) - * This method is called by the constructor when 'file_type = xda' is set - * in the input file. + * Get the time at which to sample the solution */ - void readXda(); + virtual Real solutionSampleTime() override; - /** - * Method for reading an ExodusII file, which is called when - * 'file_type = exodusII is set in the input file. - */ - void readExodusII(); - - /** - * Method for extracting value of solution based on the DOF, - * this is called by the public overloaded function that accept - * a node or element pointer. - * @param dof_index The DOF of the solution desired - * @return The solution at the given DOF - */ - virtual Real directValue(dof_id_type dof_index) const; - - /** - * Updates the times for interpolating ExodusII data - */ - void updateExodusTimeInterpolation(); - - /** - * Updates the time indices to interpolate between for ExodusII data - */ - bool updateExodusBracketingTimeIndices(); - - /** - * A wrapper method for calling the various MeshFunctions used for reading the data - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - Real evalMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method for calling the various MeshFunctions that calls the mesh function - * functionality for evaluating discontinuous shape functions near a face (where it's multivalued) - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - std::map - evalMultiValuedMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method interfacing with the libMesh mesh function for evaluating the gradient - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - libMesh::RealGradient - evalMeshFunctionGradient(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method interfacing with the libMesh mesh function that calls the gradient - * functionality for evaluating potentially discontinuous gradients at element's faces (where it's - * multivalued) - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - std::map evalMultiValuedMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * Read block ID map from the ExodusII file - */ - void readBlockIdMapFromExodusII(); - - /// File type to read (0 = xda; 1 = ExodusII) - MooseEnum _file_type; - - /// The XDA or ExodusII file that is being read - std::string _mesh_file; - - /// The XDA file that contians the EquationSystems data (xda only) - std::string _es_file; - - /// The system name to extract from the XDA file (xda only) - std::string _system_name; - - /// A list of variables to extract from the read system - std::vector _system_variables; - - /// Stores the local index need by MeshFunction - std::map _local_variable_index; - - /// Stores names of nodal variables - std::vector _nodal_variables; - - /// Stores names of elemental variables - std::vector _elemental_variables; - - /// Stores names of scalar variables - std::vector _scalar_variables; - - /// Current ExodusII time index - int _exodus_time_index; - - /// Flag for triggering interpolation of ExodusII data - bool _interpolate_times; - - /// Pointer the libMesh::mesh object - std::unique_ptr _mesh; - - /// Pointer to the libMesh::EquationSystems object - std::unique_ptr _es; - - /// Pointer libMesh::System class storing the read solution - libMesh::System * _system; - - /// Pointer the libMesh::MeshFunction object that the read data is stored - std::unique_ptr _mesh_function; - - /// Pointer to the libMesh::ExodusII used to read the files - std::unique_ptr _exodusII_io; - - /// Pointer to the serial solution vector - std::unique_ptr> _serialized_solution; - - /// Pointer to second libMesh::EquationSystems object, used for interpolation - std::unique_ptr _es2; - - /// Pointer to a second libMesh::System object, used for interpolation - libMesh::System * _system2; - - /// Pointer to second libMesh::MeshFuntion, used for interpolation - std::unique_ptr _mesh_function2; - - /// Pointer to second serial solution, used for interpolation - std::unique_ptr> _serialized_solution2; - - /// Time in the current simulation at which the solution interpolation was last updated - Real _interpolation_time; - - /// Interpolation weight factor - Real _interpolation_factor; - - /// The times available in the ExodusII file - const std::vector * _exodus_times; - - /// Time index 1, used for interpolation - int _exodus_index1; - - /// Time index 2, used for interpolation - int _exodus_index2; - - /// Scale parameter - std::vector _scale; - - /// scale_multiplier parameter - std::vector _scale_multiplier; - - /// Translation - std::vector _translation; - - /// vector about which to rotate - RealVectorValue _rotation0_vector; - - /// angle (in degrees) which to rotate through about vector _rotation0_vector - Real _rotation0_angle; - - /// Rotation matrix that performs the "_rotation0_angle about rotation0_vector" - RealTensorValue _r0; - - /// vector about which to rotate - RealVectorValue _rotation1_vector; - - /// angle (in degrees) which to rotate through about vector _rotation1_vector - Real _rotation1_angle; - - /// Rotation matrix that performs the "_rotation1_angle about rotation1_vector" - RealTensorValue _r1; - - /// transformations (rotations, translation, scales) are performed in this order - MultiMooseEnum _transformation_order; - - /// True if initial_setup has executed - bool _initialized; - - /// Map from block ID to block names. Read from the ExodusII file - std::map _block_name_to_id; - - /// Map from block names to block IDs. Read from the ExodusII file - std::map _block_id_to_name; - - /// function parser object for the resudual and on-diagonal Jacobian + /// function parser object for transforming the solution sample time SymFunctionPtr _time_transformation; - - /// Time at which to sample the solution - Real _solution_time; - -private: - static Threads::spin_mutex _solution_user_object_mutex; }; diff --git a/framework/include/userobjects/SolutionUserObjectBase.h b/framework/include/userobjects/SolutionUserObjectBase.h new file mode 100644 index 000000000000..797a918f5982 --- /dev/null +++ b/framework/include/userobjects/SolutionUserObjectBase.h @@ -0,0 +1,514 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +// MOOSE includes +#include "GeneralUserObject.h" + +// Forward declarations +namespace libMesh +{ +class ExodusII_IO; +class EquationSystems; +class System; +class MeshFunction; +template +class NumericVector; +} + +/** + * User object that reads an existing solution from an input file and + * uses it in the current simulation. + */ +class SolutionUserObjectBase : public GeneralUserObject +{ +public: + static InputParameters validParams(); + + SolutionUserObjectBase(const InputParameters & parameters); + + /** + * Get the time at which to sample the solution + */ + virtual Real solutionSampleTime() = 0; + + /** + * When reading ExodusII files, this will update the interpolation times + */ + virtual void timestepSetup() override; + + /** + * Returns the local index for a given variable name + * @param var_name The name of the variable for which the index is located + * @return The local index of the variable + */ + unsigned int getLocalVarIndex(const std::string & var_name) const; + + /** + * Returns a value at a specific location and variable checking for multiple values and weighting + * these values to + * obtain a single unique value (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValueWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type = weightingType(), + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValue(Real t, + const Point & p, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable for cases where the solution is + * multivalued at element faces + * Use pointValue for continuous shape functions or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map + discontinuousPointValue(Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable for cases where the solution is + * multivalued at element faces + * Use pointValue for continuous shape functions or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map + discontinuousPointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable checking for multiple values and + * weighting these values to + * obtain a single unique value (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradientWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type = weightingType(), + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradient(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradient(Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable for cases where the gradient is + * multivalued (e.g. at element faces) + * Use pointValueGradient for continuous gradients or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map discontinuousPointValueGradient( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable for cases where the gradient is + * multivalued (e.g. at element faces) + * Use pointValueGradient for continuous gradients or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map discontinuousPointValueGradient( + Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Return a value directly from a Node + * @param node A pointer to the node at which a value is desired + * @param var_name The variable from which to extract a value + * @return The desired value for the given node and variable name + */ + Real directValue(const Node * node, const std::string & var_name) const; + + /** + * Return a value from the centroid of an element + * @param elem A pointer to the element at which a value is desired + * @param var_name The variable from which to extract a value + * @return The desired value for the given element and variable name + */ + Real directValue(const Elem * elem, const std::string & var_name) const; + + /** + * Returns a value of a global variable + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param var_name The variable from which to extract a value + * @return The desired value for the given variable + */ + Real scalarValue(Real t, const std::string & var_name) const; + + // Required pure virtual function (not used) + virtual void initialize() override; + + // Required pure virtual function (not used) + virtual void finalize() override; + + // Required pure virtual function (not used) + virtual void execute() override; + + /// Initialize the System and Mesh objects for the solution being read + virtual void initialSetup() override; + + const std::vector & variableNames() const; + + bool isVariableNodal(const std::string & var_name) const; + + static MooseEnum weightingType() + { + return MooseEnum("found_first=1 average=2 smallest_element_id=4 largest_element_id=8", + "found_first"); + } + + /** + * Return the spatial dimension of the mesh file + * @return The spatial dimension of the mesh file + */ + unsigned int getMeshFileDimension() const { return _mesh->spatial_dimension(); } + + /** + * Return the name of the mesh file this object read the solution from + */ + const std::string getMeshFileName() const { return _mesh_file; } + + /** + * Get the map from block name to block ID. Only works for ExodusII files. + * + * @return Map from block name to block ID + */ + const std::map & getBlockNamesToIds() const + { + return _block_name_to_id; + } + + /** + * Get the map from block id to block name. Only works for ExodusII files. + * + * @return Map from block ID to block name + */ + const std::map & getBlockIdsToNames() const + { + return _block_id_to_name; + } + + /** + * Get the type of file that was read + * @return Returns a MooseEnum that specifies the type of file read + */ + MooseEnum getSolutionFileType() const; + +protected: + /** + * Method for reading XDA mesh and equation systems file(s) + * This method is called by the constructor when 'file_type = xda' is set + * in the input file. + */ + void readXda(); + + /** + * Method for reading an ExodusII file, which is called when + * 'file_type = exodusII is set in the input file. + */ + void readExodusII(); + + /** + * Method for extracting value of solution based on the DOF, + * this is called by the public overloaded function that accept + * a node or element pointer. + * @param dof_index The DOF of the solution desired + * @return The solution at the given DOF + */ + virtual Real directValue(dof_id_type dof_index) const; + + /** + * Updates the times for interpolating ExodusII data + */ + void updateExodusTimeInterpolation(); + + /** + * Updates the time indices to interpolate between for ExodusII data + */ + bool updateExodusBracketingTimeIndices(); + + /** + * A wrapper method for calling the various MeshFunctions used for reading the data + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + Real evalMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method for calling the various MeshFunctions that calls the mesh function + * functionality for evaluating discontinuous shape functions near a face (where it's multivalued) + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + std::map + evalMultiValuedMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method interfacing with the libMesh mesh function for evaluating the gradient + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + libMesh::RealGradient + evalMeshFunctionGradient(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method interfacing with the libMesh mesh function that calls the gradient + * functionality for evaluating potentially discontinuous gradients at element's faces (where it's + * multivalued) + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + std::map evalMultiValuedMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * Read block ID map from the ExodusII file + */ + void readBlockIdMapFromExodusII(); + + /// File type to read (0 = xda; 1 = ExodusII) + MooseEnum _file_type; + + /// The XDA or ExodusII file that is being read + std::string _mesh_file; + + /// The XDA file that contians the EquationSystems data (xda only) + std::string _es_file; + + /// The system name to extract from the XDA file (xda only) + std::string _system_name; + + /// A list of variables to extract from the read system + std::vector _system_variables; + + /// Stores the local index need by MeshFunction + std::map _local_variable_index; + + /// Stores names of nodal variables + std::vector _nodal_variables; + + /// Stores names of elemental variables + std::vector _elemental_variables; + + /// Stores names of scalar variables + std::vector _scalar_variables; + + /// Current ExodusII time index + int _exodus_time_index; + + /// Flag for triggering interpolation of ExodusII data + bool _interpolate_times; + + /// Pointer the libMesh::mesh object + std::unique_ptr _mesh; + + /// Pointer to the libMesh::EquationSystems object + std::unique_ptr _es; + + /// Pointer libMesh::System class storing the read solution + libMesh::System * _system; + + /// Pointer the libMesh::MeshFunction object that the read data is stored + std::unique_ptr _mesh_function; + + /// Pointer to the libMesh::ExodusII used to read the files + std::unique_ptr _exodusII_io; + + /// Pointer to the serial solution vector + std::unique_ptr> _serialized_solution; + + /// Pointer to second libMesh::EquationSystems object, used for interpolation + std::unique_ptr _es2; + + /// Pointer to a second libMesh::System object, used for interpolation + libMesh::System * _system2; + + /// Pointer to second libMesh::MeshFuntion, used for interpolation + std::unique_ptr _mesh_function2; + + /// Pointer to second serial solution, used for interpolation + std::unique_ptr> _serialized_solution2; + + /// Time in the current simulation at which the solution interpolation was last updated + Real _interpolation_time; + + /// Interpolation weight factor + Real _interpolation_factor; + + /// The times available in the ExodusII file + const std::vector * _exodus_times; + + /// Time index 1, used for interpolation + int _exodus_index1; + + /// Time index 2, used for interpolation + int _exodus_index2; + + /// Scale parameter + std::vector _scale; + + /// scale_multiplier parameter + std::vector _scale_multiplier; + + /// Translation + std::vector _translation; + + /// vector about which to rotate + RealVectorValue _rotation0_vector; + + /// angle (in degrees) which to rotate through about vector _rotation0_vector + Real _rotation0_angle; + + /// Rotation matrix that performs the "_rotation0_angle about rotation0_vector" + RealTensorValue _r0; + + /// vector about which to rotate + RealVectorValue _rotation1_vector; + + /// angle (in degrees) which to rotate through about vector _rotation1_vector + Real _rotation1_angle; + + /// Rotation matrix that performs the "_rotation1_angle about rotation1_vector" + RealTensorValue _r1; + + /// transformations (rotations, translation, scales) are performed in this order + MultiMooseEnum _transformation_order; + + /// True if initial_setup has executed + bool _initialized; + + /// Map from block ID to block names. Read from the ExodusII file + std::map _block_name_to_id; + + /// Map from block names to block IDs. Read from the ExodusII file + std::map _block_id_to_name; + +private: + static Threads::spin_mutex _solution_user_object_mutex; +}; diff --git a/framework/src/auxkernels/SolutionAux.C b/framework/src/auxkernels/SolutionAux.C index b3330958e354..f2ac0cc1fd24 100644 --- a/framework/src/auxkernels/SolutionAux.C +++ b/framework/src/auxkernels/SolutionAux.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionAux.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", SolutionAux); @@ -39,7 +39,7 @@ SolutionAux::validParams() SolutionAux::SolutionAux(const InputParameters & parameters) : AuxKernel(parameters), - _solution_object(getUserObject("solution")), + _solution_object(getUserObject("solution")), _direct(getParam("direct")), _scale_factor(getParam("scale_factor")), _add_factor(getParam("add_factor")) diff --git a/framework/src/auxkernels/SolutionScalarAux.C b/framework/src/auxkernels/SolutionScalarAux.C index 2bb3acc175c1..24fa66066828 100644 --- a/framework/src/auxkernels/SolutionScalarAux.C +++ b/framework/src/auxkernels/SolutionScalarAux.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionScalarAux.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", SolutionScalarAux); @@ -35,7 +35,7 @@ SolutionScalarAux::validParams() SolutionScalarAux::SolutionScalarAux(const InputParameters & parameters) : AuxScalarKernel(parameters), - _solution_object(getUserObject("solution")), + _solution_object(getUserObject("solution")), _scale_factor(getParam("scale_factor")), _add_factor(getParam("add_factor")) { diff --git a/framework/src/functions/Axisymmetric2D3DSolutionFunction.C b/framework/src/functions/Axisymmetric2D3DSolutionFunction.C index e88636b85cc0..e3a42e5f325c 100644 --- a/framework/src/functions/Axisymmetric2D3DSolutionFunction.C +++ b/framework/src/functions/Axisymmetric2D3DSolutionFunction.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "Axisymmetric2D3DSolutionFunction.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", Axisymmetric2D3DSolutionFunction); @@ -107,7 +107,7 @@ Axisymmetric2D3DSolutionFunction::initialSetup() // Get a pointer to the SolutionUserObject. A pointer is used because the UserObject is not // available during the // construction of the function - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); // If 'from_variable' is not specified, get the value from the SolutionUserObject if (_var_names.size() == 0) diff --git a/framework/src/functions/SolutionFunction.C b/framework/src/functions/SolutionFunction.C index e5dd7937468a..f67b14333b8e 100644 --- a/framework/src/functions/SolutionFunction.C +++ b/framework/src/functions/SolutionFunction.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionFunction.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" registerMooseObject("MooseApp", SolutionFunction); @@ -57,7 +57,7 @@ SolutionFunction::initialSetup() // Get a pointer to the SolutionUserObject. A pointer is used because the UserObject is not // available during the // construction of the function - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); std::string var_name; diff --git a/framework/src/ics/ScalarSolutionIC.C b/framework/src/ics/ScalarSolutionIC.C index 09323ca8a206..5e3243ce53f8 100644 --- a/framework/src/ics/ScalarSolutionIC.C +++ b/framework/src/ics/ScalarSolutionIC.C @@ -8,7 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "ScalarSolutionIC.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" registerMooseObject("MooseApp", ScalarSolutionIC); @@ -34,7 +34,7 @@ ScalarSolutionIC::validParams() ScalarSolutionIC::ScalarSolutionIC(const InputParameters & parameters) : ScalarInitialCondition(parameters), - _solution_object(getUserObject("solution_uo")), + _solution_object(getUserObject("solution_uo")), _solution_object_var_name(getParam("from_variable")) { } diff --git a/framework/src/ics/SolutionIC.C b/framework/src/ics/SolutionIC.C index ab6d86352e44..90448860a672 100644 --- a/framework/src/ics/SolutionIC.C +++ b/framework/src/ics/SolutionIC.C @@ -8,7 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "SolutionIC.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" #include "SystemBase.h" @@ -31,7 +31,7 @@ SolutionIC::validParams() SolutionIC::SolutionIC(const InputParameters & parameters) : InitialCondition(parameters), - _solution_object(getUserObject("solution_uo")), + _solution_object(getUserObject("solution_uo")), _solution_object_var_name(getParam("from_variable")) { } diff --git a/framework/src/userobjects/SolutionUserObject.C b/framework/src/userobjects/SolutionUserObject.C index c4dafb809664..a2bb3345fdd8 100644 --- a/framework/src/userobjects/SolutionUserObject.C +++ b/framework/src/userobjects/SolutionUserObject.C @@ -36,58 +36,9 @@ InputParameters SolutionUserObject::validParams() { // Get the input parameters from the parent class - InputParameters params = GeneralUserObject::validParams(); + InputParameters params = SolutionUserObjectBase::validParams(); params += FunctionParserUtils::validParams(); - // Add required parameters - params.addRequiredParam( - "mesh", "The name of the mesh file (must be xda/xdr or exodusII file)."); - params.addParam>( - "system_variables", - std::vector(), - "The name of the nodal and elemental variables from the file you want to use for values"); - - // When using XDA/XDR files the following must be defined - params.addParam( - "es", - "", - "The name of the file holding the equation system info in xda/xdr format (xda/xdr only)."); - params.addParam( - "system", - "nl0", - "The name of the system to pull values out of (xda/xdr only). The default name for the " - "nonlinear system is 'nl0', auxiliary system is 'aux0'"); - - // When using ExodusII a specific time is extracted - params.addParam("timestep", - "Index of the single timestep used or \"LATEST\" for " - "the last timestep (exodusII only). If not supplied, " - "time interpolation will occur."); - - // Add ability to perform coordinate transformation: scale, factor - params.addParam>( - "scale", std::vector(LIBMESH_DIM, 1), "Scale factor for points in the simulation"); - params.addParam>("scale_multiplier", - std::vector(LIBMESH_DIM, 1), - "Scale multiplying factor for points in the simulation"); - params.addParam>("translation", - std::vector(LIBMESH_DIM, 0), - "Translation factors for x,y,z coordinates of the simulation"); - params.addParam("rotation0_vector", - RealVectorValue(0, 0, 1), - "Vector about which to rotate points of the simulation."); - params.addParam( - "rotation0_angle", - 0.0, - "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector."); - params.addParam("rotation1_vector", - RealVectorValue(0, 0, 1), - "Vector about which to rotate points of the simulation."); - params.addParam( - "rotation1_angle", - 0.0, - "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector."); - params.addCustomTypeParam( "time_transformation", "t", @@ -95,86 +46,12 @@ SolutionUserObject::validParams() "Expression to transform from current simulation time to time at " "which to sample the solution."); - // following lines build the default_transformation_order - MultiMooseEnum default_transformation_order( - "rotation0 translation scale rotation1 scale_multiplier", "translation scale"); - params.addParam( - "transformation_order", - default_transformation_order, - "The order to perform the operations in. Define R0 to be the rotation matrix encoded by " - "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the " - "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, " - "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then " - "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObject at point x " - "in the simulation are the variable values at point p in the mesh."); - params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle " - "rotation1_angle transformation_order", - "Coordinate system transformation"); - params.addClassDescription("Reads a variable from a mesh in one simulation to another"); return params; } -// Static mutex definition -Threads::spin_mutex SolutionUserObject::_solution_user_object_mutex; - SolutionUserObject::SolutionUserObject(const InputParameters & parameters) - : GeneralUserObject(parameters), - FunctionParserUtils(parameters), - _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")), - _mesh_file(getParam("mesh")), - _es_file(getParam("es")), - _system_name(getParam("system")), - _system_variables(getParam>("system_variables")), - _exodus_time_index(-1), - _interpolate_times(false), - _system(nullptr), - _system2(nullptr), - _interpolation_time(0.0), - _interpolation_factor(0.0), - _exodus_times(nullptr), - _exodus_index1(-1), - _exodus_index2(-1), - _scale(getParam>("scale")), - _scale_multiplier(getParam>("scale_multiplier")), - _translation(getParam>("translation")), - _rotation0_vector(getParam("rotation0_vector")), - _rotation0_angle(getParam("rotation0_angle")), - _r0(RealTensorValue()), - _rotation1_vector(getParam("rotation1_vector")), - _rotation1_angle(getParam("rotation1_angle")), - _r1(RealTensorValue()), - _transformation_order(getParam("transformation_order")), - _initialized(false) + : SolutionUserObjectBase(parameters), FunctionParserUtils(parameters) { - // form rotation matrices with the specified angles - Real halfPi = libMesh::pi / 2.0; - Real a; - Real b; - - a = std::cos(halfPi * -_rotation0_angle / 90.0); - b = std::sin(halfPi * -_rotation0_angle / 90.0); - // the following is an anticlockwise rotation about z - RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1); - // form the rotation matrix that will take rotation0_vector to the z axis - RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector); - // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points - // back - _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z); - - a = std::cos(halfPi * -_rotation1_angle / 90.0); - b = std::sin(halfPi * -_rotation1_angle / 90.0); - // the following is an anticlockwise rotation about z - RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1); - // form the rotation matrix that will take rotation1_vector to the z axis - RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector); - // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points - // back - _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z); - - if (isParamValid("timestep") && getParam("timestep") == "-1") - mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply " - "remove this parameter altogether for interpolation"); - // setup parsed expression for the time transformation _time_transformation = std::make_shared(); setParserFeatureFlags(_time_transformation); @@ -184,1100 +61,13 @@ SolutionUserObject::SolutionUserObject(const InputParameters & parameters) if (_time_transformation->Parse(expression, "t") >= 0) mooseError("Invalid parsed function\n", expression, "\n", _time_transformation->ErrorMsg()); - // only parameter is time + // the only parameter is time _func_params.resize(1); } -SolutionUserObject::~SolutionUserObject() {} - -void -SolutionUserObject::readXda() -{ - if (!isParamSetByUser("es")) - paramError("es", "Equation system file (.xda or .xdr) should have been specified"); - - // Check that the required files exist - MooseUtils::checkFileReadable(_es_file); - MooseUtils::checkFileReadable(_mesh_file); - - // Read the libmesh::mesh from the xda file - _mesh->read(_mesh_file); - - // Create the libmesh::EquationSystems - _es = std::make_unique(*_mesh); - - // Use new read syntax (binary) - if (_file_type == "xdr") - _es->read(_es_file, - libMesh::DECODE, - EquationSystems::READ_HEADER | EquationSystems::READ_DATA | - EquationSystems::READ_ADDITIONAL_DATA); - - // Use new read syntax - else if (_file_type == "xda") - _es->read(_es_file, - libMesh::READ, - EquationSystems::READ_HEADER | EquationSystems::READ_DATA | - EquationSystems::READ_ADDITIONAL_DATA); - - // This should never occur, just in case produce an error - else - mooseError("Failed to determine proper read method for XDA/XDR equation system file: ", - _es_file); - - // Update and store the EquationSystems name locally - _es->update(); - _system = &_es->get_system(_system_name); -} - -void -SolutionUserObject::readExodusII() -{ - // Define a default system name - if (_system_name == "") - _system_name = "SolutionUserObjectSystem"; - - // Read the Exodus file - _exodusII_io = std::make_unique(*_mesh); - _exodusII_io->read(_mesh_file); - readBlockIdMapFromExodusII(); - _exodus_times = &_exodusII_io->get_time_steps(); - - if (isParamValid("timestep")) - { - std::string s_timestep = getParam("timestep"); - int n_steps = _exodusII_io->get_num_time_steps(); - if (s_timestep == "LATEST") - _exodus_time_index = n_steps; - else - { - std::istringstream ss(s_timestep); - if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps) - mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer " - "less than ", - n_steps, - ", received ", - s_timestep); - } - } - else - // Interpolate between times rather than using values from a set timestep - _interpolate_times = true; - - // Check that the number of time steps is valid - int num_exo_times = _exodus_times->size(); - if (num_exo_times == 0) - mooseError("In SolutionUserObject, exodus file contains no timesteps."); - - // Account for parallel mesh - if (dynamic_cast(_mesh.get())) - { - _mesh->allow_renumbering(true); - _mesh->prepare_for_use(); - } - else - { - _mesh->allow_renumbering(false); - _mesh->prepare_for_use(); - } - - // Create EquationSystems object for solution - _es = std::make_unique(*_mesh); - _es->add_system(_system_name); - _system = &_es->get_system(_system_name); - - // Get the variable name lists as set; these need to be sets to perform set_intersection - const std::vector & all_nodal(_exodusII_io->get_nodal_var_names()); - const std::vector & all_elemental(_exodusII_io->get_elem_var_names()); - const std::vector & all_scalar(_exodusII_io->get_global_var_names()); - - // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if - // provided - // This function could be called more than once, so clear the member variables so we don't keep - // adding to the vectors - _nodal_variables.clear(); - _elemental_variables.clear(); - _scalar_variables.clear(); - if (!_system_variables.empty()) - { - for (const auto & var_name : _system_variables) - { - if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end()) - _nodal_variables.push_back(var_name); - if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end()) - _elemental_variables.push_back(var_name); - if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end()) - // Check if the scalar matches any field variables, and ignore the var if it does. This - // means its a Postprocessor. - if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == - _nodal_variables.end() && - std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == - _elemental_variables.end()) - _scalar_variables.push_back(var_name); - } - } - else - { - _nodal_variables = all_nodal; - _elemental_variables = all_elemental; - - for (auto var_name : all_scalar) - // Check if the scalar matches any field variables, and ignore the var if it does. This means - // its a Postprocessor. - if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == - _nodal_variables.end() && - std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == - _elemental_variables.end()) - _scalar_variables.push_back(var_name); - } - - // Add the variables to the system - for (const auto & var_name : _nodal_variables) - _system->add_variable(var_name, FIRST); - - for (const auto & var_name : _elemental_variables) - _system->add_variable(var_name, CONSTANT, MONOMIAL); - - for (const auto & var_name : _scalar_variables) - _system->add_variable(var_name, FIRST, SCALAR); - - // Initialize the equations systems - _es->init(); - - // Interpolate times - if (_interpolate_times) - { - // Create a second equation system - _es2 = std::make_unique(*_mesh); - _es2->add_system(_system_name); - _system2 = &_es2->get_system(_system_name); - - // Add the variables to the system - for (const auto & var_name : _nodal_variables) - _system2->add_variable(var_name, FIRST); - - for (const auto & var_name : _elemental_variables) - _system2->add_variable(var_name, CONSTANT, MONOMIAL); - - for (const auto & var_name : _scalar_variables) - _system2->add_variable(var_name, FIRST, SCALAR); - - // Initialize - _es2->init(); - - // Update the times for interpolation (initially start at 0) - updateExodusBracketingTimeIndices(); - - // Copy the solutions from the first system - for (const auto & var_name : _nodal_variables) - { - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); - _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - } - - for (const auto & var_name : _elemental_variables) - { - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); - _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - } - - if (_scalar_variables.size() > 0) - { - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); - _exodusII_io->copy_scalar_solution( - *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); - } - - // Update the systems - _system->update(); - _es->update(); - _system2->update(); - _es2->update(); - } - - // Non-interpolated times - else - { - if (_exodus_time_index > num_exo_times) - mooseError("In SolutionUserObject, timestep = ", - _exodus_time_index, - ", but there are only ", - num_exo_times, - " time steps."); - - // Copy the values from the ExodusII file - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index); - - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index); - - if (!_scalar_variables.empty()) - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_time_index); - - // Update the equations systems - _system->update(); - _es->update(); - } -} - Real -SolutionUserObject::directValue(const Node * node, const std::string & var_name) const +SolutionUserObject::solutionSampleTime() { - // Get the libmesh variable and system numbers - unsigned int var_num = _system->variable_number(var_name); - unsigned int sys_num = _system->number(); - - // Get the node id and associated dof - dof_id_type node_id = node->id(); - const Node & sys_node = _system->get_mesh().node_ref(node_id); - mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0, - "Variable " << var_name << " has no DoFs on node " << sys_node.id()); - dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0); - - // Return the desired value for the dof - return directValue(dof_id); -} - -Real -SolutionUserObject::directValue(const Elem * elem, const std::string & var_name) const -{ - // Get the libmesh variable and system numbers - unsigned int var_num = _system->variable_number(var_name); - unsigned int sys_num = _system->number(); - - // Get the element id and associated dof - dof_id_type elem_id = elem->id(); - const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id); - mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0, - "Variable " << var_name << " has no DoFs on element " << sys_elem.id()); - dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0); - - // Return the desired value - return directValue(dof_id); -} - -void -SolutionUserObject::initialize() -{ -} - -void -SolutionUserObject::finalize() -{ -} - -void -SolutionUserObject::timestepSetup() -{ - // Update time interpolation for ExodusII solution - if (_file_type == 1 && _interpolate_times) - updateExodusTimeInterpolation(); -} - -void -SolutionUserObject::execute() -{ -} - -void -SolutionUserObject::initialSetup() -{ - - // Make sure this only happens once - if (_initialized) - return; - - // Create a libmesh::Mesh object for storing the loaded data. - // Several aspects of SolutionUserObject won't work with a DistributedMesh: - // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel. - // .) We don't know if directValue will be used, which may request - // a value on a Node we don't have. - // We force the Mesh used here to be a ReplicatedMesh. - _mesh = std::make_unique(_communicator); - - // ExodusII mesh file supplied - if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) || - MooseUtils::hasExtension(_mesh_file, "exo", true)) - { - _file_type = "exodusII"; - readExodusII(); - } - - // XDA mesh file supplied - else if (MooseUtils::hasExtension(_mesh_file, "xda")) - { - _file_type = "xda"; - readXda(); - } - - // XDR mesh file supplied - else if (MooseUtils::hasExtension(_mesh_file, "xdr")) - { - _file_type = "xdr"; - readXda(); - } - - // Produce an error for an unknown file type - else - mooseError("In SolutionUserObject, invalid file type: only .xda, .xdr, .exo and .e supported"); - - // Initialize the serial solution vector - _serialized_solution = NumericVector::build(_communicator); - _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL); - - // Pull down a full copy of this vector on every processor so we can get values in parallel - _system->solution->localize(*_serialized_solution); - - // Vector of variable numbers to apply the MeshFunction to - std::vector var_nums; - - // If no variables were given, use all of them - if (_system_variables.empty()) - { - _system->get_all_variable_numbers(var_nums); - for (const auto & var_num : var_nums) - _system_variables.push_back(_system->variable_name(var_num)); - } - - // Otherwise, gather the numbers for the variables given - else - { - for (const auto & var_name : _system_variables) - var_nums.push_back(_system->variable_number(var_name)); - } - - // Create the MeshFunction for working with the solution data - _mesh_function = std::make_unique( - *_es, *_serialized_solution, _system->get_dof_map(), var_nums); - _mesh_function->init(); - - // Tell the MeshFunctions that we might be querying them outside the - // mesh, so we can handle any errors at the MOOSE rather than at the - // libMesh level. - DenseVector default_values; - _mesh_function->enable_out_of_mesh_mode(default_values); - - // Build second MeshFunction for interpolation - if (_interpolate_times) - { - // Need to pull down a full copy of this vector on every processor so we can get values in - // parallel - _serialized_solution2 = NumericVector::build(_communicator); - _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL); - _system2->solution->localize(*_serialized_solution2); - - // Create the MeshFunction for the second copy of the data - _mesh_function2 = std::make_unique( - *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums); - _mesh_function2->init(); - _mesh_function2->enable_out_of_mesh_mode(default_values); - } - - // Populate the MeshFunction variable index - for (unsigned int i = 0; i < _system_variables.size(); ++i) - { - std::string name = _system_variables[i]; - _local_variable_index[name] = i; - } - - // Set initialization flag - _initialized = true; -} - -MooseEnum -SolutionUserObject::getSolutionFileType() const -{ - return _file_type; -} - -void -SolutionUserObject::updateExodusTimeInterpolation() -{ - if (_t != _interpolation_time) - { - if (updateExodusBracketingTimeIndices()) - { - - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); - if (_scalar_variables.size() > 0) - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); - - _system->update(); - _es->update(); - _system->solution->localize(*_serialized_solution); - - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - if (_scalar_variables.size() > 0) - _exodusII_io->copy_scalar_solution( - *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); - - _system2->update(); - _es2->update(); - _system2->solution->localize(*_serialized_solution2); - } - _interpolation_time = _t; - } -} - -bool -SolutionUserObject::updateExodusBracketingTimeIndices() -{ - if (_file_type != 1) - mooseError( - "In SolutionUserObject, getTimeInterpolationData only applicable for exodusII file type"); - - int old_index1 = _exodus_index1; - int old_index2 = _exodus_index2; - - int num_exo_times = _exodus_times->size(); - _func_params[0] = _t; - Real solution_time = evaluate(_time_transformation); - - if (solution_time < (*_exodus_times)[0]) - { - _exodus_index1 = 0; - _exodus_index2 = 0; - _interpolation_factor = 0.0; - } - else - { - for (int i = 0; i < num_exo_times - 1; ++i) - { - if (solution_time <= (*_exodus_times)[i + 1]) - { - _exodus_index1 = i; - _exodus_index2 = i + 1; - _interpolation_factor = - (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); - break; - } - else if (i == num_exo_times - 2) - { - _exodus_index1 = num_exo_times - 1; - _exodus_index2 = num_exo_times - 1; - _interpolation_factor = 1.0; - break; - } - } - } - - bool indices_modified(false); - - if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2) - indices_modified = true; - - return indices_modified; -} - -unsigned int -SolutionUserObject::getLocalVarIndex(const std::string & var_name) const -{ - // Extract the variable index for the MeshFunction(s) - std::map::const_iterator it = _local_variable_index.find(var_name); - if (it == _local_variable_index.end()) - mooseError("Value requested for nonexistent variable '", - var_name, - "' in the '", - name(), - "' SolutionUserObject.\nSystem selected: ", - _system_name, - "\nAvailable variables:\n", - ConsoleUtils::formatString(Moose::stringify(_system_variables), "")); - return it->second; -} - -Real -SolutionUserObject::pointValueWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type, - const std::set * subdomain_ids) const -{ - // first check if the FE type is continuous because in that case the value is - // unique and we can take a short cut, the default weighting_type found_first also - // shortcuts out - auto family = - _fe_problem - .getVariable( - _tid, var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD) - .feType() - .family; - - if (weighting_type == 1 || - (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC)) - return pointValue(t, p, var_name, subdomain_ids); - - // the shape function is discontinuous so we need to compute a suitable unique value - std::map values = discontinuousPointValue(t, p, var_name); - switch (weighting_type) - { - case 2: - { - Real average = 0.0; - for (auto & v : values) - average += v.second; - return average / Real(values.size()); - } - case 4: - { - Real smallest_elem_id_value = std::numeric_limits::max(); - dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); - for (auto & v : values) - if (v.first->id() < smallest_elem_id) - { - smallest_elem_id = v.first->id(); - smallest_elem_id_value = v.second; - } - return smallest_elem_id_value; - } - case 8: - { - Real largest_elem_id_value = std::numeric_limits::lowest(); - dof_id_type largest_elem_id = 0; - for (auto & v : values) - if (v.first->id() > largest_elem_id) - { - largest_elem_id = v.first->id(); - largest_elem_id_value = v.second; - } - return largest_elem_id_value; - } - } - - mooseError( - "SolutionUserObject::pointValueWrapper reaches line that it should not be able to reach."); - return 0.0; -} - -Real -SolutionUserObject::pointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return pointValue(t, p, local_var_index, subdomain_ids); -} - -Real -SolutionUserObject::pointValue(Real libmesh_dbg_var(t), - const Point & p, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // Create copy of point - Point pt(p); - - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids); - val = val + (val2 - val) * _interpolation_factor; - } - - return val; -} - -std::map -SolutionUserObject::discontinuousPointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return discontinuousPointValue(t, p, local_var_index, subdomain_ids); -} - -std::map -SolutionUserObject::discontinuousPointValue(Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - std::map map = - evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - std::map map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2); - - if (map.size() != map2.size()) - mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size"); - - // construct the interpolated map - for (auto & k : map) - { - if (map2.find(k.first) == map2.end()) - mooseError( - "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys"); - Real val = k.second; - Real val2 = map2[k.first]; - map[k.first] = val + (val2 - val) * _interpolation_factor; - } - } - - return map; -} - -RealGradient -SolutionUserObject::pointValueGradientWrapper( - Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type, - const std::set * subdomain_ids) const -{ - // the default weighting_type found_first shortcuts out - if (weighting_type == 1) - return pointValueGradient(t, p, var_name, subdomain_ids); - - // the shape function is discontinuous so we need to compute a suitable unique value - std::map values = - discontinuousPointValueGradient(t, p, var_name, subdomain_ids); - switch (weighting_type) - { - case 2: - { - RealGradient average = RealGradient(0.0, 0.0, 0.0); - for (auto & v : values) - average += v.second; - return average / Real(values.size()); - } - case 4: - { - RealGradient smallest_elem_id_value; - dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); - for (auto & v : values) - if (v.first->id() < smallest_elem_id) - { - smallest_elem_id = v.first->id(); - smallest_elem_id_value = v.second; - } - return smallest_elem_id_value; - } - case 8: - { - RealGradient largest_elem_id_value; - dof_id_type largest_elem_id = 0; - for (auto & v : values) - if (v.first->id() > largest_elem_id) - { - largest_elem_id = v.first->id(); - largest_elem_id_value = v.second; - } - return largest_elem_id_value; - } - } - - mooseError("SolutionUserObject::pointValueGradientWrapper reaches line that it should not be " - "able to reach."); - return RealGradient(0.0, 0.0, 0.0); -} - -RealGradient -SolutionUserObject::pointValueGradient(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return pointValueGradient(t, p, local_var_index, subdomain_ids); -} - -RealGradient -SolutionUserObject::pointValueGradient(Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids); - val = val + (val2 - val) * _interpolation_factor; - } - - return val; -} - -std::map -SolutionUserObject::discontinuousPointValueGradient( - Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids); -} - -std::map -SolutionUserObject::discontinuousPointValueGradient( - Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - std::map map = - evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - std::map map2 = - evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - if (map.size() != map2.size()) - mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size"); - - // construct the interpolated map - for (auto & k : map) - { - if (map2.find(k.first) == map2.end()) - mooseError( - "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys"); - RealGradient val = k.second; - RealGradient val2 = map2[k.first]; - map[k.first] = val + (val2 - val) * _interpolation_factor; - } - } - - return map; -} - -Real -SolutionUserObject::directValue(dof_id_type dof_index) const -{ - Real val = (*_serialized_solution)(dof_index); - if (_file_type == 1 && _interpolate_times) - { - Real val2 = (*_serialized_solution2)(dof_index); - val = val + (val2 - val) * _interpolation_factor; - } - return val; -} - -Real -SolutionUserObject::evalMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - DenseVector output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - (*_mesh_function)(p, 0.0, output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - (*_mesh_function2)(p, 0.0, output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - return output(local_var_index); -} - -std::map -SolutionUserObject::evalMultiValuedMeshFunction( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::map> temporary_output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (temporary_output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - - // Fill the actual map that is returned - std::map output; - for (auto & k : temporary_output) - { - mooseAssert(k.second.size() > local_var_index, - "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index " - << local_var_index << " does not exist"); - output[k.first] = k.second(local_var_index); - } - - return output; -} - -RealGradient -SolutionUserObject::evalMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::vector output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->gradient(p, 0.0, output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->gradient(p, 0.0, output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - return output[local_var_index]; -} - -std::map -SolutionUserObject::evalMultiValuedMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::map> temporary_output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (temporary_output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - - // Fill the actual map that is returned - std::map output; - for (auto & k : temporary_output) - { - mooseAssert(k.second.size() > local_var_index, - "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index " - << local_var_index << " does not exist"); - output[k.first] = k.second[local_var_index]; - } - - return output; -} - -const std::vector & -SolutionUserObject::variableNames() const -{ - return _system_variables; -} - -bool -SolutionUserObject::isVariableNodal(const std::string & var_name) const -{ - return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) != - _nodal_variables.end(); -} - -Real -SolutionUserObject::scalarValue(Real /*t*/, const std::string & var_name) const -{ - unsigned int var_num = _system->variable_number(var_name); - const DofMap & dof_map = _system->get_dof_map(); - std::vector dofs; - dof_map.SCALAR_dof_indices(dofs, var_num); - // We can handle only FIRST order scalar variables - return directValue(dofs[0]); -} - -void -SolutionUserObject::readBlockIdMapFromExodusII() -{ -#ifdef LIBMESH_HAVE_EXODUS_API - libMesh::ExodusII_IO_Helper & exio_helper = _exodusII_io->get_exio_helper(); - const auto & id_to_block = exio_helper.id_to_block_names; - _block_name_to_id.clear(); - _block_id_to_name.clear(); - for (const auto & it : id_to_block) - { - _block_name_to_id[it.second] = it.first; - _block_id_to_name[it.first] = it.second; - } -#endif + return evaluate(_time_transformation); } diff --git a/framework/src/userobjects/SolutionUserObjectBase.C b/framework/src/userobjects/SolutionUserObjectBase.C new file mode 100644 index 000000000000..722afa09b599 --- /dev/null +++ b/framework/src/userobjects/SolutionUserObjectBase.C @@ -0,0 +1,1265 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "SolutionUserObjectBase.h" + +// MOOSE includes +#include "ConsoleUtils.h" +#include "MooseError.h" +#include "MooseMesh.h" +#include "MooseUtils.h" +#include "MooseVariableFE.h" +#include "RotationMatrix.h" +#include "Function.h" + +// libMesh includes +#include "libmesh/equation_systems.h" +#include "libmesh/mesh_function.h" +#include "libmesh/numeric_vector.h" +#include "libmesh/nonlinear_implicit_system.h" +#include "libmesh/transient_system.h" +#include "libmesh/parallel_mesh.h" +#include "libmesh/serial_mesh.h" +#include "libmesh/exodusII_io.h" +#include "libmesh/exodusII_io_helper.h" +#include "libmesh/enum_xdr_mode.h" + +InputParameters +SolutionUserObjectBase::validParams() +{ + // Get the input parameters from the parent class + InputParameters params = GeneralUserObject::validParams(); + + // Add required parameters + params.addRequiredParam( + "mesh", "The name of the mesh file (must be xda/xdr or exodusII file)."); + params.addParam>( + "system_variables", + std::vector(), + "The name of the nodal and elemental variables from the file you want to use for values"); + + // When using XDA/XDR files the following must be defined + params.addParam( + "es", + "", + "The name of the file holding the equation system info in xda/xdr format (xda/xdr only)."); + params.addParam( + "system", + "nl0", + "The name of the system to pull values out of (xda/xdr only). The default name for the " + "nonlinear system is 'nl0', auxiliary system is 'aux0'"); + + // When using ExodusII a specific time is extracted + params.addParam("timestep", + "Index of the single timestep used or \"LATEST\" for " + "the last timestep (exodusII only). If not supplied, " + "time interpolation will occur."); + + // Add ability to perform coordinate transformation: scale, factor + params.addParam>( + "scale", std::vector(LIBMESH_DIM, 1), "Scale factor for points in the simulation"); + params.addParam>("scale_multiplier", + std::vector(LIBMESH_DIM, 1), + "Scale multiplying factor for points in the simulation"); + params.addParam>("translation", + std::vector(LIBMESH_DIM, 0), + "Translation factors for x,y,z coordinates of the simulation"); + params.addParam("rotation0_vector", + RealVectorValue(0, 0, 1), + "Vector about which to rotate points of the simulation."); + params.addParam( + "rotation0_angle", + 0.0, + "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector."); + params.addParam("rotation1_vector", + RealVectorValue(0, 0, 1), + "Vector about which to rotate points of the simulation."); + params.addParam( + "rotation1_angle", + 0.0, + "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector."); + + // following lines build the default_transformation_order + MultiMooseEnum default_transformation_order( + "rotation0 translation scale rotation1 scale_multiplier", "translation scale"); + params.addParam( + "transformation_order", + default_transformation_order, + "The order to perform the operations in. Define R0 to be the rotation matrix encoded by " + "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the " + "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, " + "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then " + "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObjectBase at " + "point x " + "in the simulation are the variable values at point p in the mesh."); + params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle " + "rotation1_angle transformation_order", + "Coordinate system transformation"); + params.addClassDescription("Reads a variable from a mesh in one simulation to another"); + return params; +} + +// Static mutex definition +Threads::spin_mutex SolutionUserObjectBase::_solution_user_object_mutex; + +SolutionUserObjectBase::SolutionUserObjectBase(const InputParameters & parameters) + : GeneralUserObject(parameters), + _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")), + _mesh_file(getParam("mesh")), + _es_file(getParam("es")), + _system_name(getParam("system")), + _system_variables(getParam>("system_variables")), + _exodus_time_index(-1), + _interpolate_times(false), + _system(nullptr), + _system2(nullptr), + _interpolation_time(0.0), + _interpolation_factor(0.0), + _exodus_times(nullptr), + _exodus_index1(-1), + _exodus_index2(-1), + _scale(getParam>("scale")), + _scale_multiplier(getParam>("scale_multiplier")), + _translation(getParam>("translation")), + _rotation0_vector(getParam("rotation0_vector")), + _rotation0_angle(getParam("rotation0_angle")), + _r0(RealTensorValue()), + _rotation1_vector(getParam("rotation1_vector")), + _rotation1_angle(getParam("rotation1_angle")), + _r1(RealTensorValue()), + _transformation_order(getParam("transformation_order")), + _initialized(false) +{ + // form rotation matrices with the specified angles + Real halfPi = libMesh::pi / 2.0; + Real a; + Real b; + + a = std::cos(halfPi * -_rotation0_angle / 90.0); + b = std::sin(halfPi * -_rotation0_angle / 90.0); + // the following is an anticlockwise rotation about z + RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1); + // form the rotation matrix that will take rotation0_vector to the z axis + RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector); + // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points + // back + _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z); + + a = std::cos(halfPi * -_rotation1_angle / 90.0); + b = std::sin(halfPi * -_rotation1_angle / 90.0); + // the following is an anticlockwise rotation about z + RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1); + // form the rotation matrix that will take rotation1_vector to the z axis + RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector); + // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points + // back + _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z); + + if (isParamValid("timestep") && getParam("timestep") == "-1") + mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply " + "remove this parameter altogether for interpolation"); +} + +void +SolutionUserObjectBase::readXda() +{ + if (!isParamSetByUser("es")) + paramError("es", "Equation system file (.xda or .xdr) should have been specified"); + + // Check that the required files exist + MooseUtils::checkFileReadable(_es_file); + MooseUtils::checkFileReadable(_mesh_file); + + // Read the libmesh::mesh from the xda file + _mesh->read(_mesh_file); + + // Create the libmesh::EquationSystems + _es = std::make_unique(*_mesh); + + // Use new read syntax (binary) + if (_file_type == "xdr") + _es->read(_es_file, + libMesh::DECODE, + EquationSystems::READ_HEADER | EquationSystems::READ_DATA | + EquationSystems::READ_ADDITIONAL_DATA); + + // Use new read syntax + else if (_file_type == "xda") + _es->read(_es_file, + libMesh::READ, + EquationSystems::READ_HEADER | EquationSystems::READ_DATA | + EquationSystems::READ_ADDITIONAL_DATA); + + // This should never occur, just in case produce an error + else + mooseError("Failed to determine proper read method for XDA/XDR equation system file: ", + _es_file); + + // Update and store the EquationSystems name locally + _es->update(); + _system = &_es->get_system(_system_name); +} + +void +SolutionUserObjectBase::readExodusII() +{ + // Define a default system name + if (_system_name == "") + _system_name = "SolutionUserObjectSystem"; + + // Read the Exodus file + _exodusII_io = std::make_unique(*_mesh); + _exodusII_io->read(_mesh_file); + readBlockIdMapFromExodusII(); + _exodus_times = &_exodusII_io->get_time_steps(); + + if (isParamValid("timestep")) + { + std::string s_timestep = getParam("timestep"); + int n_steps = _exodusII_io->get_num_time_steps(); + if (s_timestep == "LATEST") + _exodus_time_index = n_steps; + else + { + std::istringstream ss(s_timestep); + if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps) + mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer " + "less than ", + n_steps, + ", received ", + s_timestep); + } + } + else + // Interpolate between times rather than using values from a set timestep + _interpolate_times = true; + + // Check that the number of time steps is valid + int num_exo_times = _exodus_times->size(); + if (num_exo_times == 0) + mooseError("In SolutionUserObjectBase, exodus file contains no timesteps."); + + // Account for parallel mesh + if (dynamic_cast(_mesh.get())) + { + _mesh->allow_renumbering(true); + _mesh->prepare_for_use(); + } + else + { + _mesh->allow_renumbering(false); + _mesh->prepare_for_use(); + } + + // Create EquationSystems object for solution + _es = std::make_unique(*_mesh); + _es->add_system(_system_name); + _system = &_es->get_system(_system_name); + + // Get the variable name lists as set; these need to be sets to perform set_intersection + const std::vector & all_nodal(_exodusII_io->get_nodal_var_names()); + const std::vector & all_elemental(_exodusII_io->get_elem_var_names()); + const std::vector & all_scalar(_exodusII_io->get_global_var_names()); + + // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if + // provided + // This function could be called more than once, so clear the member variables so we don't keep + // adding to the vectors + _nodal_variables.clear(); + _elemental_variables.clear(); + _scalar_variables.clear(); + if (!_system_variables.empty()) + { + for (const auto & var_name : _system_variables) + { + if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end()) + _nodal_variables.push_back(var_name); + if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end()) + _elemental_variables.push_back(var_name); + if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end()) + // Check if the scalar matches any field variables, and ignore the var if it does. This + // means its a Postprocessor. + if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == + _nodal_variables.end() && + std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == + _elemental_variables.end()) + _scalar_variables.push_back(var_name); + } + } + else + { + _nodal_variables = all_nodal; + _elemental_variables = all_elemental; + + for (auto var_name : all_scalar) + // Check if the scalar matches any field variables, and ignore the var if it does. This means + // its a Postprocessor. + if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == + _nodal_variables.end() && + std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == + _elemental_variables.end()) + _scalar_variables.push_back(var_name); + } + + // Add the variables to the system + for (const auto & var_name : _nodal_variables) + _system->add_variable(var_name, FIRST); + + for (const auto & var_name : _elemental_variables) + _system->add_variable(var_name, CONSTANT, MONOMIAL); + + for (const auto & var_name : _scalar_variables) + _system->add_variable(var_name, FIRST, SCALAR); + + // Initialize the equations systems + _es->init(); + + // Interpolate times + if (_interpolate_times) + { + // Create a second equation system + _es2 = std::make_unique(*_mesh); + _es2->add_system(_system_name); + _system2 = &_es2->get_system(_system_name); + + // Add the variables to the system + for (const auto & var_name : _nodal_variables) + _system2->add_variable(var_name, FIRST); + + for (const auto & var_name : _elemental_variables) + _system2->add_variable(var_name, CONSTANT, MONOMIAL); + + for (const auto & var_name : _scalar_variables) + _system2->add_variable(var_name, FIRST, SCALAR); + + // Initialize + _es2->init(); + + // Update the times for interpolation (initially start at 0) + updateExodusBracketingTimeIndices(); + + // Copy the solutions from the first system + for (const auto & var_name : _nodal_variables) + { + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); + _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + } + + for (const auto & var_name : _elemental_variables) + { + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); + _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + } + + if (_scalar_variables.size() > 0) + { + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); + _exodusII_io->copy_scalar_solution( + *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); + } + + // Update the systems + _system->update(); + _es->update(); + _system2->update(); + _es2->update(); + } + + // Non-interpolated times + else + { + if (_exodus_time_index > num_exo_times) + mooseError("In SolutionUserObjectBase, timestep = ", + _exodus_time_index, + ", but there are only ", + num_exo_times, + " time steps."); + + // Copy the values from the ExodusII file + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index); + + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index); + + if (!_scalar_variables.empty()) + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_time_index); + + // Update the equations systems + _system->update(); + _es->update(); + } +} + +Real +SolutionUserObjectBase::directValue(const Node * node, const std::string & var_name) const +{ + // Get the libmesh variable and system numbers + unsigned int var_num = _system->variable_number(var_name); + unsigned int sys_num = _system->number(); + + // Get the node id and associated dof + dof_id_type node_id = node->id(); + const Node & sys_node = _system->get_mesh().node_ref(node_id); + mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0, + "Variable " << var_name << " has no DoFs on node " << sys_node.id()); + dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0); + + // Return the desired value for the dof + return directValue(dof_id); +} + +Real +SolutionUserObjectBase::directValue(const Elem * elem, const std::string & var_name) const +{ + // Get the libmesh variable and system numbers + unsigned int var_num = _system->variable_number(var_name); + unsigned int sys_num = _system->number(); + + // Get the element id and associated dof + dof_id_type elem_id = elem->id(); + const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id); + mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0, + "Variable " << var_name << " has no DoFs on element " << sys_elem.id()); + dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0); + + // Return the desired value + return directValue(dof_id); +} + +void +SolutionUserObjectBase::initialize() +{ +} + +void +SolutionUserObjectBase::finalize() +{ +} + +void +SolutionUserObjectBase::timestepSetup() +{ + // Update time interpolation for ExodusII solution + if (_file_type == 1 && _interpolate_times) + updateExodusTimeInterpolation(); +} + +void +SolutionUserObjectBase::execute() +{ +} + +void +SolutionUserObjectBase::initialSetup() +{ + + // Make sure this only happens once + if (_initialized) + return; + + // Create a libmesh::Mesh object for storing the loaded data. + // Several aspects of SolutionUserObjectBase won't work with a DistributedMesh: + // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel. + // .) We don't know if directValue will be used, which may request + // a value on a Node we don't have. + // We force the Mesh used here to be a ReplicatedMesh. + _mesh = std::make_unique(_communicator); + + // ExodusII mesh file supplied + if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) || + MooseUtils::hasExtension(_mesh_file, "exo", true)) + { + _file_type = "exodusII"; + readExodusII(); + } + + // XDA mesh file supplied + else if (MooseUtils::hasExtension(_mesh_file, "xda")) + { + _file_type = "xda"; + readXda(); + } + + // XDR mesh file supplied + else if (MooseUtils::hasExtension(_mesh_file, "xdr")) + { + _file_type = "xdr"; + readXda(); + } + + // Produce an error for an unknown file type + else + mooseError( + "In SolutionUserObjectBase, invalid file type: only .xda, .xdr, .exo and .e supported"); + + // Initialize the serial solution vector + _serialized_solution = NumericVector::build(_communicator); + _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL); + + // Pull down a full copy of this vector on every processor so we can get values in parallel + _system->solution->localize(*_serialized_solution); + + // Vector of variable numbers to apply the MeshFunction to + std::vector var_nums; + + // If no variables were given, use all of them + if (_system_variables.empty()) + { + _system->get_all_variable_numbers(var_nums); + for (const auto & var_num : var_nums) + _system_variables.push_back(_system->variable_name(var_num)); + } + + // Otherwise, gather the numbers for the variables given + else + { + for (const auto & var_name : _system_variables) + var_nums.push_back(_system->variable_number(var_name)); + } + + // Create the MeshFunction for working with the solution data + _mesh_function = std::make_unique( + *_es, *_serialized_solution, _system->get_dof_map(), var_nums); + _mesh_function->init(); + + // Tell the MeshFunctions that we might be querying them outside the + // mesh, so we can handle any errors at the MOOSE rather than at the + // libMesh level. + DenseVector default_values; + _mesh_function->enable_out_of_mesh_mode(default_values); + + // Build second MeshFunction for interpolation + if (_interpolate_times) + { + // Need to pull down a full copy of this vector on every processor so we can get values in + // parallel + _serialized_solution2 = NumericVector::build(_communicator); + _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL); + _system2->solution->localize(*_serialized_solution2); + + // Create the MeshFunction for the second copy of the data + _mesh_function2 = std::make_unique( + *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums); + _mesh_function2->init(); + _mesh_function2->enable_out_of_mesh_mode(default_values); + } + + // Populate the MeshFunction variable index + for (unsigned int i = 0; i < _system_variables.size(); ++i) + { + std::string name = _system_variables[i]; + _local_variable_index[name] = i; + } + + // Set initialization flag + _initialized = true; +} + +MooseEnum +SolutionUserObjectBase::getSolutionFileType() const +{ + return _file_type; +} + +void +SolutionUserObjectBase::updateExodusTimeInterpolation() +{ + if (_t != _interpolation_time) + { + if (updateExodusBracketingTimeIndices()) + { + + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); + if (_scalar_variables.size() > 0) + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); + + _system->update(); + _es->update(); + _system->solution->localize(*_serialized_solution); + + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + if (_scalar_variables.size() > 0) + _exodusII_io->copy_scalar_solution( + *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); + + _system2->update(); + _es2->update(); + _system2->solution->localize(*_serialized_solution2); + } + _interpolation_time = _t; + } +} + +bool +SolutionUserObjectBase::updateExodusBracketingTimeIndices() +{ + if (_file_type != 1) + mooseError("In SolutionUserObjectBase, getTimeInterpolationData only applicable for exodusII " + "file type"); + + int old_index1 = _exodus_index1; + int old_index2 = _exodus_index2; + + int num_exo_times = _exodus_times->size(); + + const auto solution_time = solutionSampleTime(); + + if (solution_time < (*_exodus_times)[0]) + { + _exodus_index1 = 0; + _exodus_index2 = 0; + _interpolation_factor = 0.0; + } + else + { + for (int i = 0; i < num_exo_times - 1; ++i) + { + if (solution_time <= (*_exodus_times)[i + 1]) + { + _exodus_index1 = i; + _exodus_index2 = i + 1; + _interpolation_factor = + (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); + break; + } + else if (i == num_exo_times - 2) + { + _exodus_index1 = num_exo_times - 1; + _exodus_index2 = num_exo_times - 1; + _interpolation_factor = 1.0; + break; + } + } + } + + bool indices_modified(false); + + if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2) + indices_modified = true; + + return indices_modified; +} + +unsigned int +SolutionUserObjectBase::getLocalVarIndex(const std::string & var_name) const +{ + // Extract the variable index for the MeshFunction(s) + std::map::const_iterator it = _local_variable_index.find(var_name); + if (it == _local_variable_index.end()) + mooseError("Value requested for nonexistent variable '", + var_name, + "' in the '", + name(), + "' SolutionUserObjectBase.\nSystem selected: ", + _system_name, + "\nAvailable variables:\n", + ConsoleUtils::formatString(Moose::stringify(_system_variables), "")); + return it->second; +} + +Real +SolutionUserObjectBase::pointValueWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type, + const std::set * subdomain_ids) const +{ + // first check if the FE type is continuous because in that case the value is + // unique and we can take a short cut, the default weighting_type found_first also + // shortcuts out + auto family = + _fe_problem + .getVariable( + _tid, var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD) + .feType() + .family; + + if (weighting_type == 1 || + (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC)) + return pointValue(t, p, var_name, subdomain_ids); + + // the shape function is discontinuous so we need to compute a suitable unique value + std::map values = discontinuousPointValue(t, p, var_name); + switch (weighting_type) + { + case 2: + { + Real average = 0.0; + for (auto & v : values) + average += v.second; + return average / Real(values.size()); + } + case 4: + { + Real smallest_elem_id_value = std::numeric_limits::max(); + dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); + for (auto & v : values) + if (v.first->id() < smallest_elem_id) + { + smallest_elem_id = v.first->id(); + smallest_elem_id_value = v.second; + } + return smallest_elem_id_value; + } + case 8: + { + Real largest_elem_id_value = std::numeric_limits::lowest(); + dof_id_type largest_elem_id = 0; + for (auto & v : values) + if (v.first->id() > largest_elem_id) + { + largest_elem_id = v.first->id(); + largest_elem_id_value = v.second; + } + return largest_elem_id_value; + } + } + + mooseError("SolutionUserObjectBase::pointValueWrapper reaches line that it should not be able to " + "reach."); + return 0.0; +} + +Real +SolutionUserObjectBase::pointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return pointValue(t, p, local_var_index, subdomain_ids); +} + +Real +SolutionUserObjectBase::pointValue(Real libmesh_dbg_var(t), + const Point & p, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // Create copy of point + Point pt(p); + + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids); + val = val + (val2 - val) * _interpolation_factor; + } + + return val; +} + +std::map +SolutionUserObjectBase::discontinuousPointValue( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return discontinuousPointValue(t, p, local_var_index, subdomain_ids); +} + +std::map +SolutionUserObjectBase::discontinuousPointValue( + Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + std::map map = + evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + std::map map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2); + + if (map.size() != map2.size()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size"); + + // construct the interpolated map + for (auto & k : map) + { + if (map2.find(k.first) == map2.end()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys"); + Real val = k.second; + Real val2 = map2[k.first]; + map[k.first] = val + (val2 - val) * _interpolation_factor; + } + } + + return map; +} + +RealGradient +SolutionUserObjectBase::pointValueGradientWrapper( + Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type, + const std::set * subdomain_ids) const +{ + // the default weighting_type found_first shortcuts out + if (weighting_type == 1) + return pointValueGradient(t, p, var_name, subdomain_ids); + + // the shape function is discontinuous so we need to compute a suitable unique value + std::map values = + discontinuousPointValueGradient(t, p, var_name, subdomain_ids); + switch (weighting_type) + { + case 2: + { + RealGradient average = RealGradient(0.0, 0.0, 0.0); + for (auto & v : values) + average += v.second; + return average / Real(values.size()); + } + case 4: + { + RealGradient smallest_elem_id_value; + dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); + for (auto & v : values) + if (v.first->id() < smallest_elem_id) + { + smallest_elem_id = v.first->id(); + smallest_elem_id_value = v.second; + } + return smallest_elem_id_value; + } + case 8: + { + RealGradient largest_elem_id_value; + dof_id_type largest_elem_id = 0; + for (auto & v : values) + if (v.first->id() > largest_elem_id) + { + largest_elem_id = v.first->id(); + largest_elem_id_value = v.second; + } + return largest_elem_id_value; + } + } + + mooseError("SolutionUserObjectBase::pointValueGradientWrapper reaches line that it should not be " + "able to reach."); + return RealGradient(0.0, 0.0, 0.0); +} + +RealGradient +SolutionUserObjectBase::pointValueGradient(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return pointValueGradient(t, p, local_var_index, subdomain_ids); +} + +RealGradient +SolutionUserObjectBase::pointValueGradient(Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids); + val = val + (val2 - val) * _interpolation_factor; + } + + return val; +} + +std::map +SolutionUserObjectBase::discontinuousPointValueGradient( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids); +} + +std::map +SolutionUserObjectBase::discontinuousPointValueGradient( + Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + std::map map = + evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + std::map map2 = + evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + if (map.size() != map2.size()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size"); + + // construct the interpolated map + for (auto & k : map) + { + if (map2.find(k.first) == map2.end()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys"); + RealGradient val = k.second; + RealGradient val2 = map2[k.first]; + map[k.first] = val + (val2 - val) * _interpolation_factor; + } + } + + return map; +} + +Real +SolutionUserObjectBase::directValue(dof_id_type dof_index) const +{ + Real val = (*_serialized_solution)(dof_index); + if (_file_type == 1 && _interpolate_times) + { + Real val2 = (*_serialized_solution2)(dof_index); + val = val + (val2 - val) * _interpolation_factor; + } + return val; +} + +Real +SolutionUserObjectBase::evalMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + DenseVector output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + (*_mesh_function)(p, 0.0, output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + (*_mesh_function2)(p, 0.0, output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + return output(local_var_index); +} + +std::map +SolutionUserObjectBase::evalMultiValuedMeshFunction( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::map> temporary_output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (temporary_output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + + // Fill the actual map that is returned + std::map output; + for (auto & k : temporary_output) + { + mooseAssert( + k.second.size() > local_var_index, + "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index " + << local_var_index << " does not exist"); + output[k.first] = k.second(local_var_index); + } + + return output; +} + +RealGradient +SolutionUserObjectBase::evalMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::vector output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->gradient(p, 0.0, output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->gradient(p, 0.0, output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + return output[local_var_index]; +} + +std::map +SolutionUserObjectBase::evalMultiValuedMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::map> temporary_output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (temporary_output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + + // Fill the actual map that is returned + std::map output; + for (auto & k : temporary_output) + { + mooseAssert( + k.second.size() > local_var_index, + "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index " + << local_var_index << " does not exist"); + output[k.first] = k.second[local_var_index]; + } + + return output; +} + +const std::vector & +SolutionUserObjectBase::variableNames() const +{ + return _system_variables; +} + +bool +SolutionUserObjectBase::isVariableNodal(const std::string & var_name) const +{ + return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) != + _nodal_variables.end(); +} + +Real +SolutionUserObjectBase::scalarValue(Real /*t*/, const std::string & var_name) const +{ + unsigned int var_num = _system->variable_number(var_name); + const DofMap & dof_map = _system->get_dof_map(); + std::vector dofs; + dof_map.SCALAR_dof_indices(dofs, var_num); + // We can handle only FIRST order scalar variables + return directValue(dofs[0]); +} + +void +SolutionUserObjectBase::readBlockIdMapFromExodusII() +{ +#ifdef LIBMESH_HAVE_EXODUS_API + libMesh::ExodusII_IO_Helper & exio_helper = _exodusII_io->get_exio_helper(); + const auto & id_to_block = exio_helper.id_to_block_names; + _block_name_to_id.clear(); + _block_id_to_name.clear(); + for (const auto & it : id_to_block) + { + _block_name_to_id[it.second] = it.first; + _block_id_to_name[it.first] = it.second; + } +#endif +} diff --git a/modules/optimization/include/userobjects/AdjointSolutionUserObject.h b/modules/optimization/include/userobjects/AdjointSolutionUserObject.h index 73eaa9e3ce07..e7bf89190842 100644 --- a/modules/optimization/include/userobjects/AdjointSolutionUserObject.h +++ b/modules/optimization/include/userobjects/AdjointSolutionUserObject.h @@ -9,19 +9,22 @@ #pragma once -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" -class AdjointSolutionUserObject : public SolutionUserObject +class AdjointSolutionUserObject : public SolutionUserObjectBase { public: static InputParameters validParams(); AdjointSolutionUserObject(const InputParameters & parameters); + virtual Real solutionSampleTime() override; + /** * Skipping parent class initialSetup since it will be called in timestepSetup */ virtual void initialSetup() override {} + /** * This will read a the files again if they have been re-written from optimization iteration */ @@ -30,8 +33,7 @@ class AdjointSolutionUserObject : public SolutionUserObject protected: /// Mapping between adjoint simulation time and adjoint simulation time const Real & _reverse_time_end; + /// The system time of the last instance the file was loaded std::time_t _file_mod_time; - /// The forward simulation time for last instance the solution was updated - Real _actual_interpolation_time; }; diff --git a/modules/optimization/src/userobjects/AdjointSolutionUserObject.C b/modules/optimization/src/userobjects/AdjointSolutionUserObject.C index 4998efa50286..63da95366324 100644 --- a/modules/optimization/src/userobjects/AdjointSolutionUserObject.C +++ b/modules/optimization/src/userobjects/AdjointSolutionUserObject.C @@ -9,6 +9,9 @@ #include "AdjointSolutionUserObject.h" +#include "libmesh/mesh_function.h" +#include "libmesh/exodusII_io.h" + #include registerMooseObject("OptimizationApp", AdjointSolutionUserObject); @@ -16,7 +19,7 @@ registerMooseObject("OptimizationApp", AdjointSolutionUserObject); InputParameters AdjointSolutionUserObject::validParams() { - InputParameters params = SolutionUserObject::validParams(); + InputParameters params = SolutionUserObjectBase::validParams(); params.addClassDescription( "Reads a variable from a mesh in one simulation to another specifically for loading forward " "solution in adjoint simulation during inverse optimization."); @@ -34,10 +37,9 @@ AdjointSolutionUserObject::validParams() } AdjointSolutionUserObject::AdjointSolutionUserObject(const InputParameters & parameters) - : SolutionUserObject(parameters), + : SolutionUserObjectBase(parameters), _reverse_time_end(getParam("reverse_time_end")), - _file_mod_time(std::numeric_limits::min()), - _actual_interpolation_time(0.0) + _file_mod_time(std::numeric_limits::min()) { if (!MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true)) paramError("mesh", @@ -45,6 +47,12 @@ AdjointSolutionUserObject::AdjointSolutionUserObject(const InputParameters & par "solution is written and read from exodus format."); } +Real +AdjointSolutionUserObject::solutionSampleTime() +{ + return _reverse_time_end - _t + _dt; +} + void AdjointSolutionUserObject::timestepSetup() { @@ -63,24 +71,11 @@ AdjointSolutionUserObject::timestepSetup() _es2.reset(); // Read the exodus file - SolutionUserObject::initialSetup(); + SolutionUserObjectBase::initialSetup(); // Make sure to communicate what solution was actually loaded - _actual_interpolation_time = 0.0; _interpolation_time = -_dt; } - // Update time interpolation for ExodusII solution - const Real actual_time = _reverse_time_end - _t + _dt; - if (_actual_interpolation_time != actual_time) - { - // This will make sure the solutions are updated for the current forward time - _interpolation_time = _actual_interpolation_time; - // Update based on the forward simulation time - updateExodusTimeInterpolation(actual_time); - // This avoids errors in public functions when the time is given as the adjoint simulation time - _interpolation_time = _t; - // Ensures we only update the solution when necessary - _actual_interpolation_time = actual_time; - } + SolutionUserObjectBase::timestepSetup(); } diff --git a/modules/phase_field/src/userobjects/SolutionRasterizer.C b/modules/phase_field/src/userobjects/SolutionRasterizer.C index ee27603eeb43..5942485ac4eb 100644 --- a/modules/phase_field/src/userobjects/SolutionRasterizer.C +++ b/modules/phase_field/src/userobjects/SolutionRasterizer.C @@ -10,6 +10,8 @@ #include "SolutionRasterizer.h" #include +#include "libmesh/mesh_function.h" +#include "libmesh/exodusII_io.h" registerMooseObject("PhaseFieldApp", SolutionRasterizer); diff --git a/test/include/postprocessors/TestDiscontinuousValuePP.h b/test/include/postprocessors/TestDiscontinuousValuePP.h index b7b4aae37f14..b8fb198c9db2 100644 --- a/test/include/postprocessors/TestDiscontinuousValuePP.h +++ b/test/include/postprocessors/TestDiscontinuousValuePP.h @@ -12,7 +12,7 @@ #include "GeneralPostprocessor.h" // Forward Declarations -class SolutionUserObject; +class SolutionUserObjectBase; /** * Compute the value of a variable or the gradient at a specified location. @@ -52,5 +52,5 @@ class TestDiscontinuousValuePP : public GeneralPostprocessor unsigned int _gradient_component; /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; }; diff --git a/test/src/postprocessors/TestDiscontinuousValuePP.C b/test/src/postprocessors/TestDiscontinuousValuePP.C index c67ca976e306..2effbf3b14a7 100644 --- a/test/src/postprocessors/TestDiscontinuousValuePP.C +++ b/test/src/postprocessors/TestDiscontinuousValuePP.C @@ -9,7 +9,7 @@ // MOOSE includes #include "TestDiscontinuousValuePP.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseTestApp", TestDiscontinuousValuePP); @@ -46,7 +46,7 @@ TestDiscontinuousValuePP::TestDiscontinuousValuePP(const InputParameters & param void TestDiscontinuousValuePP::initialSetup() { - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); } Real