225 bool is_dirichlet,
bool is_inhomog,
227 Real bscalar)
noexcept
229 const Real dxi = dxinv[0];
230 const Real dyi = dxinv[1];
234 for (
int j = lo.y; j <= hi.y; ++j) {
236 for (
int i = lo.x; i <= hi.x; ++i) {
237 if (flag(i,j,0).isRegular())
239 Ax(i,j,0,0) += bscalar*(dxi*(fx(i+1,j ,0,0) - fx(i,j,0,0))
240 + dyi*(fy(i ,j+1,0,0) - fy(i,j,0,0)));
241 Ax(i,j,0,1) += bscalar*(dxi*(fx(i+1,j ,0,1) - fx(i,j,0,1))
242 + dyi*(fy(i ,j+1,0,1) - fy(i,j,0,1)));
244 else if (flag(i,j,0).isSingleValued())
246 Real fxm_0 = fx(i,j,0,0);
247 Real fxm_1 = fx(i,j,0,1);
248 if (apx(i,j,0) > Real(0.0) && apx(i,j,0) < Real(1.0)) {
249 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
250 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
251 fxm_0 = (Real(1.0)-fracy)*fxm_0 + fracy*fx(i,jj,0,0);
252 fxm_1 = (Real(1.0)-fracy)*fxm_1 + fracy*fx(i,jj,0,1);
255 Real fxp_0 = fx(i+1,j,0,0);
256 Real fxp_1 = fx(i+1,j,0,1);
257 if (apx(i+1,j,0) > Real(0.0) && apx(i+1,j,0) < Real(1.0)) {
258 int jj = j + (fcx(i+1,j,0,0) >= Real(0.0) ? 1 : -1);
259 Real fracy = (ccm(i,jj,0) || ccm(i+1,jj,0)) ? std::abs(fcx(i+1,j,0,0)) : Real(0.0);
260 fxp_0 = (Real(1.0)-fracy)*fxp_0 + fracy*fx(i+1,jj,0,0);
261 fxp_1 = (Real(1.0)-fracy)*fxp_1 + fracy*fx(i+1,jj,0,1);
264 Real fym_0 = fy(i,j,0,0);
265 Real fym_1 = fy(i,j,0,1);
266 if (apy(i,j,0) > Real(0.0) && apy(i,j,0) < Real(1.0)) {
267 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
268 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
269 fym_0 = (Real(1.0)-fracx)*fym_0 + fracx*fy(ii,j,0,0);
270 fym_1 = (Real(1.0)-fracx)*fym_1 + fracx*fy(ii,j,0,1);
273 Real fyp_0 = fy(i,j+1,0,0);
274 Real fyp_1 = fy(i,j+1,0,1);
275 if (apy(i,j+1,0) > Real(0.0) && apy(i,j+1,0) < Real(1.0)) {
276 int ii = i + (fcy(i,j+1,0,0) >= Real(0.0) ? 1 : -1);
277 Real fracx = (ccm(ii,j,0) || ccm(ii,j+1,0)) ? std::abs(fcy(i,j+1,0,0)) : Real(0.0);
278 fyp_0 = (Real(1.0)-fracx)*fyp_0 + fracx*fy(ii,j+1,0,0);
279 fyp_1 = (Real(1.0)-fracx)*fyp_1 + fracx*fy(ii,j+1,0,1);
282 Real kappa = vol(i,j,0);
283 Real feb_0 = Real(0.0);
284 Real feb_1 = Real(0.0);
287 Real dapx = apx(i+1,j,0)-apx(i,j,0);
288 Real dapy = apy(i,j+1,0)-apy(i,j,0);
289 Real anorminv = Real(1.0)/std::sqrt(dapx*dapx+dapy*dapy);
290 Real anrmx = -dapx * anorminv;
291 Real anrmy = -dapy * anorminv;
293 Real velb_0 = 0, velb_1 = 0;
296 velb_0 = velb(i,j,0,0);
297 velb_1 = velb(i,j,0,1);
300 Real dx_eb =
amrex::max(Real(0.3), (kappa*kappa-Real(0.25))/(Real(2.)*kappa));
301 Real dg = dx_eb /
amrex::max(std::abs(anrmx),std::abs(anrmy));
302 Real dginv = Real(1.0)/dg;
303 Real gx = bc(i,j,0,0) - dg*anrmx;
304 Real gy = bc(i,j,0,1) - dg*anrmy;
305 int isx = (anrmx > Real(0.0)) ? 1 : -1;
306 int isy = (anrmy > Real(0.0)) ? 1 : -1;
313 Real oneggg = Real(1.0)+gx+gy+gxy;
315 Real velg = oneggg * vel(i ,j ,0,0)
316 + (-gy - gxy) * vel(i ,jj,0,0)
317 + (-gx - gxy) * vel(ii,j ,0,0)
318 + gxy * vel(ii,jj,0,0);
319 Real dudn = (velb_0-velg) * dginv;
321 velg = oneggg * vel(i ,j ,0,1)
322 + (-gy - gxy) * vel(i ,jj,0,1)
323 + (-gx - gxy) * vel(ii,j ,0,1)
324 + gxy * vel(ii,jj,0,1);
325 Real dvdn = (velb_1-velg) * dginv;
328 Real dudx = dudn * anrmx;
329 Real dudy = dudn * anrmy;
330 Real dvdx = dvdn * anrmx;
331 Real dvdy = dvdn * anrmy;
332 Real divu = dudx + dvdy;
333 Real xi = kapb(i,j,0);
334 Real mu = etab(i,j,0);
335 Real tautmp = (xi-(2./3.)*mu)*divu;
337 Real tauxx = mu*dudx + tautmp;
338 Real tauyx = mu*dvdx;
339 Real tauyy = mu*dvdy + tautmp;
340 Real tauxy = mu*dudy;
342 feb_0 = dxi*(dapx*tauxx + dapy*tauyx);
343 feb_1 = dxi*(dapx*tauxy + dapy*tauyy);
346 Real volinv = bscalar / kappa;
347 Ax(i,j,0,0) += volinv * (dxi*(apx(i,j,0)*fxm_0-apx(i+1,j,0)*fxp_0)
348 +dyi*(apy(i,j,0)*fym_0-apy(i,j+1,0)*fyp_0)
350 Ax(i,j,0,1) += volinv * (dxi*(apx(i,j,0)*fxm_1-apx(i+1,j,0)*fxp_1)
351 +dyi*(apy(i,j,0)*fym_1-apy(i,j+1,0)*fyp_1)
656 const Real dxi = dxinv[0];
657 const Real dyi = dxinv[1];
662 for (
int j = lo.y; j <= hi.y; ++j) {
664 for (
int i = lo.x; i <= hi.x; ++i) {
665 if (apx(i,j,0) == Real(0.0))
667 fx(i,j,0,0) = Real(0.0);
668 fx(i,j,0,1) = Real(0.0);
669 fx(i,j,0,2) = Real(0.0);
670 fx(i,j,0,3) = Real(0.0);
674 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
675 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
676 if (apx(i,j,0) < Real(1.0)) {
677 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
678 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
679 dudx = (Real(1.0)-fracy)*dudx
680 + fracy *(vel(i,jj,0,0) - vel(i-1,jj,0,0))*dxi;
681 dvdx = (Real(1.0)-fracy)*dvdx
682 + fracy *(vel(i,jj,0,1) - vel(i-1,jj,0,1))*dxi;
685 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
686 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
687 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
688 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
692 whi,wlo,jhip,jhim,jlop,jlom);
694 whi,wlo,jhip,jhim,jlop,jlom);
712 const Real dxi = dxinv[0];
713 const Real dyi = dxinv[1];
718 for (
int j = lo.y; j <= hi.y; ++j) {
720 for (
int i = lo.x; i <= hi.x; ++i) {
721 if (apy(i,j,0) == 0.) {
722 fy(i,j,0,0) = Real(0.0);
723 fy(i,j,0,1) = Real(0.0);
724 fy(i,j,0,2) = Real(0.0);
725 fy(i,j,0,3) = Real(0.0);
729 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
730 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
731 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
732 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
736 whi,wlo,ihip,ihim,ilop,ilom);
738 whi,wlo,ihip,ihim,ilop,ilom);
740 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
741 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
742 if (apy(i,j,0) < Real(1.0)) {
743 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
744 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
745 dudy = (Real(1.0)-fracx)*dudy
746 + fracx *(vel(ii,j,0,0) - vel(ii,j-1,0,0))*dyi;
747 dvdy = (Real(1.0)-fracx)*dvdy
748 + fracx *(vel(ii,j,0,1) - vel(ii,j-1,0,1))*dyi;
771 0,AMREX_SPACEDIM>
const& bct,
772 Dim3 const& dlo,
Dim3 const& dhi)
noexcept
774 const Real dxi = dxinv[0];
775 const Real dyi = dxinv[1];
780 for (
int j = lo.y; j <= hi.y; ++j) {
781 for (
int i = lo.x; i <= hi.x; ++i) {
782 if (apx(i,j,0) == Real(0.0))
784 fx(i,j,0,0) = Real(0.0);
785 fx(i,j,0,1) = Real(0.0);
786 fx(i,j,0,2) = Real(0.0);
787 fx(i,j,0,3) = Real(0.0);
791 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
792 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
793 if (apx(i,j,0) < Real(1.0)) {
794 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
795 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
796 dudx = (Real(1.0)-fracy)*dudx
797 + fracy *(vel(i,jj,0,0) - vel(i-1,jj,0,0))*dxi;
798 dvdx = (Real(1.0)-fracy)*dvdx
799 + fracy *(vel(i,jj,0,1) - vel(i-1,jj,0,1))*dxi;
802 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
803 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
804 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
805 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
809 bvxlo,bvxhi,bct,dlo,dhi,
810 whi,wlo,jhip,jhim,jlop,jlom);
812 bvxlo,bvxhi,bct,dlo,dhi,
813 whi,wlo,jhip,jhim,jlop,jlom);
834 0,AMREX_SPACEDIM>
const& bct,
835 Dim3 const& dlo,
Dim3 const& dhi)
noexcept
837 const Real dxi = dxinv[0];
838 const Real dyi = dxinv[1];
843 for (
int j = lo.y; j <= hi.y; ++j) {
844 for (
int i = lo.x; i <= hi.x; ++i) {
845 if (apy(i,j,0) == 0.) {
846 fy(i,j,0,0) = Real(0.0);
847 fy(i,j,0,1) = Real(0.0);
848 fy(i,j,0,2) = Real(0.0);
849 fy(i,j,0,3) = Real(0.0);
853 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
854 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
855 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
856 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
860 bvylo,bvyhi,bct,dlo,dhi,
861 whi,wlo,ihip,ihim,ilop,ilom);
863 bvylo,bvyhi,bct,dlo,dhi,
864 whi,wlo,ihip,ihim,ilop,ilom);
866 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
867 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
868 if (apy(i,j,0) < Real(1.0)) {
869 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
870 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
871 dudy = (Real(1.0)-fracx)*dudy
872 + fracx *(vel(ii,j,0,0) - vel(ii,j-1,0,0))*dyi;
873 dvdy = (Real(1.0)-fracx)*dvdy
874 + fracx *(vel(ii,j,0,1) - vel(ii,j-1,0,1))*dyi;