From 88ff6acedd89ed984f26e3ba560186cfd0c1f44b Mon Sep 17 00:00:00 2001 From: Peter German Date: Wed, 27 Sep 2023 10:18:58 -0600 Subject: [PATCH 1/8] Add modifications to be able to add FV interface kernel to multiple objects. (#25599) --- framework/include/fviks/FVInterfaceKernel.h | 36 ++++-- framework/src/fviks/FVInterfaceKernel.C | 80 ++++++++++--- framework/src/problems/FEProblemBase.C | 33 ++++- test/include/executioners/SteadySolve2.h | 1 + test/src/executioners/SteadySolve2.C | 14 ++- .../fviks/diffusion/gold/multisystem_out.e | Bin 0 -> 58616 bytes test/tests/fviks/diffusion/multisystem.i | 113 ++++++++++++++++++ test/tests/fviks/diffusion/tests | 8 ++ 8 files changed, 255 insertions(+), 30 deletions(-) create mode 100644 test/tests/fviks/diffusion/gold/multisystem_out.e create mode 100644 test/tests/fviks/diffusion/multisystem.i diff --git a/framework/include/fviks/FVInterfaceKernel.h b/framework/include/fviks/FVInterfaceKernel.h index 2e32073fef06..c632cf95a882 100644 --- a/framework/include/fviks/FVInterfaceKernel.h +++ b/framework/include/fviks/FVInterfaceKernel.h @@ -100,10 +100,14 @@ class FVInterfaceKernel : public MooseObject, const std::set & sub2() const { return _subdomain2; } /** - * @return The system associated with this object. Either an undisplaced or displaced nonlinear - * system + * @return The system associated with the first variable */ - const SystemBase & sys() const { return _sys; } + const SystemBase & sys() const { return _sys1; } + + /** + * @return The system associated with the second variable + */ + const SystemBase & sys2() const { return _sys2; } /** * @return Whether the \p FaceInfo element is on the 1st side of the interface @@ -133,13 +137,16 @@ class FVInterfaceKernel : public MooseObject, /** * Process the provided residual given \p var_num and whether this is on the neighbor side */ - void addResidual(Real resid, unsigned int var_num, bool neighbor); + void addResidual(Real resid, unsigned int var_num, Assembly & assembly, bool neighbor); using TaggingInterface::addJacobian; /** * Process the derivatives for the provided residual and dof index */ - void addJacobian(const ADReal & resid, dof_id_type dof_index, Real scaling_factor); + void addJacobian(const ADReal & resid, + dof_id_type dof_index, + Assembly & assembly, + Real scaling_factor); /** * @return A structure that contains information about the face info element and skewness @@ -192,16 +199,25 @@ class FVInterfaceKernel : public MooseObject, /// The SubProblem SubProblem & _subproblem; - /// the system object - SystemBase & _sys; + /// the system object for variable 1 + SystemBase & _sys1; - /// The Assembly object - Assembly & _assembly; + /// the system object for variable 2 + SystemBase & _sys2; -private: + /// Variable on one side of the interface MooseVariableFV & _var1; + + /// Variable on the other side of the interface MooseVariableFV & _var2; + /// The Assembly object for system 1 + Assembly & _assembly1; + + /// The Assembly object for system 2 + Assembly & _assembly2; + +private: std::set _subdomain1; std::set _subdomain2; diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index 36bda2e08522..16c661134de4 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -98,13 +98,19 @@ FVInterfaceKernel::FVInterfaceKernel(const InputParameters & parameters) ADFunctorInterface(this), _tid(getParam("_tid")), _subproblem(*getCheckedPointerParam("_subproblem")), - _sys(*getCheckedPointerParam("_sys")), - _assembly(_subproblem.assembly(_tid, _sys.number())), - _var1(_sys.getFVVariable(_tid, getParam("variable1"))), - _var2(_sys.getFVVariable(_tid, - isParamValid("variable2") - ? getParam("variable2") - : getParam("variable1"))), + _sys1(_subproblem.getVariable(_tid, getParam("variable1")).sys()), + _sys2(_subproblem + .getVariable(_tid, + isParamValid("variable2") ? getParam("variable2") + : getParam("variable1")) + .sys()), + _var1(_sys1.getFVVariable(_tid, getParam("variable1"))), + _var2(_sys2.getFVVariable(_tid, + isParamValid("variable2") + ? getParam("variable2") + : getParam("variable1"))), + _assembly1(_subproblem.assembly(_tid, _sys1.number())), + _assembly2(_subproblem.assembly(_tid, _sys2.number())), _mesh(_subproblem.mesh()) { if (getParam("use_displaced_mesh")) @@ -157,9 +163,12 @@ FVInterfaceKernel::setupData(const FaceInfo & fi) } void -FVInterfaceKernel::addResidual(const Real resid, const unsigned int var_num, const bool neighbor) +FVInterfaceKernel::addResidual(const Real resid, + const unsigned int var_num, + Assembly & assembly, + const bool neighbor) { - neighbor ? prepareVectorTagNeighbor(_assembly, var_num) : prepareVectorTag(_assembly, var_num); + neighbor ? prepareVectorTagNeighbor(assembly, var_num) : prepareVectorTag(assembly, var_num); _local_re(0) = resid; accumulateTaggedLocalResidual(); } @@ -167,9 +176,10 @@ FVInterfaceKernel::addResidual(const Real resid, const unsigned int var_num, con void FVInterfaceKernel::addJacobian(const ADReal & resid, const dof_id_type dof_index, + Assembly & assembly, const Real scaling_factor) { - addJacobian(_assembly, + addJacobian(assembly, std::array{{resid}}, std::array{{dof_index}}, scaling_factor); @@ -183,10 +193,29 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi) const auto var_elem_num = _elem_is_one ? _var1.number() : _var2.number(); const auto var_neigh_num = _elem_is_one ? _var2.number() : _var1.number(); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual()); - addResidual(r, var_elem_num, false); - addResidual(-r, var_neigh_num, true); + bool is_var1_system = + (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()); + + // std::cout << "Adding this " + + // if (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()) + // { + // addResidual(_elem_is_one ? r : -r, _var1.number(), _assembly1, _elem_is_one ? true : false); + // addResidual(-r, var_neigh_num, assembly_neigh, true); + // } + // else + // addResidual(-r, var_neigh_num, assembly_neigh, true); + + // if (&_sys1 == &_sys2) + // addResidual(r, var_neigh_num, assembly_neigh, true); + + addResidual(r, var_elem_num, assembly_elem, false); + addResidual(-r, var_neigh_num, assembly_neigh, true); } void @@ -208,12 +237,35 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi) const auto elem_scaling_factor = _elem_is_one ? _var1.scalingFactor() : _var2.scalingFactor(); const auto neigh_scaling_factor = _elem_is_one ? _var2.scalingFactor() : _var1.scalingFactor(); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); + bool is_var1_system = + (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()); + + // if (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()) + // { + // addResidualsAndJacobian(_assembly1, + // std::array{{_elem_is_one ? r : -r}}, + // _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), + // _var1.scalingFactor()); + // addResidualsAndJacobian( + // assembly_neigh, std::array{{-r}}, neigh_dof_indices, neigh_scaling_factor); + // } + // else + // addResidualsAndJacobian( + // assembly_neigh, std::array{{r}}, neigh_dof_indices, neigh_scaling_factor); + + // if (&_sys1 == &_sys2) + // addResidualsAndJacobian( + // assembly_neigh, std::array{{r}}, neigh_dof_indices, neigh_scaling_factor); + addResidualsAndJacobian( - _assembly, std::array{{r}}, elem_dof_indices, elem_scaling_factor); + _assembly1, std::array{{r}}, elem_dof_indices, elem_scaling_factor); addResidualsAndJacobian( - _assembly, std::array{{-r}}, neigh_dof_indices, neigh_scaling_factor); + _assembly1, std::array{{-r}}, neigh_dof_indices, neigh_scaling_factor); } Moose::ElemArg diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index 7bd8ff2f71e5..f51cc8587fcc 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -176,7 +176,7 @@ FEProblemBase::validParams() "True to skip the NonlinearSystem check for work to do (e.g. Make sure " "that there are variables to solve for)."); params.addParam("allow_initial_conditions_with_restart", - false, + true, "True to allow the user to specify initial conditions when restarting. " "Initial conditions can override any restarted field"); @@ -2929,7 +2929,25 @@ FEProblemBase::addFVInterfaceKernel(const std::string & fv_ik_name, const std::string & name, InputParameters & parameters) { + const auto nl_sys_num1 = + determineNonlinearSystem(parameters.varName("variable1", name), true).first + ? determineNonlinearSystem(parameters.varName("variable1", name), true).second + : (unsigned int)0; + + const auto nl_sys_num2 = + determineNonlinearSystem(parameters.varName("variable2", name), true).first + ? determineNonlinearSystem(parameters.varName("variable2", name), true).second + : (unsigned int)0; + + parameters.set("_sys") = _nl[nl_sys_num1].get(); addObject(fv_ik_name, name, parameters); + + // if (nl_sys_num1 != nl_sys_num2) + // { + // InputParameters params2 = parameters; + // params2.set("_sys") = _nl[nl_sys_num2].get(); + // addObject(fv_ik_name, name + "_var2", params2); + // } } // InterfaceKernels //// @@ -3558,6 +3576,7 @@ FEProblemBase::swapBackMaterialsNeighbor(THREAD_ID tid) void FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::string & object_name) { + const bool system_already_set = parameters.get("_sys"); const auto nl_sys_num = parameters.isParamValid("variable") && determineNonlinearSystem(parameters.varName("variable", object_name)).first @@ -3568,7 +3587,8 @@ FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::st parameters.get("use_displaced_mesh")) { parameters.set("_subproblem") = _displaced_problem.get(); - parameters.set("_sys") = &_displaced_problem->nlSys(nl_sys_num); + if (!system_already_set) + parameters.set("_sys") = &_displaced_problem->nlSys(nl_sys_num); } else { @@ -3580,7 +3600,8 @@ FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::st parameters.set("use_displaced_mesh") = false; parameters.set("_subproblem") = this; - parameters.set("_sys") = _nl[nl_sys_num].get(); + if (!system_already_set) + parameters.set("_sys") = _nl[nl_sys_num].get(); } } @@ -6019,12 +6040,16 @@ FEProblemBase::computeResidualAndJacobian(const NumericVector & soln, { // vector tags { + _current_nl_sys->associateVectorToTag(residual, _current_nl_sys->residualVectorTag()); const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL); _fe_vector_tags.clear(); for (const auto & residual_vector_tag : residual_vector_tags) - _fe_vector_tags.insert(residual_vector_tag._id); + // We filter out tags which do not have associated vectors in the current nonlinear + // system. This is essential to be able to use system-dependent residual tags. + if (_current_nl_sys->hasVector(residual_vector_tag._id)) + _fe_vector_tags.insert(residual_vector_tag._id); setCurrentResidualVectorTags(_fe_vector_tags); } diff --git a/test/include/executioners/SteadySolve2.h b/test/include/executioners/SteadySolve2.h index 98731794e7bd..3403d1a8d759 100644 --- a/test/include/executioners/SteadySolve2.h +++ b/test/include/executioners/SteadySolve2.h @@ -42,4 +42,5 @@ class SteadySolve2 : public Executioner const unsigned int _first_nl_sys; const unsigned int _second_nl_sys; bool _last_solve_converged; + unsigned int _number_of_iterations; }; diff --git a/test/src/executioners/SteadySolve2.C b/test/src/executioners/SteadySolve2.C index 8afd511da4b6..8c98875c1d79 100644 --- a/test/src/executioners/SteadySolve2.C +++ b/test/src/executioners/SteadySolve2.C @@ -27,6 +27,10 @@ SteadySolve2::validParams() "The first nonlinear system to solve"); params.addRequiredParam("second_nl_sys_to_solve", "The second nonlinear system to solve"); + params.addRangeCheckedParam("number_of_iterations", + 1, + "number_of_iterations>0", + "The number of iterations between the two systems."); return params; } @@ -37,7 +41,8 @@ SteadySolve2::SteadySolve2(const InputParameters & parameters) _time_step(_problem.timeStep()), _time(_problem.time()), _first_nl_sys(_problem.nlSysNum(getParam("first_nl_sys_to_solve"))), - _second_nl_sys(_problem.nlSysNum(getParam("second_nl_sys_to_solve"))) + _second_nl_sys(_problem.nlSysNum(getParam("second_nl_sys_to_solve"))), + _number_of_iterations(getParam("number_of_iterations")) { _fixed_point_solve->setInnerSolve(_feproblem_solve); @@ -107,7 +112,12 @@ SteadySolve2::execute() return _problem.nlConverged(sys_num); }; - const bool converged = solve(_first_nl_sys) && solve(_second_nl_sys); + bool converged = true; + + for (unsigned int i = 0; i < _number_of_iterations; i++) + { + converged = converged && solve(_first_nl_sys) && solve(_second_nl_sys); + } if (!converged) { diff --git a/test/tests/fviks/diffusion/gold/multisystem_out.e b/test/tests/fviks/diffusion/gold/multisystem_out.e new file mode 100644 index 0000000000000000000000000000000000000000..e9f1645a2fd77fe78164a9ba29f545bd97b849ff GIT binary patch literal 58616 zcmeHQ3A7wlnSLQ`2|-kXh-{?<5|Edj1kkp4NnXN}ERYu>dr{q8-LI1Fs!nzF%S+@i zaE{A};~>Zx1w4v?j@xitILe?Xf}`lL#Bl>1G3>IMJrI&G-+z~?>h9`(eJd^GP5Pem zzpkoR_0_-Kd;h!Cy>sWy-><8yi)bw92++X=%P}0YXmxds!fDqS-!584*|$n`#?wJg zrD$aBBAtx}jiGbL&01yte>~fDFkQ0>7TuHujpgS?y0DVnMDILQK4+9IUpz~A?0VVG z;=NVZ-zG3@hhU~Loay*`n7}bgmWOw;e-M84{d8wnSIqtFA9VjX!r|{N69_3FdW~+Z zy*5K!J4Eu}ut(W_k!$Q7^laoBJsRSXu8cv`!!_cC>Cbj@E>3&toa|Dt9n$T2pmCss zK!=zm+ZeRGvh6xJ?K%K7s%O-F*W&k}2G`%6E4Zd_J9(qz+Kz9mcD*d*UQ5{f)AOK% zL1S`u0XE?Oi2RTKe|XUxvWsMQ_=*AM^8NrW2hVFb9tZGQGvPd@Um^2KY#eqS;M;zI zUWi2UzOz^<_;z^>JgsQB6@Rk$-y!&)1{%x!aP@#I9wHuyWI7NB#52<4U&o&@Ryuq! z&gov#jZDOgxPQzobd38;o?DW3r+Z1q^xXbDKQKN=_szo)%KIAXD^x58&d#xIg+ny0_K${SNmYm~r7HLkymD zLitoC@IHDHKN+(gzx>`c>`Ie+4<-Had)E$k@7t&14zdxMf_#a5M%q)#G1i!G7YfVG zLWMgI-9zuEgkw6Tbod4ML;8+5>(v>IPg#I*!SbaUI9X_>>Rmy=0>?Ao=bwteo=Kl(>r7JqPyA zThPl&IF>I)eHQm49sd&k94zqY--Ji?4CUzrc<;S}o@4K?mWO;@CXz?I{{(u+!R-EW zSf0xp@k`H%XD57-?yW8R8{B&f-79cOF4F0tAeo}Lx8tJ39q!(*)4c*$)+=w5p!kjA zemjCY+`ZqSd$}J3dN`~B?wyf)CBJWPaPM7_dnFI=XmBr;S%*?4VL^eY=x~jE;+^6d z=@TvH;|cfTxpxejg>ry7%}SHwA`vdxh5T10RU<3Tb3=ZBB-lator zmt-KKB#-5&1Kc?Pu#xi1$57 z@b_J-KPBG#6n~Naw0QqBkmkjYAFn?n_;@x)DX%{#_<4@rzTgi1Z-TGCk$97~>CX%P zUSMp0NPj`_$zJ>7$Mu&4zb`#V=^_1P!S^c+?g;%A!T;YGdza`Ngd7{hdo~DpHi-9a zWNJ8D-zeU@QOLJRynmCBbCcj>vygW)nG5`E7IJS9d~Ff(ZxQ@$6?)hz_++o$dY`_H z@jG7M#`KcYw+X(t3;k>t{BIX}dR55rs*vMVA;+sij@N`7uL(I`6LRbja_nGod{*DV zuI&JXO<3K+q*YX|Ec|KGR$ zKW+J|EsKpuw8Nb`7#$6CG%&I>K>d)^|40KR)Nf6Fqcp}sV*OU7#aD?*<(OdJpJm&@rH6LC1lP2Wg=9g5C!@0W<-0BIqR0M9?JA z$)Hm}?+2X<>IO{)O#w{>Q8zMm2UEZHOwehdS)kdV(?Mr|&IEk`^g&Px^dXQAIt$bT zngf~(ng{9y%?B+2oek;(eHe5OXd!43XfbFB=v>fwprxRG&@#~Zpyi;CfGz-C2)YP# zG3XM|rJ&0|2Iz8-2}*-9pe)D&jRNf_Xf%F_$Tmb|A0o1mw3D=zv=@y7kZO7wTvahtUw6nCew70amw7ayuw7>KL=?BsmG+>iI5WP>( z3HT+NAm~K=5|NJ(k)KFkk^UlmM*5BP9qB*Phom1#Uy}YLeMGRU>rSD7smw7ATWrrLSsuhxAYBm(m}l zA4>m|ekc7+`kC}E=~vR9q#sHDk$xlnMfwT(i1ZEV6Vex?4@lcfn@d|u8%x_tn@U?s z8%oR zo;P{&2c$qDWc{tJT{*XRS%0P;PSa*go;GXptkyuHcl9sqJ-0C+6|{$ghUv{}&FnNw$;Idv*i-*EW7Lsq7OxfKPAVPX7nDDx$*$(kG zXvKihrIa?wagD5%GckjeLbPEL0ouN87L1(d77aGXmL|8-9BcN9l(w*ES#Li@P%Wi6 zKuk;u*-auc=P>JsD3K#|z=FFPbF)=a*h;b5Ry}iH0E-y(_9LzNn5tbjSi#)7j z_xqNaUDNLt2CWvNF&HJ2u}UZ#9L)0`LZvcfWyRMWwr_c?_C!rZJy9*CTELQ>G4YW| zA}sbo;WlP=&~!3!d0M?ZBtN9!KVy_LCO(acIw9<$4v0`6!Ix#o`q{HkpjA*&r>&gp zS>*HhAjwE)$}OpqLJ+;%Wg;_t0Af<_>ZQRo_*Qx-G68iVhT#w+v4@>x}Hxu({ z{3T5zXCBJoyU$_Vr~~AbreI>YCIh5$ zwXIf;DX1eHss>B+h8lA*wyr-jx*!IBeaiiD>{<6;hNWtgKq=vG#G}5H_VP?)}KJB{Vswr38^n1>2`JVL~yZN9snwsJ>H1 zX;@?(k!2q<2{Kl>?254vb;phZ);Ud*oWffZyk%kfh8%-KZ(rVal2{qD8`7915!GAt zSk-Mtz-Yz%<#Xni+cHkU-L#d_3hMfz_wv4Y#zA8r=RKTFhl5oM*ftLbH|9Oze>Cml z$y4YH9HY`6%a9K@9i?Yu^fe|C9K^~%fK5l~)k~XUM(Nc{8vr&N^HGm&G_7%m1N4#7 z<}aUTdv<1^VEJ?A#_BD4k#vY?*e%e21q-mT3K+#;a^6D{9kTIx(_2H}bt$E(ebz+) z=*cb>FuRLSNUMQN!eP^n=Ojk4hET5{obf$-NUAR2GZBXju-X=T1ZWpgHD@pdEiNSy zxq+}8KBe(B*rLAH?3kXSphj9cYhWhJK!jD*T8|#BYcNLe5=37F$b-x%PH?O~xB$Sa zycR(5gOrAzKnf?CDK7FqDzPU|X-{DtfoTy? z!ZZ@<6tz#qNj6#1bm|mt=_Hj@r%nlT2^}S=Mu(|WRQnwE_Bv%mv_|_JOm;$Fq;-zv zIaZ;aP1CqaNA()3Mzb~s80?jAYBgG;7uj^KzDyr$6v|Tsf2y)IxSw)Yr(Py4go)S7 zxm<-7Z4PYEK!LbSPrli0J+&nqh{KF)<#N;LC?=&e`Y5sWL=X&cdJGs{N>XAoKn!o8 zj3}X*;t#ly$>=F-FxSLQyHDaadw;M*3ws+Ff`Vfkqm@pKpo-^FX)Do}p=QlSIv@@^(HChUO#5W?r8>SmZ1j~doa;nifg^UJuSTU- zaGiA3iM|4lVFyWqrD38k)!8infDxZJJNiLH{N(cQOq_VPuF zFU-VX-Dz*~OpNL{Q)o$W5^zTcZ#r|P35BR*nk zhmCpR_`J>>bnU9XFV#s<8m4YzEe8y%b3iujIpi`)p9AZi3n8kH(0@2trJACr+uq%` z4^BSe)7j;`Ak2xS?nLxjGwUSj6&a>?uBFeGgM3I_?2GZCg$2nPr6fifIh|>-QYY$6 z+pD9mz{A^h5cYkwXZ}tX8KRZGuO^{Lbqr2`olu0d5T;HjGOYiTvfk`aWLhT_X>J!X z3LhwJFi$c`1ltKkNUdQ?GE=wx1BIPXB(w;QI>T_yO9uRxg(5sNx2Ic#43In&5u8NX zCcsW8LRttDg(7p@?3($sQcYA@4UXx_y}wsvUu@B1^@^YcAKI56Va16C_Kd$wPYu(e zS08yC0WEsWFbckBt=FNt$ent*9rJBs7GIObNiu3CM|>nk)o-WIA-+V2T#Z!p?dfzN z1MF!+Bj5h7Q~=Nj~Sad&*nzYZP&Q#)w)l(u{zB@8o6+(3l`Uj==hSg;^_uIFMw(ITI# zOYB5UjKg8#1_D$aQg8Y(2?4(4`&Kj^CL|738$Py^DO9qSfyoY7WgAQFv9@He2$A|w zCSPi=28L%l`l}Aewu3ODb0tDF2B>Nj4p_UwUI40Y_V2dZ4l)lZH*LJHb>pc8SDz=VL~6RY{@e*GBDqjfTn7xa~PJ` zK@rVM?4SS#XzP}yU@s%fb1~Ti`TBXuP)v@y?wvB9O)|1g2{Lsh7BxQST)mn(^92P7m?+;Ly zU+n5v>VUMEF1E1D)$zvd3@fWlGKG~LENxb{Yl>CrG)bxeTAa2}puI_Ax4ueg^Lpp^ zoWHQu&Ec@naabmVe1%8q%lpS0H$cQ{{E9l{fMPXHM2cAAKOb3w4((+ND{HozmAf>} zTWjzF#~g|?O-iGUq+3rB!1&$a)cr8vFes&~Txbsfu`dZ$FRem@Y19E(e8_`9WA%XL z7}c?fdU=SYegk0xX9*8fIc(j0E#ApSEAqvB+qxY{w&6!BDBfHk%k;3fVG1~@Wt;#z zqZOouFrnWDwcOxt@--)k~{jPr?Qda zkE2BXs^%01xHbY%dzb{UX>^lBhgC~f)N?rTw5iw(qg5@^^LD#+Xt!I3O=Ac3R)E}S zqE>n0JPBag#MY{-vC&oCw!8pUZ-ApsKY&eJ(Y4~nF`0u{E{ z?OWH3be6o50w>X;OBP?9N~+owQ(8YZz3<1aS92_XwPiW|(re2>)gRR(J>C8SmTfFv z%~)`j(qiI`C`&cR^kC7%HbdCcIrdA%By~#L%R|RnE&GV;D}$7#zQh?HUI?^SeUPO1 zDw)_jvQfuKl4`W%*+nei9dDH+9PaH^K$}bKeQ<2FCTs8c=|7+))n+U2xr{XkOyu^~ z%WW39;y`yd4zP>*0t_OfMYM=8vwSkXQ^q(PEE79ZV_vl@)z@7CBa0E zIru3))hLxL>}VFZy(&dH>XAbPkT&JTo)cBkl6ri|HsLDK0*zHb@q+o8)jRGnfZ8fI}|98js}z6;%v<-w?kwa4m4$-23s z*+4ylE$ds<+rO;$T%&jKyoMc3H1XHXRvKwAiQD8b z>mX)NIH^nZoRTlL*(iCZQYyId`7V)B2GAUq^2PkBi`?ctxzKw^i$r4^ni%frl6v8= zq8z&fv&2_-Y;4GG7qOWx?K8tyRjguH9k6NhJ4x2#t!2Z`G4IBZoa19f;UMd%Te*a_ z{jqtlS&Hs2BsG|X!w?YD&<5*VQr*fmfLLU_2;`z1z6PM;X2i9C%-SKM`uBS1AaJlc zHs7WkCc(&R9z77AWikg%AaWE?ZjTxmIYEidImEE~05;w<9pWz9W32Aj95CFRy-sR1 zp)sV+L=FLKq)y5qvvt9(-+@ydkowufj8S)_p7q#9bF>%g0H$lYbvwi*;BbwD9i>rr zg{~51Q+;s3w=(Umqz3(l3$PKV>LL$i&8UdA$4LOaT9D{@0>EOU3XkEc0*-*DnV9DP zrZ&w`*7{b}_XXl4u5C8VkU#mM>b$S^2TRR%$|Me3pDIWyz^tm84d|m@jbjAuizE_3>6pMr~|s`6t-=3SR%ACWMFH0zd9;b*CC&Y7eAYknAt6#PaX`KWW)bf`5+ z@nv})_Tfdaw8G4|X-rKr5G$^=JeR_3_-KaDb#KkaY%Vd|Qx@Goe4|?IqaOGm9jIqL z{`E2e%3M-~xF>-8+s7)nBpmW(d1N6_B`z0PwW}&OXrwdcx{M}~TQ!EP3(!T7cU;dh z$gTt%7+WB!`=MPdT3H)8uCT`7O$^(qDauqpHEku^1P>)E8hK{k5=D!Yxb>1Xs7BF? zYhfc@Dy-um)pZyEW-aP}ww#PC02!C+fR1H)hEt($D@xb~5K~UEO(<=%#-&!#A4!c+ z1W0SjO{&p}6>V_107N?VES~sCU~>i&iGi01t9Qm$uY)5OqfMo zWCqB6$4B+fYE|xn39Sr6)OS#Yi0U~7^X#2+$;B6EvDrmkWSSrz^ofL=$s=_RQN*N) zzanaT8{1^hGfkV|>TNv&tZ7Bj zTVzTia?>xE)$=kYCwp0|VRsc8;C)N`Tm7hss!_}(#+m`Ll9Q4Kzez-H*hV4Ll(xu4 zVI+NT{_>>=<(n9so{>gjAG@<*KGxysTQtB)ZnF#nj!iJC{uK{6Xa+h>A~GTdvyopX zsP%?vF{{@p$jc_`P9tL>qQl5f*2brNAD@oC|D$pc3 z#b50jX#5m|3B}k8ql5^GH7$~e%obCrG#Fu4cbhiPs%mWig#TfYLwlA1qGyR^6 zzWKW|yl3^FXE%-8@zdM&Umy5`v&TQWL4V}>TaR9Oz{UFieC(q8ue{aMe?N1m*ZS}O z(*O9_H$Jyw_t*44m0$UffBnFb`d>yb_}n8`zf1q?C-0m;dEIUL6J0l@mychsKRNs9 zPu}*6S^87T->_-Jf@8n`dwpa2>1#gty}R{|>yEuQ zGjo4^)0sCvf8Ffw>YHvg4*11AKhrlKls{rL_}gsXdF0vWy-nZztNG(@xU5g#a?&wd zo_qOk`j)FsUpMEC7xgWFyzxt4yjgUL3e?(v|Pg zcaJyT`{l#%d)kCgpE|Ku-@WLrhaUOG&HC=kj(+}AFCU}t&MkWUqMIJT@8WOwJo05+ zKdxW@^jC1-kAK_$@WXiSlvA&~@g(5nF8k>V-vI8a6-VCx5xnP|#h*A~D&D(g+#zFU zf{&HAer+fNzP|U)3ySxE&tG1wy)XxSf9r}jUfu*bE?d2I=po3pe&MKVehWFLzBXaP z{*e3Jf4X#D8hZH9^O?VYQQvx9=JfOW&egX(_`v1opE5 Date: Wed, 27 Sep 2023 10:30:49 -0600 Subject: [PATCH 2/8] Adapt other interface kernels as well. (#25599) --- .../FVScalarLagrangeMultiplierInterface.h | 3 +++ .../FVScalarLagrangeMultiplierInterface.C | 21 ++++++++++++------- .../FVOnlyAddDiffusionToOneSideOfInterface.C | 14 +++++++++---- 3 files changed, 27 insertions(+), 11 deletions(-) diff --git a/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h b/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h index c407a2ef20c8..5be24aefb238 100644 --- a/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h +++ b/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h @@ -33,4 +33,7 @@ class FVScalarLagrangeMultiplierInterface : public FVInterfaceKernel /// The Lagrange Multiplier value const ADVariableValue & _lambda; + + /// Lambda assembly + Assembly & _lambda_assembly; }; diff --git a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C index 595bb7200905..d8f2ee099fed 100644 --- a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C +++ b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C @@ -26,7 +26,8 @@ FVScalarLagrangeMultiplierInterface::FVScalarLagrangeMultiplierInterface( const InputParameters & params) : FVInterfaceKernel(params), _lambda_var(*getScalarVar("lambda", 0)), - _lambda(adCoupledScalarValue("lambda")) + _lambda(adCoupledScalarValue("lambda")), + _lambda_assembly(_subproblem.assembly(_tid, _lambda_var.sys().number())) { } @@ -38,17 +39,20 @@ FVScalarLagrangeMultiplierInterface::computeResidual(const FaceInfo & fi) const auto var_elem_num = _elem_is_one ? var1().number() : var2().number(); const auto var_neigh_num = _elem_is_one ? var2().number() : var1().number(); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + const auto r = MetaPhysicL::raw_value(_lambda[0]) * fi.faceArea() * fi.faceCoord() * (2 * _elem_is_one - 1); // Primal residual - addResidual(r, var_elem_num, false); - addResidual(-r, var_neigh_num, true); + addResidual(r, var_elem_num, assembly_elem, false); + addResidual(-r, var_neigh_num, assembly_neigh, true); // LM residual. We may not have any actual ScalarKernels in our simulation so we need to manually // make sure the scalar residuals get cached for later addition const auto lm_r = MetaPhysicL::raw_value(computeQpResidual()) * fi.faceArea() * fi.faceCoord(); - addResiduals(_assembly, + addResiduals(_lambda_assembly, std::array{{lm_r}}, _lambda_var.dofIndices(), _lambda_var.scalingFactor()); @@ -67,17 +71,20 @@ FVScalarLagrangeMultiplierInterface::computeJacobian(const FaceInfo & fi) const auto elem_scaling_factor = _elem_is_one ? var1().scalingFactor() : var2().scalingFactor(); const auto neigh_scaling_factor = _elem_is_one ? var2().scalingFactor() : var1().scalingFactor(); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + // Primal const auto primal_r = _lambda[0] * fi.faceArea() * fi.faceCoord() * (2 * _elem_is_one - 1); addResidualsAndJacobian( - _assembly, std::array{{primal_r}}, elem_dof_indices, elem_scaling_factor); + assembly_elem, std::array{{primal_r}}, elem_dof_indices, elem_scaling_factor); addResidualsAndJacobian( - _assembly, std::array{{-primal_r}}, neigh_dof_indices, neigh_scaling_factor); + assembly_neigh, std::array{{-primal_r}}, neigh_dof_indices, neigh_scaling_factor); // LM const auto lm_r = computeQpResidual() * fi.faceArea() * fi.faceCoord(); mooseAssert(_lambda_var.dofIndices().size() == 1, "We should only have one dof"); - addResidualsAndJacobian(_assembly, + addResidualsAndJacobian(_lambda_assembly, std::array{{lm_r}}, _lambda_var.dofIndices(), _lambda_var.scalingFactor()); diff --git a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C index c98cafe4e0a5..9cd539132328 100644 --- a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C +++ b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C @@ -39,10 +39,13 @@ FVOnlyAddDiffusionToOneSideOfInterface::computeResidual(const FaceInfo & fi) const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual()); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + if (_elem_is_one) - addResidual(r, var_elem_num, false); + addResidual(r, var_elem_num, assembly_elem, false); else - addResidual(-r, var_neigh_num, true); + addResidual(-r, var_neigh_num, assembly_neigh, true); } void @@ -65,14 +68,17 @@ FVOnlyAddDiffusionToOneSideOfInterface::computeJacobian(const FaceInfo & fi) mooseAssert((elem_dof_indices.size() == 1) && (neigh_dof_indices.size() == 1), "We're currently built to use CONSTANT MONOMIALS"); + auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; + auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; + const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); if (_elem_is_one) addResidualsAndJacobian( - _assembly, std::array{{r}}, elem_dof_indices, elem_scaling_factor); + assembly_elem, std::array{{r}}, elem_dof_indices, elem_scaling_factor); else addResidualsAndJacobian( - _assembly, std::array{{-r}}, neigh_dof_indices, neighbor_scaling_factor); + assembly_neigh, std::array{{-r}}, neigh_dof_indices, neighbor_scaling_factor); } ADReal From fe88b3f652dd1d406173fbbd5787451b2fc42782 Mon Sep 17 00:00:00 2001 From: Peter German Date: Wed, 27 Sep 2023 12:41:57 -0600 Subject: [PATCH 3/8] Smooth residual and jacobian computation for FVIKs. (#25599) --- framework/include/fviks/FVInterfaceKernel.h | 33 ++----- .../FVScalarLagrangeMultiplierInterface.h | 3 - framework/src/fviks/FVInterfaceKernel.C | 96 +++++-------------- .../FVScalarLagrangeMultiplierInterface.C | 23 ++--- framework/src/problems/FEProblemBase.C | 21 ++-- .../FVOnlyAddDiffusionToOneSideOfInterface.C | 16 ++-- 6 files changed, 55 insertions(+), 137 deletions(-) diff --git a/framework/include/fviks/FVInterfaceKernel.h b/framework/include/fviks/FVInterfaceKernel.h index c632cf95a882..f27171f08905 100644 --- a/framework/include/fviks/FVInterfaceKernel.h +++ b/framework/include/fviks/FVInterfaceKernel.h @@ -100,14 +100,10 @@ class FVInterfaceKernel : public MooseObject, const std::set & sub2() const { return _subdomain2; } /** - * @return The system associated with the first variable + * @return The system associated with this object. Either an undisplaced or displaced nonlinear + * system */ - const SystemBase & sys() const { return _sys1; } - - /** - * @return The system associated with the second variable - */ - const SystemBase & sys2() const { return _sys2; } + const SystemBase & sys() const { return _var1.sys(); } /** * @return Whether the \p FaceInfo element is on the 1st side of the interface @@ -137,16 +133,13 @@ class FVInterfaceKernel : public MooseObject, /** * Process the provided residual given \p var_num and whether this is on the neighbor side */ - void addResidual(Real resid, unsigned int var_num, Assembly & assembly, bool neighbor); + void addResidual(Real resid, unsigned int var_num, bool neighbor); using TaggingInterface::addJacobian; /** * Process the derivatives for the provided residual and dof index */ - void addJacobian(const ADReal & resid, - dof_id_type dof_index, - Assembly & assembly, - Real scaling_factor); + void addJacobian(const ADReal & resid, dof_id_type dof_index, Real scaling_factor); /** * @return A structure that contains information about the face info element and skewness @@ -199,23 +192,11 @@ class FVInterfaceKernel : public MooseObject, /// The SubProblem SubProblem & _subproblem; - /// the system object for variable 1 - SystemBase & _sys1; - - /// the system object for variable 2 - SystemBase & _sys2; - - /// Variable on one side of the interface MooseVariableFV & _var1; - - /// Variable on the other side of the interface MooseVariableFV & _var2; - /// The Assembly object for system 1 - Assembly & _assembly1; - - /// The Assembly object for system 2 - Assembly & _assembly2; + /// The Assembly object + Assembly & _assembly; private: std::set _subdomain1; diff --git a/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h b/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h index 5be24aefb238..c407a2ef20c8 100644 --- a/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h +++ b/framework/include/fviks/FVScalarLagrangeMultiplierInterface.h @@ -33,7 +33,4 @@ class FVScalarLagrangeMultiplierInterface : public FVInterfaceKernel /// The Lagrange Multiplier value const ADVariableValue & _lambda; - - /// Lambda assembly - Assembly & _lambda_assembly; }; diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index 16c661134de4..206e2e2dbfcc 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -98,19 +98,19 @@ FVInterfaceKernel::FVInterfaceKernel(const InputParameters & parameters) ADFunctorInterface(this), _tid(getParam("_tid")), _subproblem(*getCheckedPointerParam("_subproblem")), - _sys1(_subproblem.getVariable(_tid, getParam("variable1")).sys()), - _sys2(_subproblem + _var1(_subproblem.getVariable(_tid, getParam("variable1")) + .sys() + .getFVVariable(_tid, getParam("variable1"))), + _var2(_subproblem .getVariable(_tid, isParamValid("variable2") ? getParam("variable2") : getParam("variable1")) - .sys()), - _var1(_sys1.getFVVariable(_tid, getParam("variable1"))), - _var2(_sys2.getFVVariable(_tid, - isParamValid("variable2") - ? getParam("variable2") - : getParam("variable1"))), - _assembly1(_subproblem.assembly(_tid, _sys1.number())), - _assembly2(_subproblem.assembly(_tid, _sys2.number())), + .sys() + .getFVVariable(_tid, + isParamValid("variable2") + ? getParam("variable2") + : getParam("variable1"))), + _assembly(_subproblem.assembly(_tid, _var1.sys().number())), _mesh(_subproblem.mesh()) { if (getParam("use_displaced_mesh")) @@ -163,12 +163,9 @@ FVInterfaceKernel::setupData(const FaceInfo & fi) } void -FVInterfaceKernel::addResidual(const Real resid, - const unsigned int var_num, - Assembly & assembly, - const bool neighbor) +FVInterfaceKernel::addResidual(const Real resid, const unsigned int var_num, const bool neighbor) { - neighbor ? prepareVectorTagNeighbor(assembly, var_num) : prepareVectorTag(assembly, var_num); + neighbor ? prepareVectorTagNeighbor(_assembly, var_num) : prepareVectorTag(_assembly, var_num); _local_re(0) = resid; accumulateTaggedLocalResidual(); } @@ -176,10 +173,9 @@ FVInterfaceKernel::addResidual(const Real resid, void FVInterfaceKernel::addJacobian(const ADReal & resid, const dof_id_type dof_index, - Assembly & assembly, const Real scaling_factor) { - addJacobian(assembly, + addJacobian(_assembly, std::array{{resid}}, std::array{{dof_index}}, scaling_factor); @@ -190,32 +186,11 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi) { setupData(fi); - const auto var_elem_num = _elem_is_one ? _var1.number() : _var2.number(); - const auto var_neigh_num = _elem_is_one ? _var2.number() : _var1.number(); - - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; - const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual()); - bool is_var1_system = - (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()); - - // std::cout << "Adding this " - - // if (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()) - // { - // addResidual(_elem_is_one ? r : -r, _var1.number(), _assembly1, _elem_is_one ? true : false); - // addResidual(-r, var_neigh_num, assembly_neigh, true); - // } - // else - // addResidual(-r, var_neigh_num, assembly_neigh, true); - - // if (&_sys1 == &_sys2) - // addResidual(r, var_neigh_num, assembly_neigh, true); - - addResidual(r, var_elem_num, assembly_elem, false); - addResidual(-r, var_neigh_num, assembly_neigh, true); + addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true); + if (_var1.sys().number() == _var2.sys().number()) + addResidual(_elem_is_one ? -r : r, _var2.number(), _elem_is_one ? true : false); } void @@ -232,40 +207,21 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi) const auto & elem_dof_indices = _elem_is_one ? _var1.dofIndices() : _var2.dofIndices(); const auto & neigh_dof_indices = _elem_is_one ? _var2.dofIndicesNeighbor() : _var1.dofIndicesNeighbor(); + mooseAssert((elem_dof_indices.size() == 1) && (neigh_dof_indices.size() == 1), "We're currently built to use CONSTANT MONOMIALS"); - const auto elem_scaling_factor = _elem_is_one ? _var1.scalingFactor() : _var2.scalingFactor(); - const auto neigh_scaling_factor = _elem_is_one ? _var2.scalingFactor() : _var1.scalingFactor(); - - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); - bool is_var1_system = - (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()); - - // if (_subproblem.systemBaseNonlinear(_subproblem.currentNlSysNum()).number() == _sys1.number()) - // { - // addResidualsAndJacobian(_assembly1, - // std::array{{_elem_is_one ? r : -r}}, - // _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), - // _var1.scalingFactor()); - // addResidualsAndJacobian( - // assembly_neigh, std::array{{-r}}, neigh_dof_indices, neigh_scaling_factor); - // } - // else - // addResidualsAndJacobian( - // assembly_neigh, std::array{{r}}, neigh_dof_indices, neigh_scaling_factor); - - // if (&_sys1 == &_sys2) - // addResidualsAndJacobian( - // assembly_neigh, std::array{{r}}, neigh_dof_indices, neigh_scaling_factor); - - addResidualsAndJacobian( - _assembly1, std::array{{r}}, elem_dof_indices, elem_scaling_factor); - addResidualsAndJacobian( - _assembly1, std::array{{-r}}, neigh_dof_indices, neigh_scaling_factor); + addResidualsAndJacobian(_assembly, + std::array{{_elem_is_one ? r : -r}}, + _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), + _var1.scalingFactor()); + if (_var1.sys().number() == _var2.sys().number()) + addResidualsAndJacobian(_assembly, + std::array{{_elem_is_one ? -r : r}}, + _elem_is_one ? _var2.dofIndicesNeighbor() : _var2.dofIndices(), + _var2.scalingFactor()); } Moose::ElemArg diff --git a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C index d8f2ee099fed..5f66c6fe7c7b 100644 --- a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C +++ b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C @@ -26,9 +26,10 @@ FVScalarLagrangeMultiplierInterface::FVScalarLagrangeMultiplierInterface( const InputParameters & params) : FVInterfaceKernel(params), _lambda_var(*getScalarVar("lambda", 0)), - _lambda(adCoupledScalarValue("lambda")), - _lambda_assembly(_subproblem.assembly(_tid, _lambda_var.sys().number())) + _lambda(adCoupledScalarValue("lambda")) { + if (_var1.sys().number() != _var2.sys().number()) + mooseError(this->type(), " does not support multiple nonlinear systems!"); } void @@ -39,20 +40,17 @@ FVScalarLagrangeMultiplierInterface::computeResidual(const FaceInfo & fi) const auto var_elem_num = _elem_is_one ? var1().number() : var2().number(); const auto var_neigh_num = _elem_is_one ? var2().number() : var1().number(); - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; - const auto r = MetaPhysicL::raw_value(_lambda[0]) * fi.faceArea() * fi.faceCoord() * (2 * _elem_is_one - 1); // Primal residual - addResidual(r, var_elem_num, assembly_elem, false); - addResidual(-r, var_neigh_num, assembly_neigh, true); + addResidual(r, var_elem_num, false); + addResidual(-r, var_neigh_num, true); // LM residual. We may not have any actual ScalarKernels in our simulation so we need to manually // make sure the scalar residuals get cached for later addition const auto lm_r = MetaPhysicL::raw_value(computeQpResidual()) * fi.faceArea() * fi.faceCoord(); - addResiduals(_lambda_assembly, + addResiduals(_assembly, std::array{{lm_r}}, _lambda_var.dofIndices(), _lambda_var.scalingFactor()); @@ -71,20 +69,17 @@ FVScalarLagrangeMultiplierInterface::computeJacobian(const FaceInfo & fi) const auto elem_scaling_factor = _elem_is_one ? var1().scalingFactor() : var2().scalingFactor(); const auto neigh_scaling_factor = _elem_is_one ? var2().scalingFactor() : var1().scalingFactor(); - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; - // Primal const auto primal_r = _lambda[0] * fi.faceArea() * fi.faceCoord() * (2 * _elem_is_one - 1); addResidualsAndJacobian( - assembly_elem, std::array{{primal_r}}, elem_dof_indices, elem_scaling_factor); + _assembly, std::array{{primal_r}}, elem_dof_indices, elem_scaling_factor); addResidualsAndJacobian( - assembly_neigh, std::array{{-primal_r}}, neigh_dof_indices, neigh_scaling_factor); + _assembly, std::array{{-primal_r}}, neigh_dof_indices, neigh_scaling_factor); // LM const auto lm_r = computeQpResidual() * fi.faceArea() * fi.faceCoord(); mooseAssert(_lambda_var.dofIndices().size() == 1, "We should only have one dof"); - addResidualsAndJacobian(_lambda_assembly, + addResidualsAndJacobian(_assembly, std::array{{lm_r}}, _lambda_var.dofIndices(), _lambda_var.scalingFactor()); diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index f51cc8587fcc..9e6f36ea5cb6 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -2929,25 +2929,18 @@ FEProblemBase::addFVInterfaceKernel(const std::string & fv_ik_name, const std::string & name, InputParameters & parameters) { - const auto nl_sys_num1 = + /// Checking the nonlinear system number of variable1 + /// If variable1 and variable2 live on different systems, the interface kernel + /// is added to the system of variable1 alone. Another interface kernel needs to be + /// created with flipped variables for the system of variable2 + const auto nl_sys_num = determineNonlinearSystem(parameters.varName("variable1", name), true).first ? determineNonlinearSystem(parameters.varName("variable1", name), true).second : (unsigned int)0; - const auto nl_sys_num2 = - determineNonlinearSystem(parameters.varName("variable2", name), true).first - ? determineNonlinearSystem(parameters.varName("variable2", name), true).second - : (unsigned int)0; - - parameters.set("_sys") = _nl[nl_sys_num1].get(); + /// This is set because the Attribute system filters IKs based on this + parameters.set("_sys") = _nl[nl_sys_num].get(); addObject(fv_ik_name, name, parameters); - - // if (nl_sys_num1 != nl_sys_num2) - // { - // InputParameters params2 = parameters; - // params2.set("_sys") = _nl[nl_sys_num2].get(); - // addObject(fv_ik_name, name + "_var2", params2); - // } } // InterfaceKernels //// diff --git a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C index 9cd539132328..d40663c794f6 100644 --- a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C +++ b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C @@ -27,6 +27,8 @@ FVOnlyAddDiffusionToOneSideOfInterface::FVOnlyAddDiffusionToOneSideOfInterface( const InputParameters & params) : FVInterfaceKernel(params), _coeff2(getFunctor("coeff2")) { + if (_var1.sys().number() != _var2.sys().number()) + mooseError(this->type(), " does not support multiple nonlinear systems!"); } void @@ -39,13 +41,10 @@ FVOnlyAddDiffusionToOneSideOfInterface::computeResidual(const FaceInfo & fi) const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual()); - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; - if (_elem_is_one) - addResidual(r, var_elem_num, assembly_elem, false); + addResidual(r, var_elem_num, false); else - addResidual(-r, var_neigh_num, assembly_neigh, true); + addResidual(-r, var_neigh_num, true); } void @@ -68,17 +67,14 @@ FVOnlyAddDiffusionToOneSideOfInterface::computeJacobian(const FaceInfo & fi) mooseAssert((elem_dof_indices.size() == 1) && (neigh_dof_indices.size() == 1), "We're currently built to use CONSTANT MONOMIALS"); - auto & assembly_elem = _elem_is_one ? _assembly1 : _assembly2; - auto & assembly_neigh = _elem_is_one ? _assembly2 : _assembly1; - const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); if (_elem_is_one) addResidualsAndJacobian( - assembly_elem, std::array{{r}}, elem_dof_indices, elem_scaling_factor); + _assembly, std::array{{r}}, elem_dof_indices, elem_scaling_factor); else addResidualsAndJacobian( - assembly_neigh, std::array{{-r}}, neigh_dof_indices, neighbor_scaling_factor); + _assembly, std::array{{-r}}, neigh_dof_indices, neighbor_scaling_factor); } ADReal From 4b67dcb60c6311336ba911232b775a4aef15efa3 Mon Sep 17 00:00:00 2001 From: Peter German Date: Wed, 27 Sep 2023 13:49:56 -0600 Subject: [PATCH 4/8] Modify convectioncorreclationik. (#25599) --- framework/src/problems/FEProblemBase.C | 8 ++--- .../fviks/FVConvectionCorrelationInterface.C | 30 +++++++++++-------- test/src/executioners/SteadySolve2.C | 11 ++++--- test/tests/fviks/diffusion/tests | 2 +- 4 files changed, 26 insertions(+), 25 deletions(-) diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index 9e6f36ea5cb6..a7901dc9ce5d 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -176,7 +176,7 @@ FEProblemBase::validParams() "True to skip the NonlinearSystem check for work to do (e.g. Make sure " "that there are variables to solve for)."); params.addParam("allow_initial_conditions_with_restart", - true, + false, "True to allow the user to specify initial conditions when restarting. " "Initial conditions can override any restarted field"); @@ -6033,16 +6033,12 @@ FEProblemBase::computeResidualAndJacobian(const NumericVector & soln, { // vector tags { - _current_nl_sys->associateVectorToTag(residual, _current_nl_sys->residualVectorTag()); const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL); _fe_vector_tags.clear(); for (const auto & residual_vector_tag : residual_vector_tags) - // We filter out tags which do not have associated vectors in the current nonlinear - // system. This is essential to be able to use system-dependent residual tags. - if (_current_nl_sys->hasVector(residual_vector_tag._id)) - _fe_vector_tags.insert(residual_vector_tag._id); + _fe_vector_tags.insert(residual_vector_tag._id); setCurrentResidualVectorTags(_fe_vector_tags); } diff --git a/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C b/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C index 18875f45f390..3e037e74c6e4 100644 --- a/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C +++ b/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C @@ -43,16 +43,20 @@ FVConvectionCorrelationInterface::FVConvectionCorrelationInterface(const InputPa mooseError( "The bulk distance should be specified or 'wall_cell_is_bulk' should be set to true for " "the FVTwoVarConvectionCorrelationInterface"); - if ("wraps_" + var1().name() != _temp_fluid.functorName()) - paramError(NS::T_fluid, "variable1 must be equal to T_fluid parameter."); - if ("wraps_" + var2().name() != _temp_solid.functorName()) - paramError(NS::T_solid, "variable2 must be equal to T_solid parameter."); } ADReal FVConvectionCorrelationInterface::computeQpResidual() { - const Elem * current_elem = elemIsOne() ? &_face_info->elem() : _face_info->neighborPtr(); + // If variable1 is fluid and variable 1 is on elem or + // if variable2 is fluid and variable 2 is on elem + // the fluid element will be elem otherwise it is the neighbor + bool var1_is_fluid = ("wraps_" + var1().name() == _temp_fluid.functorName()); + const Elem * elem_on_fluid_side = + (elemIsOne() && var1_is_fluid) || (!elemIsOne() && !var1_is_fluid) + ? &_face_info->elem() + : _face_info->neighborPtr(); + const Elem * bulk_elem; const auto state = determineState(); if (!_use_wall_cell) @@ -60,29 +64,31 @@ FVConvectionCorrelationInterface::computeQpResidual() Point p = _face_info->faceCentroid(); Point du = Point(MetaPhysicL::raw_value(_normal)); du *= _bulk_distance; - if (elemIsOne()) + // The normal always points outwards from the elem (towards the neighbor) + if (elem_on_fluid_side == &_face_info->elem()) p -= du; else p += du; bulk_elem = (*_pl)(p); } else - bulk_elem = current_elem; + bulk_elem = elem_on_fluid_side; mooseAssert(bulk_elem, "The element at bulk_distance from the wall was not found in the mesh. " "Increase the number of ghost layers with the 'ghost_layers' parameter."); - mooseAssert(var1().hasBlocks(bulk_elem->subdomain_id()), + mooseAssert((var1_is_fluid ? var1() : var2()).hasBlocks(bulk_elem->subdomain_id()), "The fluid temperature is not defined at bulk_distance from the wall."); - const auto face_arg_side1 = singleSidedFaceArg(var1(), _face_info); - const auto face_arg_side2 = singleSidedFaceArg(var2(), _face_info); + const auto fluid_side = singleSidedFaceArg(var1_is_fluid ? var1() : var2(), _face_info); + const auto solid_side = singleSidedFaceArg(var1_is_fluid ? var2() : var1(), _face_info); + const auto bulk_elem_arg = makeElemArg(bulk_elem); /// We make sure that the gradient*normal part is addressed auto multipler = _normal * (_face_info->faceCentroid() - bulk_elem->vertex_average()) > 0 ? 1 : -1; - return multipler * _htc(face_arg_side1, state) * - (_temp_fluid(bulk_elem_arg, state) - _temp_solid(face_arg_side2, state)); + return multipler * _htc(fluid_side, state) * + (_temp_fluid(bulk_elem_arg, state) - _temp_solid(solid_side, state)); } diff --git a/test/src/executioners/SteadySolve2.C b/test/src/executioners/SteadySolve2.C index 8c98875c1d79..53cd593fd578 100644 --- a/test/src/executioners/SteadySolve2.C +++ b/test/src/executioners/SteadySolve2.C @@ -117,12 +117,11 @@ SteadySolve2::execute() for (unsigned int i = 0; i < _number_of_iterations; i++) { converged = converged && solve(_first_nl_sys) && solve(_second_nl_sys); - } - - if (!converged) - { - _console << "Aborting as solve did not converge" << std::endl; - break; + if (!converged) + { + _console << "Aborting as solve did not converge" << std::endl; + break; + } } _problem.computeIndicators(); diff --git a/test/tests/fviks/diffusion/tests b/test/tests/fviks/diffusion/tests index 73ee6e78077c..4dcda059f759 100644 --- a/test/tests/fviks/diffusion/tests +++ b/test/tests/fviks/diffusion/tests @@ -20,7 +20,7 @@ "Outputs/file_base='harmonic'" [] [diffusion-multisystem] - issues = '#22356 #8780' + issues = '#25599' design = 'FVDiffusionInterface.md Problem/index.md' requirement = 'The system shall be able to solve a block-restricted diffusion problem where the variables live on different nonlinear systems.' type = Exodiff From dc527f60db54e1d3943dfdcdfbbcde39618e7029 Mon Sep 17 00:00:00 2001 From: Peter German Date: Wed, 27 Sep 2023 13:54:01 -0600 Subject: [PATCH 5/8] Remove unused variables, clean up source code. (#25599) --- framework/include/fviks/FVInterfaceKernel.h | 6 ------ framework/src/fviks/FVInterfaceKernel.C | 7 ------- framework/src/fviks/FVScalarLagrangeMultiplierInterface.C | 2 +- test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C | 2 +- 4 files changed, 2 insertions(+), 15 deletions(-) diff --git a/framework/include/fviks/FVInterfaceKernel.h b/framework/include/fviks/FVInterfaceKernel.h index f27171f08905..ff6ad444a88b 100644 --- a/framework/include/fviks/FVInterfaceKernel.h +++ b/framework/include/fviks/FVInterfaceKernel.h @@ -99,12 +99,6 @@ class FVInterfaceKernel : public MooseObject, */ const std::set & sub2() const { return _subdomain2; } - /** - * @return The system associated with this object. Either an undisplaced or displaced nonlinear - * system - */ - const SystemBase & sys() const { return _var1.sys(); } - /** * @return Whether the \p FaceInfo element is on the 1st side of the interface */ diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index 206e2e2dbfcc..cf5c0ce883fc 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -204,13 +204,6 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi) { setupData(fi); - const auto & elem_dof_indices = _elem_is_one ? _var1.dofIndices() : _var2.dofIndices(); - const auto & neigh_dof_indices = - _elem_is_one ? _var2.dofIndicesNeighbor() : _var1.dofIndicesNeighbor(); - - mooseAssert((elem_dof_indices.size() == 1) && (neigh_dof_indices.size() == 1), - "We're currently built to use CONSTANT MONOMIALS"); - const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); addResidualsAndJacobian(_assembly, diff --git a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C index 5f66c6fe7c7b..63e7db8422d9 100644 --- a/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C +++ b/framework/src/fviks/FVScalarLagrangeMultiplierInterface.C @@ -28,7 +28,7 @@ FVScalarLagrangeMultiplierInterface::FVScalarLagrangeMultiplierInterface( _lambda_var(*getScalarVar("lambda", 0)), _lambda(adCoupledScalarValue("lambda")) { - if (_var1.sys().number() != _var2.sys().number()) + if (var1().sys().number() != var2().sys().number()) mooseError(this->type(), " does not support multiple nonlinear systems!"); } diff --git a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C index d40663c794f6..125c10a174f3 100644 --- a/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C +++ b/test/src/fviks/FVOnlyAddDiffusionToOneSideOfInterface.C @@ -27,7 +27,7 @@ FVOnlyAddDiffusionToOneSideOfInterface::FVOnlyAddDiffusionToOneSideOfInterface( const InputParameters & params) : FVInterfaceKernel(params), _coeff2(getFunctor("coeff2")) { - if (_var1.sys().number() != _var2.sys().number()) + if (var1().sys().number() != var2().sys().number()) mooseError(this->type(), " does not support multiple nonlinear systems!"); } From 7bf27f2b3fe77e7c9b8823300a02f747445c8e0b Mon Sep 17 00:00:00 2001 From: Peter German Date: Wed, 27 Sep 2023 17:10:10 -0600 Subject: [PATCH 6/8] Add Alex's idea on setting the systems. (#25599) --- framework/include/problems/FEProblemBase.h | 13 +++++--- framework/src/fviks/FVInterfaceKernel.C | 20 +++++++---- framework/src/problems/FEProblemBase.C | 33 +++++++------------ .../fviks/FVConvectionCorrelationInterface.h | 3 ++ .../fviks/FVConvectionCorrelationInterface.C | 12 +++---- 5 files changed, 43 insertions(+), 38 deletions(-) diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h index 81325f4b6d00..7b058d71cf9a 100644 --- a/framework/include/problems/FEProblemBase.h +++ b/framework/include/problems/FEProblemBase.h @@ -851,13 +851,15 @@ class FEProblemBase : public SubProblem, public Restartable * @param name Name for the object to be created * @param parameters InputParameters for the object * @param threaded Whether or not to create n_threads copies of the object + * @param var_param_name The name of the parameter on the object which holds the primary variable. * @return A vector of shared_ptrs to the added objects */ template std::vector> addObject(const std::string & type, const std::string & name, InputParameters & parameters, - const bool threaded = true); + const bool threaded = true, + const std::string & var_param_name = "variable"); // Postprocessors ///// virtual void addPostprocessor(const std::string & pp_name, @@ -2295,7 +2297,9 @@ class FEProblemBase : public SubProblem, public Restartable * * This is needed due to header includes/forward declaration issues */ - void addObjectParamsHelper(InputParameters & params, const std::string & object_name); + void addObjectParamsHelper(InputParameters & params, + const std::string & object_name, + const std::string & var_param_name = "variable"); #ifdef LIBMESH_ENABLE_AMR Adaptivity _adaptivity; @@ -2585,12 +2589,13 @@ std::vector> FEProblemBase::addObject(const std::string & type, const std::string & name, InputParameters & parameters, - const bool threaded) + const bool threaded, + const std::string & var_param_name) { parallel_object_only(); // Add the _subproblem and _sys parameters depending on use_displaced_mesh - addObjectParamsHelper(parameters, name); + addObjectParamsHelper(parameters, name, var_param_name); const auto n_threads = threaded ? libMesh::n_threads() : 1; std::vector> objects(n_threads); diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index cf5c0ce883fc..553f9235f687 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -188,8 +188,11 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi) const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual()); - addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true); - if (_var1.sys().number() == _var2.sys().number()) + // If the two variables belong to two different nonlinear systems, we only contribute to the one + // which is being assembled right now + if (_var1.sys().number() == _subproblem.currentNlSysNum()) + addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true); + if (_var2.sys().number() == _subproblem.currentNlSysNum()) addResidual(_elem_is_one ? -r : r, _var2.number(), _elem_is_one ? true : false); } @@ -206,11 +209,14 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi) const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual(); - addResidualsAndJacobian(_assembly, - std::array{{_elem_is_one ? r : -r}}, - _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), - _var1.scalingFactor()); - if (_var1.sys().number() == _var2.sys().number()) + // If the two variables belong to two different nonlinear systems, we only contribute to the one + // which is being assembled right now + if (_var1.sys().number() == _subproblem.currentNlSysNum()) + addResidualsAndJacobian(_assembly, + std::array{{_elem_is_one ? r : -r}}, + _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), + _var1.scalingFactor()); + if (_var2.sys().number() == _subproblem.currentNlSysNum()) addResidualsAndJacobian(_assembly, std::array{{_elem_is_one ? -r : r}}, _elem_is_one ? _var2.dofIndicesNeighbor() : _var2.dofIndices(), diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index a7901dc9ce5d..23daba5757c6 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -2929,18 +2929,10 @@ FEProblemBase::addFVInterfaceKernel(const std::string & fv_ik_name, const std::string & name, InputParameters & parameters) { - /// Checking the nonlinear system number of variable1 - /// If variable1 and variable2 live on different systems, the interface kernel - /// is added to the system of variable1 alone. Another interface kernel needs to be - /// created with flipped variables for the system of variable2 - const auto nl_sys_num = - determineNonlinearSystem(parameters.varName("variable1", name), true).first - ? determineNonlinearSystem(parameters.varName("variable1", name), true).second - : (unsigned int)0; - - /// This is set because the Attribute system filters IKs based on this - parameters.set("_sys") = _nl[nl_sys_num].get(); - addObject(fv_ik_name, name, parameters); + /// We assume that variable1 and variable2 can live on different systems, in this case + /// the user needs to create two interface kernels with flipped variables and parameters + addObject( + fv_ik_name, name, parameters, /*threaded=*/true, /*variable_param_name=*/"variable1"); } // InterfaceKernels //// @@ -3567,21 +3559,21 @@ FEProblemBase::swapBackMaterialsNeighbor(THREAD_ID tid) } void -FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::string & object_name) +FEProblemBase::addObjectParamsHelper(InputParameters & parameters, + const std::string & object_name, + const std::string & var_param_name) { - const bool system_already_set = parameters.get("_sys"); const auto nl_sys_num = - parameters.isParamValid("variable") && - determineNonlinearSystem(parameters.varName("variable", object_name)).first - ? determineNonlinearSystem(parameters.varName("variable", object_name)).second + parameters.isParamValid(var_param_name) && + determineNonlinearSystem(parameters.varName(var_param_name, object_name)).first + ? determineNonlinearSystem(parameters.varName(var_param_name, object_name)).second : (unsigned int)0; if (_displaced_problem && parameters.have_parameter("use_displaced_mesh") && parameters.get("use_displaced_mesh")) { parameters.set("_subproblem") = _displaced_problem.get(); - if (!system_already_set) - parameters.set("_sys") = &_displaced_problem->nlSys(nl_sys_num); + parameters.set("_sys") = &_displaced_problem->nlSys(nl_sys_num); } else { @@ -3593,8 +3585,7 @@ FEProblemBase::addObjectParamsHelper(InputParameters & parameters, const std::st parameters.set("use_displaced_mesh") = false; parameters.set("_subproblem") = this; - if (!system_already_set) - parameters.set("_sys") = _nl[nl_sys_num].get(); + parameters.set("_sys") = _nl[nl_sys_num].get(); } } diff --git a/modules/navier_stokes/include/fviks/FVConvectionCorrelationInterface.h b/modules/navier_stokes/include/fviks/FVConvectionCorrelationInterface.h index 21c9b8f809c4..af5416a3d2a7 100644 --- a/modules/navier_stokes/include/fviks/FVConvectionCorrelationInterface.h +++ b/modules/navier_stokes/include/fviks/FVConvectionCorrelationInterface.h @@ -35,4 +35,7 @@ class FVConvectionCorrelationInterface : public FVInterfaceKernel /// libmesh object to find points in the mesh std::unique_ptr _pl; + + /// Boolean to see if variable1 is the fluid + const bool _var1_is_fluid; }; diff --git a/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C b/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C index 3e037e74c6e4..b0e6aadac3f0 100644 --- a/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C +++ b/modules/navier_stokes/src/fviks/FVConvectionCorrelationInterface.C @@ -37,7 +37,8 @@ FVConvectionCorrelationInterface::FVConvectionCorrelationInterface(const InputPa _htc(getFunctor("h")), _bulk_distance(getParam("bulk_distance")), _use_wall_cell(getParam("wall_cell_is_bulk")), - _pl(mesh().getPointLocator()) + _pl(mesh().getPointLocator()), + _var1_is_fluid("wraps_" + var1().name() == _temp_fluid.functorName()) { if (!_use_wall_cell && (_bulk_distance < 0)) mooseError( @@ -51,9 +52,8 @@ FVConvectionCorrelationInterface::computeQpResidual() // If variable1 is fluid and variable 1 is on elem or // if variable2 is fluid and variable 2 is on elem // the fluid element will be elem otherwise it is the neighbor - bool var1_is_fluid = ("wraps_" + var1().name() == _temp_fluid.functorName()); const Elem * elem_on_fluid_side = - (elemIsOne() && var1_is_fluid) || (!elemIsOne() && !var1_is_fluid) + (elemIsOne() && _var1_is_fluid) || (!elemIsOne() && !_var1_is_fluid) ? &_face_info->elem() : _face_info->neighborPtr(); @@ -77,11 +77,11 @@ FVConvectionCorrelationInterface::computeQpResidual() mooseAssert(bulk_elem, "The element at bulk_distance from the wall was not found in the mesh. " "Increase the number of ghost layers with the 'ghost_layers' parameter."); - mooseAssert((var1_is_fluid ? var1() : var2()).hasBlocks(bulk_elem->subdomain_id()), + mooseAssert((_var1_is_fluid ? var1() : var2()).hasBlocks(bulk_elem->subdomain_id()), "The fluid temperature is not defined at bulk_distance from the wall."); - const auto fluid_side = singleSidedFaceArg(var1_is_fluid ? var1() : var2(), _face_info); - const auto solid_side = singleSidedFaceArg(var1_is_fluid ? var2() : var1(), _face_info); + const auto fluid_side = singleSidedFaceArg(_var1_is_fluid ? var1() : var2(), _face_info); + const auto solid_side = singleSidedFaceArg(_var1_is_fluid ? var2() : var1(), _face_info); const auto bulk_elem_arg = makeElemArg(bulk_elem); From 275e912496ca21bd2cbc094ea57ac0ff2097c6d3 Mon Sep 17 00:00:00 2001 From: Peter German Date: Thu, 28 Sep 2023 13:20:04 -0600 Subject: [PATCH 7/8] Transform condition into assertion. (#25599) --- framework/src/fviks/FVInterfaceKernel.C | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index 553f9235f687..0d9cc4e8e907 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -190,9 +190,10 @@ FVInterfaceKernel::computeResidual(const FaceInfo & fi) // If the two variables belong to two different nonlinear systems, we only contribute to the one // which is being assembled right now - if (_var1.sys().number() == _subproblem.currentNlSysNum()) - addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true); - if (_var2.sys().number() == _subproblem.currentNlSysNum()) + mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(), + "The interface kernel should contribute to the system which variable1 belongs to!"); + addResidual(_elem_is_one ? r : -r, _var1.number(), _elem_is_one ? false : true); + if (_var1.sys().number() == _var2.sys().number()) addResidual(_elem_is_one ? -r : r, _var2.number(), _elem_is_one ? true : false); } @@ -211,12 +212,13 @@ FVInterfaceKernel::computeJacobian(const FaceInfo & fi) // If the two variables belong to two different nonlinear systems, we only contribute to the one // which is being assembled right now - if (_var1.sys().number() == _subproblem.currentNlSysNum()) - addResidualsAndJacobian(_assembly, - std::array{{_elem_is_one ? r : -r}}, - _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), - _var1.scalingFactor()); - if (_var2.sys().number() == _subproblem.currentNlSysNum()) + mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(), + "The interface kernel should contribute to the system which variable1 belongs to!"); + addResidualsAndJacobian(_assembly, + std::array{{_elem_is_one ? r : -r}}, + _elem_is_one ? _var1.dofIndices() : _var1.dofIndicesNeighbor(), + _var1.scalingFactor()); + if (_var1.sys().number() == _var2.sys().number()) addResidualsAndJacobian(_assembly, std::array{{_elem_is_one ? -r : r}}, _elem_is_one ? _var2.dofIndicesNeighbor() : _var2.dofIndices(), From e8e0c98e344983c3b9e3aab218353f7a90a5166b Mon Sep 17 00:00:00 2001 From: Peter German Date: Fri, 29 Sep 2023 10:02:47 -0600 Subject: [PATCH 8/8] Add documentation for instructions on how to set kernels up using multiple nonlinear systems. (#25599) --- .../doc/content/source/fviks/FVDiffusionInterface.md | 4 ++++ .../doc/content/syntax/FVInterfaceKernels/index.md | 12 ++++++++++++ .../source/fviks/FVConvectionCorrelationInterface.md | 4 ++++ 3 files changed, 20 insertions(+) diff --git a/framework/doc/content/source/fviks/FVDiffusionInterface.md b/framework/doc/content/source/fviks/FVDiffusionInterface.md index 5d3140728a04..f035e49afa7d 100644 --- a/framework/doc/content/source/fviks/FVDiffusionInterface.md +++ b/framework/doc/content/source/fviks/FVDiffusionInterface.md @@ -5,6 +5,10 @@ The diffusive flux is obtained from a two point gradient, and the diffusivity is interpolated to the interface. +!alert note +This kernel supports interfaces between variables which belong to different nonlinear systems. +For instructions on how to set these cases up, visit the [FVInterfaceKernels syntax page](syntax/FVInterfaceKernels/index.md). + ## Example input file syntax In this example, two diffusion problems with a source terms are solved on each side diff --git a/framework/doc/content/syntax/FVInterfaceKernels/index.md b/framework/doc/content/syntax/FVInterfaceKernels/index.md index 552f85568bcf..cbf17b215de7 100644 --- a/framework/doc/content/syntax/FVInterfaceKernels/index.md +++ b/framework/doc/content/syntax/FVInterfaceKernels/index.md @@ -28,6 +28,7 @@ side. So making use of the `subdomain` parameters, we provide a protected method called `elemIsOne()` that returns a boolean indicating whether the `FaceInfo::elem` side of the interface corresponds to the `subdomain1` side of the interface. This allows the developer to write code like the following: + ``` FVFooInterface::FVFooInterface(const InputParameters & params) : FVInterfaceKernel(params), @@ -47,5 +48,16 @@ FVFooInterface::computeQpResidual() /// Code that uses coef_elem and coef_neighbor } ``` + and have confidence that they have good data in `coef_elem` and `coef_neighbor` and have clarity about what is happening in their code. + +!alert! note +When using an FVInterfaceKernel which connects variables that belong to different nonlinear systems, +create two kernels with flipped variable and material property parameters. The reason behind this +is that the interface kernel will only contribute to the system which `variable1` belongs to. +For an example, see: + +!listing /test/tests/fviks/diffusion/multisystem.i + +!alert-end! diff --git a/modules/navier_stokes/doc/content/source/fviks/FVConvectionCorrelationInterface.md b/modules/navier_stokes/doc/content/source/fviks/FVConvectionCorrelationInterface.md index 37e10570846c..50071e380b2b 100644 --- a/modules/navier_stokes/doc/content/source/fviks/FVConvectionCorrelationInterface.md +++ b/modules/navier_stokes/doc/content/source/fviks/FVConvectionCorrelationInterface.md @@ -10,6 +10,10 @@ with $q_s$ the surface convective heat flux, $h_{correlation}$ the heat transfer defined by the correlation as a material property and $T_{solid/fluid}$ the temperature of the adjacent solid and fluid. +!alert note +This kernel supports interfaces between variables which belong to different nonlinear systems. +For instructions on how to set these cases up, visit the [FVInterfaceKernels syntax page](syntax/FVInterfaceKernels/index.md). + ## Example input file syntax In this example, a cold fluid is flowing next to a centrally-heated solid region. The heat diffuses