From c309dace45a8c0592084fbec1410bf06d4a972ec Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 18 Oct 2023 10:26:06 -0700 Subject: [PATCH] Implement FaceArg evaluation for finite element variables Refs #19420 --- framework/include/variables/MooseVariableFE.h | 14 ++-- framework/src/variables/MooseVariableFE.C | 53 ++++++++++++-- .../fe-var-for-fv-neumann/gold/test_out.e | Bin 0 -> 43236 bytes .../functors/fe-var-for-fv-neumann/test.i | 68 ++++++++++++++++++ .../functors/fe-var-for-fv-neumann/tests | 18 +++++ 5 files changed, 143 insertions(+), 10 deletions(-) create mode 100644 test/tests/functors/fe-var-for-fv-neumann/gold/test_out.e create mode 100644 test/tests/functors/fe-var-for-fv-neumann/test.i create mode 100644 test/tests/functors/fe-var-for-fv-neumann/tests diff --git a/framework/include/variables/MooseVariableFE.h b/framework/include/variables/MooseVariableFE.h index a66a1ed5cdde..f8d232d3b150 100644 --- a/framework/include/variables/MooseVariableFE.h +++ b/framework/include/variables/MooseVariableFE.h @@ -698,16 +698,22 @@ class MooseVariableFE : public MooseVariableField using NodeArg = Moose::NodeArg; using ElemPointArg = Moose::ElemPointArg; + /** + * A common method that both evaluate(FaceArg) and evaluateDot(FaceArg) can call. A value + * evaluation vs dot evaluation is delineated via the passed-in \p cache_data, e.g. if the + * passed-in cache data is the sln data member then this will return a value evaluation and if the + * cache data is the dot data member then this will return a dot evaluation + */ + ValueType + faceEvaluate(const FaceArg &, const StateArg &, const std::vector & cache_data) const; + ValueType evaluate(const ElemQpArg & elem_qp, const StateArg & state) const override final; ValueType evaluate(const ElemSideQpArg & elem_side_qp, const StateArg & state) const override final; ValueType evaluate(const ElemArg &, const StateArg &) const override final; ValueType evaluate(const ElemPointArg &, const StateArg &) const override final; ValueType evaluate(const NodeArg & node_arg, const StateArg & state) const override final; - ValueType evaluate(const FaceArg &, const StateArg &) const override final - { - mooseError("Face info functor overload not yet implemented for finite element variables"); - } + ValueType evaluate(const FaceArg &, const StateArg &) const override final; GradientType evaluateGradient(const ElemQpArg & elem_qp, const StateArg & state) const override; GradientType evaluateGradient(const ElemSideQpArg & elem_side_qp, diff --git a/framework/src/variables/MooseVariableFE.C b/framework/src/variables/MooseVariableFE.C index 8180550a156a..47f53b7409e3 100644 --- a/framework/src/variables/MooseVariableFE.C +++ b/framework/src/variables/MooseVariableFE.C @@ -1031,6 +1031,52 @@ MooseVariableFE::evaluate(const ElemArg & elem_arg, const StateArg & return _current_elem_qp_functor_sln[0]; } +template +typename MooseVariableFE::ValueType +MooseVariableFE::faceEvaluate(const FaceArg & face_arg, + const StateArg & state, + const std::vector & cache_data) const +{ + const QMonomial qrule(face_arg.fi->elem().dim() - 1, CONSTANT); + auto side_evaluate = + [this, &qrule, &state, &cache_data](const Elem * const elem, const unsigned int side) + { + // We can use whatever we want for the point argument since it won't be used + const ElemSideQpArg elem_side_qp_arg{elem, side, /*qp=*/0, &qrule, Point(0, 0, 0)}; + evaluateOnElementSide(elem_side_qp_arg, state, /*cache_eligible=*/false); + return cache_data[0]; + }; + + const auto continuity = this->getContinuity(); + const bool on_elem = !face_arg.face_side || (face_arg.face_side == face_arg.fi->elemPtr()); + const bool on_neighbor = + !face_arg.face_side || (face_arg.face_side == face_arg.fi->neighborPtr()); + if (on_neighbor) + mooseAssert( + face_arg.fi->neighborPtr(), + "If we are signaling we should evaluate on the neighbor, we better have a neighbor"); + + // Only do multiple evaluations if we are not continuous and we are on an internal face + if ((continuity != C_ZERO && continuity != C_ONE) && on_elem && on_neighbor) + return (side_evaluate(face_arg.fi->elemPtr(), face_arg.fi->elemSideID()) + + side_evaluate(face_arg.fi->neighborPtr(), face_arg.fi->neighborSideID())) / + 2; + else if (on_elem) + return side_evaluate(face_arg.fi->elemPtr(), face_arg.fi->elemSideID()); + else if (on_neighbor) + return side_evaluate(face_arg.fi->neighborPtr(), face_arg.fi->neighborSideID()); + else + mooseError( + "Attempted to evaluate a moose finite element variable on a face where it is not defined"); +} + +template +typename MooseVariableFE::ValueType +MooseVariableFE::evaluate(const FaceArg & face_arg, const StateArg & state) const +{ + return faceEvaluate(face_arg, state, _current_elem_side_qp_functor_sln); +} + template typename MooseVariableFE::ValueType MooseVariableFE::evaluate(const ElemPointArg & elem_point_arg, @@ -1214,12 +1260,7 @@ MooseVariableFE::evaluateDot(const FaceArg & face_arg, const StateAr mooseAssert(_time_integrator && _time_integrator->dt(), "A time derivative is being requested but we do not have a time integrator so we'll " "have no idea how to compute it"); - const QMonomial qrule(face_arg.fi->elem().dim() - 1, CONSTANT); - // We can use whatever we want for the point argument since it won't be used - const ElemSideQpArg elem_side_qp_arg{ - face_arg.fi->elemPtr(), face_arg.fi->elemSideID(), /*qp=*/0, &qrule, Point(0, 0, 0)}; - evaluateOnElementSide(elem_side_qp_arg, state, /*cache_eligible=*/false); - return _current_elem_side_qp_functor_dot[0]; + return faceEvaluate(face_arg, state, _current_elem_side_qp_functor_dot); } template <> diff --git a/test/tests/functors/fe-var-for-fv-neumann/gold/test_out.e b/test/tests/functors/fe-var-for-fv-neumann/gold/test_out.e new file mode 100644 index 0000000000000000000000000000000000000000..588790717cff3340929068c4783e8d07b837b44e GIT binary patch literal 43236 zcmeHQ&5z{9bsxpHH2Y;mTFLOi2?*%m1i|jk$I81aq(P5%XE&SdSG22D93>&xWU-rC zHrXQCvom`N{15mPAP#~+C&5NeK~6pd0b;;`6B$-u$AJ+ckt0R6BSrBchuoast77r1 zNlsTyYdxds!t^xBF8=COz500d>eZ`@moER(($W%rj^gtqKFa|Mtp2Fv-xVQSjJsyd`PssPI zgr(WLL<@hO_%70|JpXit;fDmX%5ez$@eIcrvlwae9|-3-Jr_6$J)|_t&(^YM$K?AV zkH+Jbo>8>wvn@N8&p00C(-IFU4W4&D^dh|W_jFHX8u%f~bm@2TIgZa`_#CsxzO}{T z#E(MUE`180M>Zb$yRYE;v6knb^@7Mw{cvE7BR@>7`%&y7j}Hj@lk`45%lI7e`~c-A zeklD$zn>i0+y01n20mLLZjcV~IDg;3^>KiI;Vc+$rxRk0#Kz^)r&2!+C_w}&=VK{L zQY)IIr^Ul#c=jv!9OZTA{xHfSQ`+>A;Za-=tq8AQ!oMS|UnGcoic4i81JOVqgylEr zAL5VWXe@b4aj9(Sy-$id%<22YRYDi<%J<*I{gE)4N{?ugf#PQLEq|Y`5x<-s#oaAE z`TkpodnB1okLrO8<@Dag!)3%Tr$=#TP45H5ebk9iJGMDIl?k;~8Ho0;iJRsB#kYvN z!=H4-J@V&y+y}eGJ$eZ-h)2WG>6H#|44y8{4*Yffx<%i0@ zj3s$BBNjnZIL5J1wq-2I2SuzCgmsk4IY$%ul!?uVCt342;T@;G)S^ylWK*|W=Lgl3 zJNO_xN%5DRC=8jCQahk@sm>k2N2Yy*_fP50SQIgj^XTTKE1SG$6Yb(P#78*!Mf^L? z<5$zCG-aUcTjJ*A+xTAQ@sqrdrLjZ?ikrRr)C#VbYyELs+_Mxn!_E6?37TO_-18;e z-QxZx#TDh7(@W(fLz&($;-QQyaZkvI_%6r2g!|)Ep=cn?#scESr?dEy7cg3}@~S6RPC57jvt zh@RKP&Fj>+i24Lltkmw~p7KEbfTRbP<+A(?rHApyi>#PA=)q*)@ymGr1VzpQaK7?? z;`%d`Z`oPUcm2A2&tEaW{#6QpW_~04{^mLl=C`x^HwCi!z3ll9JRZy+X1~An28HjK zZ)fp-ls$hZi~nx+`+M1YKgr_#G>i9s7Vm>B-iKMdpNV*j|BDisX9)z%OXrDq(6C~E zFg(l!u;BmU=Kn!!uM2t?B{1g_AQ_NkPnv9z>`pQ$$*MH>Bl(r)Of(-M{RPR$B zPI9W0J)gw&NqnBdhxA+$Q{qZ&i7(TT>BzKXdNNI!u1s5|FKLi;NLnO4Bo7l!k}gSG zojysIq({;r)0gSW^kg~`U*b}J=|lLGXS#mcWB7lFgST~vGcU6K$TJOLPo?sVaWQ6g z%3PzrV_d&^^Y*6k%-#U!PT@14F&d5R&^4|>y=!d5gUN`6X<{HL4O)%Ur_YGLQ4}T2 zBE9M-lNgonYkvPaONKH_1*Bu2cf9lG zR{Lk4b=cYU^Y*%Zer0{l^VZL;tUAnFIqTSV|NMmt8g6cF-ns2e!x_f<>D3FT*Xz-D zzrTHL^Hw`xrUuy^rigKI633_{#vQDz6!}F_zJ$7d(@Blh^Tz7RbF1smon7O_SG?MM zwYkll2~;ZqvNdeR0{JU0wKB-*xE0~kH(8NU`rmdiYQ5~ROgyx4X{fPU)Rf#t=~$l^;{w5CMM zkvxf>ac$$y=IuLe0QKihfY3`8y@MZ-Y{B#B&0eZUj}JMRR^Ty9Z62Zch9)8@#5{J5^-_!+81h8PdDWi%jBc=ODm`fafo6@BzYhU90b)lsFA#!t!eZ>>LEt zo8B&ifO_Wuz~pBC6|8THY!S_+px%8MMrT{;8JGOncZLB=pTC%wRQm=beQfex41b{gAQLQi9VTk;}D*P@UER++meKz}q2d6-kC!xcXgkjdy7JA0YuGuU{Ae;#% z-m-GL7KF`V30((}_X#D8%B3D)?6>!+pbq%%#16_O)qYK*-sv;* z8uHrH(f@^y%X~d>PE=As>5R**Ymx-G7?-IOif|CdE=KABLDdZbR?9&D?Cig`cvRE7*ZL^W35iIQAVI+cvbYIL%qLvE3~@behzD z>OD5qLs}nnti-XgS*2VFwT1w!24a1N^bJBmE@u%|pLtQtsLo?*2@O>;Q&EyYEC4h5 zz@19g^;!znEBMwNY&EiDXJ~O{Ik+9zqrPifSXU$R)>>*x<*^Xi>wrNT1zCGmXG*3R zhkD5>V2a)WS!So6ow*t=fkL6WAB*La0%@-l#tf@Sc}5FdPT>j}=+72NAVy!9TH~?M zr4>-(N}<*96bUe)$#SI?=XEiR;*sSHnR9QZN(jgTFl)toz{)7>)O_jmhvg2dq3)I}Yax z*w9?mTZ%tGnUTs_NS%p=g0!d`mmRoTe`6Sd)lgZ=TvqLE0zQ1!O0EL9f-q~XSt1IX zI++|%l;N`~fH2$1y8axcbkQ@i#p)R$L16)Cjcw;17T`;GsE1>0X?13Vqn1LIdE3MO z(#hC`HB=dlu6xGWwF_r2eBs>Mg_-gVK+kinc8RHH+#kMqUq9DQwHG;4-qjs?t(G zA@5g7%FMe}f)YV1%Ul+ogb>)K1Lv{=pHXePdxxH(P=KNy2G}na@V)AD<8nhso#&l% zSxL5`EtP}Vcbf_jumDtPc}p-$@I0(l*n=pfoq`yp_-s=4{{Wo5%23JxoTL!iTj(PS zRV{*ecV|Uon%OqM?hz1KhJj8f_VXsU4@nY_4x7m3Bg9sA`NCGg=!h z>#9Ik5p=RV4ah~tOp>6X3zP4e1RLIcN!7q0jwUpL%p_w;qk3CT?)hVjZId>LOfJPt zQNmFfs&=qR$dZwVfr@h8#1z%02g4H>b+pQJJ>$maOLuPGm~~(=448HJDR~$vvRS@T zFGI36C?nLa)yoiO9btmraKyE8NeUks4oi|KvmI2g*k!<}N~6x1 z12Vr+&oaEh^sytLR&@ETD&PT7-Cw>tU;DmURQIF@qI8AIo>S1?ZgyFr*Qqvh*0i_N zuKCL^$e^*tl6G@Ol3;gXi+cC%KQ=T-YBHr)1Eo*bmwMP?#Ui z+WAn<4!$v-($-MQa*G`c{X}C5^lT(f%vgvEV>}Yut4acHBwrs>IxM%S+&dT5wVGbn zPAPQ>_Td}JU%q=8szA)lZ?Fk=z=qwn@>GTodrawIX+!HF3Z>bv*uzUtNYBV1b+YVL!3c3|%`?R%fm>@aVE1b7H#oBhAhfG=`!Z|xzuUg*n(3SyYX zTEOlb4HnW@-Y6c<#moC#_99^kz6c(^^i`I<$REKyE~7TeVAL2E(ifFN5opij?6Q{- ziTA(swb=8x*z;)Y1FguR^hMCEx0!Wafe)+HPwns3iUXjxA@7It(ufdOcF~EDz$vyU zUdkM}xVPd6`(s8zK|GZ)bw{F5)lAC|onV4fU7*Q; zJ3CnF@~K~LWu-oO6qdR*n6l)cYMI8}xmjgJNP&hcb*2cw>e^cgQUj2D;=w6W3z6z- zmV74POr$y>mY8|uymwZ{ zz(SIYLUWY74NK0<3M9MP%;b{W=}L^sw|YHnH3t>AQ0l^K>U38xDXM>?KMZ{`xRhm8 zArCj|Stgwbw=}IoNeX9gC4$Q|yG$m17n_)Aw_bgZg?euz+-Aa4?_|_xaRZlChX6Ym z4{a-saHO8p8J07zPGO&=9e7}h>M+WC>@rQ&YYo7OAGoxqI_nR~wQ0$xmU^{u;MFsR zw!Ub7Xux92FqNq~g+-5@chRE9&O2~`&dzBPmf+P>_(rXkp19MYkN{g-bkF#`Fa3VK z!FJq}l|yOk?lg@m*%eK%b0t7$2ta=fA!)$iw!@y<`$2%J6? zlK9h7{0 z!nzO6{#zY(Pqk{ph#y#@!@-HJ4Xo)IUz&QU>p=;y*}JoPJ)! zpn)~&fLxptLZEejh@*O|vWa?sh^C(Ha)*5*M} zr$7A2Ogp4kelkMAxfHVJ)+xh*J;4FaLDaWX81T@Z`l8mL9x4sXOY~}|bfX2>*6b)X zn06?U$SG`Vb^>trIf`t0VOz7)1r#btdb=U$lp$ZBxes4MH4AQ-(?CqOxc+?g*)@6@WH$(d@x@M=jWlI}#a_O=ty3-l$XPL_yt2js&A`^=Snu#v>f9 z*AlbBFd~Z;gp@K8&SuOzK(R%s}PIWHV z0I_Rm1Y|nmmj?hUID8~^UmcKhyFi4?$HoD@?I9gcg{sZ{cf;=}Sc1iq0gPL4ZZxl; z8d)v_N{?(glQo7>K_OPohHxHcheTaj|DeHTPGR9~bt>0tnU-Z(0M>G?2FJO~DdfBk zqG>N;I$4eepsHX)(Msq_9gx&ORgB7=dRyksa>$}+2otYp>exOF3U@+BxDBgNCU4Yr zUnz6uJL+wvtvfP_O7$#nySB_;)3**_**oDBrUu}sZ3N}us(=STZ*y0AP@LOrudOfF ze9mX|9-AL|#L3S2DEA&C%SB8^8NNBR^P_ zue`8vQFwoh;jNY`4{G zjWGqQxAS%#K&Q~)JgKnJNw!Ws#Fk*zJe>!}RL2zJa;G_-4uxLV3|xmoxPHg>FQT;} zzz8Z(SsWz@F*$2N)TxKuVGfQwZJ=Yu^l`(;L`{5sP88s&&z+ z*o$Gb+N5j4=;gaFu+D>y>FvcZS}#L&3D&W(1Y$W1TLYj2bN$MIsxwF(#&OxkxmjGd zJMGa2l7kIx7>*yBaXcWF?lRD;8=vg|fj+$-S=a{3jmqbLRK?gY)LU@OUqo z|H|LnGXIUo{}c0V{`()8|H1D+G~eOzzHI&{fB$Lo$2|Vm%%AXhrukny-anh~^LU~8 O-#p%<=Kt_`fARkp`|F7S literal 0 HcmV?d00001 diff --git a/test/tests/functors/fe-var-for-fv-neumann/test.i b/test/tests/functors/fe-var-for-fv-neumann/test.i new file mode 100644 index 000000000000..f82bf2bc18d7 --- /dev/null +++ b/test/tests/functors/fe-var-for-fv-neumann/test.i @@ -0,0 +1,68 @@ +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 1 + nx = 20 + [] +[] + +[Variables] + [fe][] + [fv] + type = MooseVariableFVReal + [] +[] + +[Kernels] + [diff] + type = Diffusion + variable = fe + [] +[] + +[FVKernels] + [diff] + type = FVDiffusion + variable = fv + coeff = 1 + [] +[] + +[BCs] + [left] + type = DirichletBC + variable = fe + value = 0 + boundary = left + [] + [right] + type = DirichletBC + variable = fe + value = 1 + boundary = right + [] + [] + +[FVBCs] + [left] + type = FVDirichletBC + variable = fv + value = 0 + boundary = left + [] + [right] + type = FVFunctorNeumannBC + variable = fv + functor = fe + boundary = right + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/functors/fe-var-for-fv-neumann/tests b/test/tests/functors/fe-var-for-fv-neumann/tests new file mode 100644 index 000000000000..1280dc04720a --- /dev/null +++ b/test/tests/functors/fe-var-for-fv-neumann/tests @@ -0,0 +1,18 @@ +[Tests] + design = 'Functors/index.md' + issues = '#19420' + [exo] + type = Exodiff + input = test.i + exodiff = test_out.e + requirement = 'The system shall be able to accurately evaluate a finite element variable through the functor system in a finite volume Dirichlet boundary condition.' + [] + [jac] + type = PetscJacobianTester + ratio_tol = 1e-7 + difference_tol = 1e-5 + input = test.i + run_sim = True + requirement = 'The system shall compute the correct Jacobian when evaluating a finite element variable through the functor system in a finite volume Dirichlet boundary condition.' + [] +[]