From bcad0ef5b085cd7cca0327230a4c91b65f56d6b6 Mon Sep 17 00:00:00 2001 From: TravisMitchell Date: Wed, 20 Dec 2023 12:28:32 +1000 Subject: [PATCH 1/4] + initial fix to remove nan. todo-new field for solid, remove edge race condition --- models/multiphase/d2q9_pf_velocity/Dynamics.R | 15 +- .../multiphase/d2q9_pf_velocity/Dynamics.c.Rt | 941 ++++++++---------- models/multiphase/d2q9_pf_velocity/conf.mk | 3 +- 3 files changed, 449 insertions(+), 510 deletions(-) diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.R b/models/multiphase/d2q9_pf_velocity/Dynamics.R index d2319bb53..789c36aa9 100644 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.R +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.R @@ -86,20 +86,20 @@ if (Options$RT) { AddStage("BaseInit" , "Init_distributions" , save=Fields$group %in% c("g","h","Vel")) # iteration - AddStage("BaseIter" , "calcHydroIter" , save=Fields$group %in% c("g","h","Vel","nw") , load=DensityAll$group %in% c("g","h","Vel","nw")) # TODO: is nw needed here? + AddStage("BaseIter" , "calcHydroIter" , save=Fields$group %in% c("g","h","Vel","nw") , load=DensityAll$group %in% c("PF","g","h","Vel","nw")) # TODO: is nw needed here? AddStage("PhaseIter" , "calcPhaseFIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("g","h","Vel","nw")) - AddStage("WallIter" , "calcWallPhaseIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("nw")) + AddStage("WallIter" , "calcWallPhaseIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("nw","PF")) } AddAction("Iteration", c("BaseIter", "PhaseIter","WallIter")) AddAction("Init" , c("PhaseInit","WallInit", "WallIter","BaseInit")) # Outputs: -AddQuantity(name="Rho", unit="kg/m3") -AddQuantity(name="PhaseField",unit="1") -AddQuantity(name="U", unit="m/s",vector=T) -AddQuantity(name="NormalizedPressure", unit="Pa") -AddQuantity(name="Pressure", unit="Pa") +AddQuantity(name="Rho", unit="kg/m3") +AddQuantity(name="PhaseField", unit="1") +AddQuantity(name="U", unit="m/s",vector=T) +AddQuantity(name="NormalizedPressure", unit="1") +AddQuantity(name="Pressure", unit="Pa") AddQuantity(name="Normal", unit="1", vector=T) # Initialisation States @@ -209,5 +209,4 @@ if (Options$Outflow) { } AddNodeType(name="Solid", group="BOUNDARY") AddNodeType(name="Wall", group="BOUNDARY") -AddNodeType(name="BGK", group="COLLISION") AddNodeType(name="MRT", group="COLLISION") diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt index 96e0cd550..3484c8e2d 100755 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt @@ -114,8 +114,27 @@ CudaDeviceFunction vector_t calcGradPhi(){ gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,-1) - PhaseF(-1,-1) )/6.0; gradPhi.y = 0.0; } else if ((NodeType & NODE_BOUNDARY) == NODE_Wall){ - gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0; - gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0; + // compute gradient with one-sided diffs but first check if neighbours are nans + bool nan1 = isnan(PhaseF(1,0)); + bool nan2 = isnan(PhaseF(-1,0)); + bool nan3 = isnan(PhaseF(0,1)); + bool nan4 = isnan(PhaseF(0,-1)); + + if ((nan1 && nan2) || (!nan1 && !nan2)){ + gradPhi.x = 0.0; + } else if (nan1){ + gradPhi.x = PhaseF(-1,0) - PhaseF(0,0); + } else { + gradPhi.x = PhaseF(1,0) - PhaseF(0,0); + } + + if ((nan3 && nan4) || (!nan3 && !nan4)){ + gradPhi.y = 0.0; + } else if (nan3){ + gradPhi.y = PhaseF(0,-1) - PhaseF(0,0); + } else { + gradPhi.y = PhaseF(0,1) - PhaseF(0,0); + } } else if (NodeType & NODE_BOUNDARY){ gradPhi.x = 0.0; gradPhi.y = 0.0; @@ -125,8 +144,27 @@ CudaDeviceFunction vector_t calcGradPhi(){ } #else if ((NodeType & NODE_BOUNDARY) == NODE_Wall){ - gradPhi.x = (PhaseF(1,0) - PhaseF(-1,0))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(1,-1) - PhaseF(-1,1))/12.0; - gradPhi.y = (PhaseF(0,1) - PhaseF(0,-1))/3.0 + (PhaseF(1,1) - PhaseF(-1,-1) + PhaseF(-1,1) - PhaseF(1,-1))/12.0; + // compute gradient with one-sided diffs but first check if neighbours are nans + bool nan1 = isnan(PhaseF(1,0)); + bool nan2 = isnan(PhaseF(-1,0)); + bool nan3 = isnan(PhaseF(0,1)); + bool nan4 = isnan(PhaseF(0,-1)); + + if ((nan1 && nan2) || (!nan1 && !nan2)){ + gradPhi.x = 0.0; + } else if (nan1){ + gradPhi.x = PhaseF(-1,0) - PhaseF(0,0); + } else { + gradPhi.x = PhaseF(1,0) - PhaseF(0,0); + } + + if ((nan3 && nan4) || (!nan3 && !nan4)){ + gradPhi.y = 0.0; + } else if (nan3){ + gradPhi.y = PhaseF(0,-1) - PhaseF(0,0); + } else { + gradPhi.y = PhaseF(0,1) - PhaseF(0,0); + } } else if (NodeType & NODE_BOUNDARY){ gradPhi.x = 0.0; gradPhi.y = 0.0; @@ -285,6 +323,7 @@ CudaDeviceFunction real_t calcRho(real_t myPhaseF) CudaDeviceFunction void Init_phase() { PhaseF = PhaseField_init;// read initial PhaseField distribution from config.xml file + // Special initialisation cases if ( Period > 0 ) { if ( Wave > 0) { real_t InterfacePoint = MidPoint + Perturbation*cos(2.0*PI*X/Period + PI); @@ -304,21 +343,23 @@ CudaDeviceFunction void Init_phase() { PhaseF = 0.5 * (PhaseField_h + PhaseField_l) - 0.5 * (PhaseField_h - PhaseField_l) * BubbleType * tanh(2.0*(Ri - Radius)/W); } - if ((NodeType & NODE_BOUNDARY) == NODE_Wall) PhaseF = -999; + + // Wall initialisation to calculate normal in second stage of init action + if ((NodeType & NODE_BOUNDARY) == NODE_Wall) PhaseF = NAN; } CudaDeviceFunction void Init_wallNorm(){ if ((NodeType & NODE_BOUNDARY) == NODE_Wall) { int solidFlags[9], i; - solidFlags[0] = PhaseF(0,0) /-998; - solidFlags[1] = PhaseF(1,0) /-998; - solidFlags[2] = PhaseF(0,1) /-998; - solidFlags[3] = PhaseF(-1,0) /-998; - solidFlags[4] = PhaseF(0,-1) /-998; - solidFlags[5] = PhaseF(1,1) /-998; - solidFlags[6] = PhaseF(-1,1 )/-998; - solidFlags[7] = PhaseF(-1,-1)/-998; - solidFlags[8] = PhaseF(1,-1) /-998; + solidFlags[0] = 1; + solidFlags[1] = isnan(PhaseF(1,0)) ? 1 : 0; + solidFlags[2] = isnan(PhaseF(0,1)) ? 1 : 0; + solidFlags[3] = isnan(PhaseF(-1,0)) ? 1 : 0; + solidFlags[4] = isnan(PhaseF(0,-1)) ? 1 : 0; + solidFlags[5] = isnan(PhaseF(1,1)) ? 1 : 0; + solidFlags[6] = isnan(PhaseF(-1,1)) ? 1 : 0; + solidFlags[7] = isnan(PhaseF(-1,-1)) ? 1 : 0; + solidFlags[8] = isnan(PhaseF(1,-1)) ? 1 : 0; real_t myNorm[2]={0.0,0.0}; for (i = 0 ; i < 9; i++){ @@ -350,29 +391,25 @@ CudaDeviceFunction void Init_wallNorm(){ CudaDeviceFunction void Init_distributions(){ - // Find Gradients and normals: real_t myPhaseF = PhaseF(0,0); - vector_t gradPhi = calcGradPhi(); real_t nx, ny; - // gradPhi.z stores the magnitude of the vector - nx = gradPhi.x/gradPhi.z; + nx = gradPhi.x/gradPhi.z; //gradPhi.z stores the magnitude of the vector ny = gradPhi.y/gradPhi.z; - // Define Equilibrium, then initialise all da things U = VelocityX; V = VelocityY; real_t mag = U*U + V*V; - real_t Gamma[9], F_phi[9]; + real_t Gamma, F_phi; real_t pfavg = 0.5*(PhaseField_l+PhaseField_h); real_t tmp1 = (1.0 - 4.0*(myPhaseF - pfavg)*(myPhaseF - pfavg))/W; for (int i=0; i< 9; i++){ - Gamma[i] = calcGamma(i, U, V, mag); - F_phi[i] = calcF_phi(i, tmp1, nx, ny); - h[i] = myPhaseF * Gamma[i] - 0.5*F_phi[i]; // heq - g[i] = Gamma[i] - wf[i]; //geq + Gamma = calcGamma(i, U, V, mag); + F_phi = calcF_phi(i, tmp1, nx, ny); + h[i] = myPhaseF * Gamma - 0.5*F_phi; // heq + g[i] = Gamma - wf[i]; //geq } #ifdef OPTIONS_RT @@ -386,7 +423,6 @@ CudaDeviceFunction void Init_distributions(){ C(h_old, h) } ?> } #endif - PhaseF = ; } // ITERATION: @@ -425,26 +461,25 @@ CudaDeviceFunction void BC_Switcher () } CudaDeviceFunction void calcHydroIter() { - real_t tmpPhase = 0.0; if ((NodeType & NODE_ADDITIONALS) == NODE_Smoothing){ Init_distributions(); - } - else{ + } else { PhaseF = PhaseF(0,0); - tmpPhase = PhaseF; + real_t myPhaseF = PhaseF; - real_t rho = calcRho(myPhaseF); log_data_from_special_nodes(); + #ifdef OPTIONS_debug + real_t rho = calcRho(myPhaseF); AddToMomentumX(U*rho); AddToMomentumY(V*rho); #endif - BC_Switcher(); + BC_Switcher(); - if (PhaseF <= 0.5) { -// AddToBubbleVelocityX( (1 - PhaseF)*U ); -// AddToBubbleVelocityY( (1 - PhaseF)*V ); + if (PhaseF <= 0.5) { + AddToBubbleVelocityX( (1 - PhaseF)*U ); + AddToBubbleVelocityY( (1 - PhaseF)*V ); AddToBubbleLocationY( (1 - PhaseF)*Y ); AddToSumPhiGas( (1 - PhaseF) ); } @@ -455,15 +490,10 @@ CudaDeviceFunction void calcHydroIter() { #else if (NodeType & NODE_MRT) { CollisionMRT();} #endif - if (tmpPhase <= 0.5) { - AddToBubbleVelocityX( (1 - tmpPhase)*U ); - AddToBubbleVelocityY( (1 - tmpPhase)*V ); - } - myPhaseF = PhaseF(0,0); - rho = calcRho(myPhaseF); - #ifdef OPTIONS_debug + myPhaseF = PhaseF(0,0); + rho = calcRho(myPhaseF); AddToMomentumX_afterCol(U*rho); AddToMomentumY_afterCol(V*rho); #endif @@ -505,465 +535,382 @@ CudaDeviceFunction void calcPhaseFIter(){ } CudaDeviceFunction void calcWallPhaseIter(){ -// Contact angle is defined with respect to the high density fluid + // Contact angle is defined with respect to the high density fluid PhaseF = PhaseF(0,0); if ((NodeType & NODE_BOUNDARY) == NODE_Wall) { if (nw_x == 0 && nw_y == 0) { - PhaseF = 0.0; // TODO: why not simply PhaseF = PhaseF(0,0); + PhaseF = 0.0; + printf("EDGE CASE: Solid found with no normal at (%.1f %.1f)\n", X, Y); return; } - real_t myA, myH, myPhase; - + real_t myA, myPhase; myPhase = PhaseF_dyn(nw_x, nw_y); - myH = sqrt(nw_x*nw_x + nw_y*nw_y ); - - if (fabs(radAngle - PI/2) < 1e-4) { PhaseF = myPhase; } - else { - myA = 1.0 - 0.5*myH * (4.0/W) * cos( radAngle ); + if (isnan(myPhase)) { + printf("EDGE CASE: Solid normal pointing into solid node at (%.1f %.1f)\n", X, Y); + PhaseF=0.0; + return; + } + if (fabs(radAngle - PI/2) < 1e-4) { + PhaseF = myPhase; + } else { + myA = 1.0 - 0.5 * sqrt(nw_x*nw_x + nw_y*nw_y ) * (4.0/W) * cos( radAngle ); PhaseF = (myA - sqrt( myA*myA - 4.0 * (myA-1.0)*myPhase))/(myA-1.0) - myPhase; } - } + } + if (isnan(PhaseF)) printf("NaN at (%.1f %.1f)\n", X, Y); } #ifdef OPTIONS_CM - CudaDeviceFunction void relax_CM_hydro(real_t f_in[9], real_t tau, vector_t u) - { - real_t uxuy = u.x*u.y; - real_t ux2 = u.x*u.x; - real_t uy2 = u.y*u.y; - - real_t s_v = 1./tau; - real_t s_b = omega_bulk; - - real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; - - real_t temp[9]; real_t cm_eq[9]; - for (int i = 0; i < 9; i++) { - temp[i] = f_in[i];} - - //raw moments from density-probability functions - //[m00, m10, m01, m20, m02, m11, m21, m12, m22] - f_in[0] = m00; - f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; - f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; - f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; - f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; - f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; - f_in[6] = temp[5] + temp[6] - temp[7] - temp[8]; - f_in[7] = temp[5] - temp[6] - temp[7] + temp[8]; - f_in[8] = temp[5] + temp[6] + temp[7] + temp[8]; - - //central moments from raw moments - temp[0] = f_in[0]; - temp[1] = -f_in[0]*u.x + f_in[1]; - temp[2] = -f_in[0]*u.y + f_in[2]; - temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3]; - temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4]; - temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; - temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6]; - temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7]; - temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8]; - - //collision in central moments space - //calculate equilibrium distributions in cm space - cm_eq[0] = m00; - cm_eq[1] = -u.x*(m00 - 1.0); - cm_eq[2] = -u.y*(m00 - 1.0); - cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2; - cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2; - cm_eq[5] = uxuy*(m00 - 1.0); - cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.); - cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.); - cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2; - //collide eq: -S*(cm - cm_eq) - f_in[0] = cm_eq[0] - temp[0]; - f_in[1] = cm_eq[1] - temp[1]; - f_in[2] = cm_eq[2] - temp[2]; - f_in[3] = 0.5*(cm_eq[3] - temp[3])*(s_b + s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b - s_v); - f_in[4] = 0.5*(cm_eq[3] - temp[3])*(s_b - s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b + s_v); - f_in[5] = s_v*(cm_eq[5] - temp[5]); - f_in[6] = cm_eq[6] - temp[6]; - f_in[7] = cm_eq[7] - temp[7]; - f_in[8] = cm_eq[8] - temp[8]; - - //back to raw moments - temp[0] = f_in[0]; - temp[1] = f_in[0]*u.x + f_in[1]; - temp[2] = f_in[0]*u.y + f_in[2]; - temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3]; - temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4]; - temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; - temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6]; - temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7]; - temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8]; - - //back to density-probability functions - f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; - f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2; - f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2; - f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2; - f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2; - f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4; - f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4; - f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4; - f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4; - } +CudaDeviceFunction void relax_CM_hydro(real_t f_in[9], real_t tau, vector_t u) +{ + real_t uxuy = u.x*u.y; + real_t ux2 = u.x*u.x; + real_t uy2 = u.y*u.y; + + real_t s_v = 1./tau; + real_t s_b = omega_bulk; + + real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; + + real_t temp[9]; real_t cm_eq[9]; + for (int i = 0; i < 9; i++) { + temp[i] = f_in[i];} + + //raw moments from density-probability functions + //[m00, m10, m01, m20, m02, m11, m21, m12, m22] + f_in[0] = m00; + f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; + f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; + f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; + f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; + f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; + f_in[6] = temp[5] + temp[6] - temp[7] - temp[8]; + f_in[7] = temp[5] - temp[6] - temp[7] + temp[8]; + f_in[8] = temp[5] + temp[6] + temp[7] + temp[8]; + + //central moments from raw moments + temp[0] = f_in[0]; + temp[1] = -f_in[0]*u.x + f_in[1]; + temp[2] = -f_in[0]*u.y + f_in[2]; + temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3]; + temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4]; + temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; + temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6]; + temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7]; + temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8]; + + //collision in central moments space + //calculate equilibrium distributions in cm space + cm_eq[0] = m00; + cm_eq[1] = -u.x*(m00 - 1.0); + cm_eq[2] = -u.y*(m00 - 1.0); + cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2; + cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2; + cm_eq[5] = uxuy*(m00 - 1.0); + cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.); + cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.); + cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2; + //collide eq: -S*(cm - cm_eq) + f_in[0] = cm_eq[0] - temp[0]; + f_in[1] = cm_eq[1] - temp[1]; + f_in[2] = cm_eq[2] - temp[2]; + f_in[3] = 0.5*(cm_eq[3] - temp[3])*(s_b + s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b - s_v); + f_in[4] = 0.5*(cm_eq[3] - temp[3])*(s_b - s_v) + 0.5*(cm_eq[4] - temp[4])*(s_b + s_v); + f_in[5] = s_v*(cm_eq[5] - temp[5]); + f_in[6] = cm_eq[6] - temp[6]; + f_in[7] = cm_eq[7] - temp[7]; + f_in[8] = cm_eq[8] - temp[8]; + + //back to raw moments + temp[0] = f_in[0]; + temp[1] = f_in[0]*u.x + f_in[1]; + temp[2] = f_in[0]*u.y + f_in[2]; + temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3]; + temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4]; + temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; + temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6]; + temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7]; + temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8]; + + //back to density-probability functions + f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; + f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2; + f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2; + f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2; + f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2; + f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4; + f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4; + f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4; + f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4; +} - CudaDeviceFunction void relax_and_collide_CM_hydro(real_t f_in[9], real_t tau, vector_t u, vector_t Fhydro, real_t rho) - { - real_t uxuy = u.x*u.y; - real_t ux2 = u.x*u.x; - real_t uy2 = u.y*u.y; - - real_t s_v = 1./tau; - real_t s_b = omega_bulk; - real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; - real_t temp[9]; real_t cm_eq[9]; - for (int i = 0; i < 9; i++) { - temp[i] = f_in[i];} - //raw moments from density-probability functions - //[m00, m10, m01, m20, m02, m11, m21, m12, m22] - f_in[0] = m00; - f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; - f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; - f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; - f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; - f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; - //central moments from raw moments - temp[3] = f_in[0]*ux2 - 2.0*f_in[1]*u.x + f_in[3]; - temp[4] = f_in[0]*uy2 - 2.0*f_in[2]*u.y + f_in[4]; - temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; - //collision in central moments space - //calculate equilibrium distributions in cm space - cm_eq[0] = m00; - cm_eq[1] = -u.x*(m00 - 1.0); - cm_eq[2] = -u.y*(m00 - 1.0); - cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2; - cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2; - cm_eq[5] = uxuy*(m00 - 1.0); - cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.); - cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.); - cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2; - - //calculate forces in cm space - //He et al. scheme: get_continuous_Maxwellian_DF(dzeta=dzeta, psi=1, u=(ux, uy)) - velocity term - //collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*force_in_cm_space // SOI - f_in[0] = cm_eq[0]; - f_in[1] = 0.5*Fhydro.x/rho + cm_eq[1]; - f_in[2] = 0.5*Fhydro.y/rho + cm_eq[2]; - f_in[3] = 0.5*cm_eq[3]*(s_b + s_v) + 0.5*cm_eq[4]*(s_b - s_v) - temp[3]*(0.5*s_b + 0.5*s_v - 1.0) - 0.5*temp[4]*(s_b - s_v); - f_in[4] = 0.5*cm_eq[3]*(s_b - s_v) + 0.5*cm_eq[4]*(s_b + s_v) - 0.5*temp[3]*(s_b - s_v) - temp[4]*(0.5*s_b + 0.5*s_v - 1.0); - f_in[5] = cm_eq[5]*s_v - temp[5]*(s_v - 1.0); - f_in[6] = 1./6.*Fhydro.y/rho + cm_eq[6]; - f_in[7] = 1./6.*Fhydro.x/rho + cm_eq[7]; - f_in[8] = cm_eq[8]; - - //back to raw moments - temp[0] = f_in[0]; - temp[1] = f_in[0]*u.x + f_in[1]; - temp[2] = f_in[0]*u.y + f_in[2]; - temp[3] = f_in[0]*ux2 + 2.0*f_in[1]*u.x + f_in[3]; - temp[4] = f_in[0]*uy2 + 2.0*f_in[2]*u.y + f_in[4]; - temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; - temp[6] = f_in[0]*ux2*u.y + 2.0*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2.0*f_in[5]*u.x + f_in[6]; - temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2.0*f_in[2]*uxuy + f_in[4]*u.x + 2.0*f_in[5]*u.y + f_in[7]; - temp[8] = f_in[0]*ux2*uy2 + 2.0*f_in[1]*u.x*uy2 + 2.0*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4.0*f_in[5]*uxuy + 2.0*f_in[6]*u.y + 2.0*f_in[7]*u.x + f_in[8]; - //back to density-probability functions - f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; - f_in[1] = 0.5*temp[1] + 0.5*temp[3] - 0.5*temp[7] - 0.5*temp[8]; - f_in[2] = 0.5*temp[2] + 0.5*temp[4] - 0.5*temp[6] - 0.5*temp[8]; - f_in[3] = -0.5*temp[1] + 0.5*temp[3] + 0.5*temp[7] - 0.5*temp[8]; - f_in[4] = -0.5*temp[2] + 0.5*temp[4] + 0.5*temp[6] - 0.5*temp[8]; - f_in[5] = 0.25*temp[5] + 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8]; - f_in[6] = -0.25*temp[5] + 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8]; - f_in[7] = 0.25*temp[5] - 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8]; - f_in[8] = -0.25*temp[5] - 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8]; - } +CudaDeviceFunction void relax_and_collide_CM_hydro(real_t f_in[9], real_t tau, vector_t u, vector_t Fhydro, real_t rho) +{ + real_t uxuy = u.x*u.y; + real_t ux2 = u.x*u.x; + real_t uy2 = u.y*u.y; + + real_t s_v = 1./tau; + real_t s_b = omega_bulk; + real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; + real_t temp[9]; real_t cm_eq[9]; + for (int i = 0; i < 9; i++) { + temp[i] = f_in[i];} + //raw moments from density-probability functions + //[m00, m10, m01, m20, m02, m11, m21, m12, m22] + f_in[0] = m00; + f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; + f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; + f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; + f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; + f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; + //central moments from raw moments + temp[3] = f_in[0]*ux2 - 2.0*f_in[1]*u.x + f_in[3]; + temp[4] = f_in[0]*uy2 - 2.0*f_in[2]*u.y + f_in[4]; + temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; + //collision in central moments space + //calculate equilibrium distributions in cm space + cm_eq[0] = m00; + cm_eq[1] = -u.x*(m00 - 1.0); + cm_eq[2] = -u.y*(m00 - 1.0); + cm_eq[3] = m00*ux2 + 1./3.*m00 - ux2; + cm_eq[4] = m00*uy2 + 1./3.*m00 - uy2; + cm_eq[5] = uxuy*(m00 - 1.0); + cm_eq[6] = -u.y*(m00*ux2 + 1./3.*m00 - ux2 - 1./3.); + cm_eq[7] = -u.x*(m00*uy2 + 1./3.*m00 - uy2 - 1./3.); + cm_eq[8] = m00*ux2*uy2 + 1./3.*m00*ux2 + 1./3.*m00*uy2 + 1./9.*m00 - ux2*uy2 - 1./3.*ux2 - 1./3.*uy2; + + //calculate forces in cm space + //He et al. scheme: get_continuous_Maxwellian_DF(dzeta=dzeta, psi=1, u=(ux, uy)) - velocity term + //collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*force_in_cm_space // SOI + f_in[0] = cm_eq[0]; + f_in[1] = 0.5*Fhydro.x/rho + cm_eq[1]; + f_in[2] = 0.5*Fhydro.y/rho + cm_eq[2]; + f_in[3] = 0.5*cm_eq[3]*(s_b + s_v) + 0.5*cm_eq[4]*(s_b - s_v) - temp[3]*(0.5*s_b + 0.5*s_v - 1.0) - 0.5*temp[4]*(s_b - s_v); + f_in[4] = 0.5*cm_eq[3]*(s_b - s_v) + 0.5*cm_eq[4]*(s_b + s_v) - 0.5*temp[3]*(s_b - s_v) - temp[4]*(0.5*s_b + 0.5*s_v - 1.0); + f_in[5] = cm_eq[5]*s_v - temp[5]*(s_v - 1.0); + f_in[6] = 1./6.*Fhydro.y/rho + cm_eq[6]; + f_in[7] = 1./6.*Fhydro.x/rho + cm_eq[7]; + f_in[8] = cm_eq[8]; + + //back to raw moments + temp[0] = f_in[0]; + temp[1] = f_in[0]*u.x + f_in[1]; + temp[2] = f_in[0]*u.y + f_in[2]; + temp[3] = f_in[0]*ux2 + 2.0*f_in[1]*u.x + f_in[3]; + temp[4] = f_in[0]*uy2 + 2.0*f_in[2]*u.y + f_in[4]; + temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; + temp[6] = f_in[0]*ux2*u.y + 2.0*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2.0*f_in[5]*u.x + f_in[6]; + temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2.0*f_in[2]*uxuy + f_in[4]*u.x + 2.0*f_in[5]*u.y + f_in[7]; + temp[8] = f_in[0]*ux2*uy2 + 2.0*f_in[1]*u.x*uy2 + 2.0*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4.0*f_in[5]*uxuy + 2.0*f_in[6]*u.y + 2.0*f_in[7]*u.x + f_in[8]; + //back to density-probability functions + f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; + f_in[1] = 0.5*temp[1] + 0.5*temp[3] - 0.5*temp[7] - 0.5*temp[8]; + f_in[2] = 0.5*temp[2] + 0.5*temp[4] - 0.5*temp[6] - 0.5*temp[8]; + f_in[3] = -0.5*temp[1] + 0.5*temp[3] + 0.5*temp[7] - 0.5*temp[8]; + f_in[4] = -0.5*temp[2] + 0.5*temp[4] + 0.5*temp[6] - 0.5*temp[8]; + f_in[5] = 0.25*temp[5] + 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8]; + f_in[6] = -0.25*temp[5] + 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8]; + f_in[7] = 0.25*temp[5] - 0.25*temp[6] - 0.25*temp[7] + 0.25*temp[8]; + f_in[8] = -0.25*temp[5] - 0.25*temp[6] + 0.25*temp[7] + 0.25*temp[8]; +} - CudaDeviceFunction void relax_and_collide_CM_phase_field(real_t f_in[9], real_t tau, vector_t u, vector_t F_phi, real_t rho) - { - real_t uxuy = u.x*u.y; - real_t ux2 = u.x*u.x; - real_t uy2 = u.y*u.y; - - real_t s_v = 1./tau; - real_t s_b = omega_bulk; - - real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; - - real_t temp[9]; - - for (int i = 0; i < 9; i++) { - temp[i] = f_in[i];} - - //raw moments from density-probability functions - //[m00, m10, m01, m20, m02, m11, m21, m12, m22] - f_in[0] = m00; - f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; - f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; - - // f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; - // f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; - // f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; - // f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; - // f_in[6] = temp[5] + temp[6] - temp[7] - temp[8]; - // f_in[7] = temp[5] - temp[6] - temp[7] + temp[8]; - // f_in[8] = temp[5] + temp[6] + temp[7] + temp[8]; - - // Central moments from raw moments - // temp[0] = m00; - temp[1] = -f_in[0]*u.x + f_in[1]; - temp[2] = -f_in[0]*u.y + f_in[2]; - - // temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3]; - // temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4]; - // temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; - // temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6]; - // temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7]; - // temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8]; - - // Collision in central moments space - // Default is second order integration (trapezoidal discretization) - // collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*F_phi_cm - // relax 1st moments - f_in[0] = m00; - f_in[1] = -F_phi.x*(0.5*s_v - 1.0) - temp[1]*(s_v - 1.0); - f_in[2] = -F_phi.y*(0.5*s_v - 1.0) - temp[2]*(s_v - 1.0); - f_in[3] = 1./3.*m00; - f_in[4] = 1./3.*m00; - f_in[5] = 0; - f_in[6] = 1./6.*F_phi.y; - f_in[7] = 1./6.*F_phi.x; - f_in[8] = 1./9.*m00; - - //back to raw moments - temp[0] = f_in[0]; - temp[1] = f_in[0]*u.x + f_in[1]; - temp[2] = f_in[0]*u.y + f_in[2]; - temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3]; - temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4]; - temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; - temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6]; - temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7]; - temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8]; - - //back to density-probability functions - f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; - f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2; - f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2; - f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2; - f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2; - f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4; - f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4; - f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4; - f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4; - } +CudaDeviceFunction void relax_and_collide_CM_phase_field(real_t f_in[9], real_t tau, vector_t u, vector_t F_phi, real_t rho) +{ + real_t uxuy = u.x*u.y; + real_t ux2 = u.x*u.x; + real_t uy2 = u.y*u.y; + + real_t s_v = 1./tau; + real_t s_b = omega_bulk; + + real_t m00 = f_in[0] + f_in[1] + f_in[2] + f_in[3] + f_in[4] + f_in[5] + f_in[6] + f_in[7] + f_in[8]; + + real_t temp[9]; + + for (int i = 0; i < 9; i++) { + temp[i] = f_in[i];} + + //raw moments from density-probability functions + //[m00, m10, m01, m20, m02, m11, m21, m12, m22] + f_in[0] = m00; + f_in[1] = temp[1] - temp[3] + temp[5] - temp[6] - temp[7] + temp[8]; + f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; + + // f_in[2] = temp[2] - temp[4] + temp[5] + temp[6] - temp[7] - temp[8]; + // f_in[3] = temp[1] + temp[3] + temp[5] + temp[6] + temp[7] + temp[8]; + // f_in[4] = temp[2] + temp[4] + temp[5] + temp[6] + temp[7] + temp[8]; + // f_in[5] = temp[5] - temp[6] + temp[7] - temp[8]; + // f_in[6] = temp[5] + temp[6] - temp[7] - temp[8]; + // f_in[7] = temp[5] - temp[6] - temp[7] + temp[8]; + // f_in[8] = temp[5] + temp[6] + temp[7] + temp[8]; + + // Central moments from raw moments + // temp[0] = m00; + temp[1] = -f_in[0]*u.x + f_in[1]; + temp[2] = -f_in[0]*u.y + f_in[2]; + + // temp[3] = f_in[0]*ux2 - 2*f_in[1]*u.x + f_in[3]; + // temp[4] = f_in[0]*uy2 - 2*f_in[2]*u.y + f_in[4]; + // temp[5] = f_in[0]*uxuy - f_in[1]*u.y - f_in[2]*u.x + f_in[5]; + // temp[6] = -f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 - f_in[3]*u.y - 2*f_in[5]*u.x + f_in[6]; + // temp[7] = -f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy - f_in[4]*u.x - 2*f_in[5]*u.y + f_in[7]; + // temp[8] = f_in[0]*ux2*uy2 - 2*f_in[1]*u.x*uy2 - 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy - 2*f_in[6]*u.y - 2*f_in[7]*u.x + f_in[8]; + + // Collision in central moments space + // Default is second order integration (trapezoidal discretization) + // collide eq: (eye(9)-S)*cm + S*cm_eq + (eye(9)-S/2.)*F_phi_cm + // relax 1st moments + f_in[0] = m00; + f_in[1] = -F_phi.x*(0.5*s_v - 1.0) - temp[1]*(s_v - 1.0); + f_in[2] = -F_phi.y*(0.5*s_v - 1.0) - temp[2]*(s_v - 1.0); + f_in[3] = 1./3.*m00; + f_in[4] = 1./3.*m00; + f_in[5] = 0; + f_in[6] = 1./6.*F_phi.y; + f_in[7] = 1./6.*F_phi.x; + f_in[8] = 1./9.*m00; + + //back to raw moments + temp[0] = f_in[0]; + temp[1] = f_in[0]*u.x + f_in[1]; + temp[2] = f_in[0]*u.y + f_in[2]; + temp[3] = f_in[0]*ux2 + 2*f_in[1]*u.x + f_in[3]; + temp[4] = f_in[0]*uy2 + 2*f_in[2]*u.y + f_in[4]; + temp[5] = f_in[0]*uxuy + f_in[1]*u.y + f_in[2]*u.x + f_in[5]; + temp[6] = f_in[0]*ux2*u.y + 2*f_in[1]*uxuy + f_in[2]*ux2 + f_in[3]*u.y + 2*f_in[5]*u.x + f_in[6]; + temp[7] = f_in[0]*u.x*uy2 + f_in[1]*uy2 + 2*f_in[2]*uxuy + f_in[4]*u.x + 2*f_in[5]*u.y + f_in[7]; + temp[8] = f_in[0]*ux2*uy2 + 2*f_in[1]*u.x*uy2 + 2*f_in[2]*ux2*u.y + f_in[3]*uy2 + f_in[4]*ux2 + 4*f_in[5]*uxuy + 2*f_in[6]*u.y + 2*f_in[7]*u.x + f_in[8]; + + //back to density-probability functions + f_in[0] = temp[0] - temp[3] - temp[4] + temp[8]; + f_in[1] = temp[1]/2 + temp[3]/2 - temp[7]/2 - temp[8]/2; + f_in[2] = temp[2]/2 + temp[4]/2 - temp[6]/2 - temp[8]/2; + f_in[3] = -temp[1]/2 + temp[3]/2 + temp[7]/2 - temp[8]/2; + f_in[4] = -temp[2]/2 + temp[4]/2 + temp[6]/2 - temp[8]/2; + f_in[5] = temp[5]/4 + temp[6]/4 + temp[7]/4 + temp[8]/4; + f_in[6] = -temp[5]/4 + temp[6]/4 - temp[7]/4 + temp[8]/4; + f_in[7] = temp[5]/4 - temp[6]/4 - temp[7]/4 + temp[8]/4; + f_in[8] = -temp[5]/4 - temp[6]/4 + temp[7]/4 + temp[8]/4; +} - CudaDeviceFunction vector_t calcTotalHydrodynamicForceCM(vector_t gradPhi, real_t myPhaseF, real_t rho, real_t tau, real_t mu, real_t p) - { - // arguments: - // gradPhi - Phase field gradients - // myPhaseF - PhaseField value at cell of interest (0,0 or BC like inlet_PhaseF, moving_wall_PhaseF) - - // Fluid Properties: // - // real_t rho - fluid density - // real_t tau - relaxation parameter - // real_t mu - chemical potential - // real_t p - normalized pressure - - real_t Gamma[9], geq[9], mag; // equilibrium, pressure equilibrium, velocity magnitude - real_t F_pressure[2], F_body[2], F_mu[2], F_surf_tension[2]; // Forces - vector_t F_total_hydro; - - real_t stress[3]; // Stress tensor calculation - real_t R[9]; // Populations for MRT relaxation - - // Force Calculations - // eq 19 - real_t density_coeff = (Density_h-Density_l)/(PhaseField_h-PhaseField_l); - F_pressure[0] = (-1.0/3.0) * p * density_coeff * gradPhi.x; - F_pressure[1] = (-1.0/3.0) * p * density_coeff * gradPhi.y; - - F_body[0] = -1.0*(rho-Density_l)*BuoyancyX + rho*GravitationX; - F_body[1] = -1.0*(rho-Density_l)*BuoyancyY + rho*GravitationY; - - // eq 4 - F_surf_tension[0] = mu*gradPhi.x; - F_surf_tension[1] = mu*gradPhi.y; - - // Calculate viscous force - for (int j=0;j PhaseField_h) { - tau = tau_h + 0.5; - } else { - // Inverse update: - //tau = (C - PhaseField_l)/(PhaseField_h - PhaseField_l) * (1.0/tau_h - 1.0/tau_l) + 1.0/tau_l; - //tau = 1.0/tau + 0.5; - // Linear update: - tau = 0.5 + tau_l + C*(tau_h - tau_l); - // Viscosity update: - //DynVisc = Density_l*Viscosity_l + C * (Density_h*Viscosity_h - Density_l*Viscosity_l); - //tau = 3.0 * DynVisc / rho + 0.5; - } + #ifdef OPTIONS_debug + AddToF_pressureX(F_pressure[0]); + AddToF_pressureY(F_pressure[1]); + AddToF_bodyX(F_body[0]); + AddToF_bodyY(F_body[1]); + AddToF_surf_tensionX(F_surf_tension[0]); + AddToF_surf_tensionY(F_surf_tension[1]); + AddToF_muX(F_mu[0]); + AddToF_muY(F_mu[1]); + AddToF_total_hydroX(F_total_hydro.x); + AddToF_total_hydroY(F_total_hydro.y); + #endif - // GRADIENTS AND NORMALS - gradPhi = calcGradPhi(); - nx = gradPhi.x/gradPhi.z; - ny = gradPhi.y/gradPhi.z; - - // CALCULATE FORCES: - real_t density_coeff = (Density_h-Density_l)/(PhaseField_h-PhaseField_l); - F_pressure[0] = (-1.0/3.0) * p * density_coeff * gradPhi.x; - F_pressure[1] = (-1.0/3.0) * p * density_coeff * gradPhi.y; - F_body[0] = -1.0*(rho-Density_l)*BuoyancyX + rho*GravitationX; - F_body[1] = -1.0*(rho-Density_l)*BuoyancyY + rho*GravitationY; - - // Finding viscous force: - for (j=0;j Date: Wed, 20 Dec 2023 13:09:48 +1000 Subject: [PATCH 2/4] + add nan flag for cart grid as well for consistent behaviour --- models/multiphase/d2q9_pf_velocity/Dynamics.R | 2 +- models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt | 16 ++++++---------- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.R b/models/multiphase/d2q9_pf_velocity/Dynamics.R index 789c36aa9..602ac8191 100644 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.R +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.R @@ -88,7 +88,7 @@ if (Options$RT) { # iteration AddStage("BaseIter" , "calcHydroIter" , save=Fields$group %in% c("g","h","Vel","nw") , load=DensityAll$group %in% c("PF","g","h","Vel","nw")) # TODO: is nw needed here? AddStage("PhaseIter" , "calcPhaseFIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("g","h","Vel","nw")) - AddStage("WallIter" , "calcWallPhaseIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("nw","PF")) + AddStage("WallIter" , "calcWallPhaseIter" , save=Fields$group %in% c("PF") , load=DensityAll$group %in% c("nw","PF")) # Purposefully read/write of PF for boundary. complex geom may force RACE condition. } AddAction("Iteration", c("BaseIter", "PhaseIter","WallIter")) diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt index 3484c8e2d..3c91b8c54 100755 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt @@ -221,8 +221,7 @@ CudaDeviceFunction real_t calcMu(real_t myPhase){ } #endif - // eq 5 from "Improved locality of the phase-field lattice Boltzmann - // model for immiscible fluids at high density ratios" + // eq 5 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios" real_t pfavg = 0.5*(PhaseField_l+PhaseField_h); real_t mu = 4.0*(12.0*sigma/W)*(myPhase-PhaseField_l)*(myPhase-PhaseField_h)*(myPhase-pfavg) - (1.5 *sigma*W) * lpPhi; @@ -230,8 +229,7 @@ CudaDeviceFunction real_t calcMu(real_t myPhase){ } CudaDeviceFunction real_t calcGamma(int i, real_t u, real_t v, real_t u2mag){ - // eq 10 from "Improved locality of the phase-field lattice Boltzmann - // model for immiscible fluids at high density ratios" + // eq 10 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios" real_t gamma, tmp; tmp = (d2q9_ex[i]*u+d2q9_ey[i]*v); @@ -240,8 +238,7 @@ CudaDeviceFunction real_t calcGamma(int i, real_t u, real_t v, real_t u2mag){ } CudaDeviceFunction real_t calcF_phi(int i, real_t numerator, real_t nx, real_t ny){ - // eq 7 from "Improved locality of the phase-field lattice Boltzmann - // model for immiscible fluids at high density ratios" + // eq 7 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios" // nx: normalized PhaseField gradient in x direction // this method calculates F_phi in [0-9] lattice coord @@ -252,8 +249,7 @@ CudaDeviceFunction real_t calcF_phi(int i, real_t numerator, real_t nx, real_t n } CudaDeviceFunction vector_t calcF_phi_xy(vector_t gradPhi, real_t myPhaseF, real_t pfavg){ - // eq 7 from "Improved locality of the phase-field lattice Boltzmann - // model for immiscible fluids at high density ratios" + // eq 7 from "Improved locality of the phase-field lattice Boltzmann model for immiscible fluids at high density ratios" // this method calculates F_phi in x,y coordinates real_t nx = gradPhi.x/gradPhi.z; // GradPhi normalized in x, y direction @@ -540,8 +536,8 @@ CudaDeviceFunction void calcWallPhaseIter(){ if ((NodeType & NODE_BOUNDARY) == NODE_Wall) { if (nw_x == 0 && nw_y == 0) { - PhaseF = 0.0; - printf("EDGE CASE: Solid found with no normal at (%.1f %.1f)\n", X, Y); + PhaseF = NAN; + //printf("EDGE CASE: Solid found with no normal at (%.1f %.1f)\n", X, Y); return; } real_t myA, myPhase; From a1be45a016ce8f04cbaa12fb53680b9a125adbd4 Mon Sep 17 00:00:00 2001 From: TravisMitchell Date: Wed, 20 Dec 2023 13:15:22 +1000 Subject: [PATCH 3/4] + comment print statements --- models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt index 3c91b8c54..0b7ab7929 100755 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt @@ -543,7 +543,7 @@ CudaDeviceFunction void calcWallPhaseIter(){ real_t myA, myPhase; myPhase = PhaseF_dyn(nw_x, nw_y); if (isnan(myPhase)) { - printf("EDGE CASE: Solid normal pointing into solid node at (%.1f %.1f)\n", X, Y); + //printf("EDGE CASE: Solid normal pointing into solid node at (%.1f %.1f)\n", X, Y); PhaseF=0.0; return; } @@ -554,7 +554,6 @@ CudaDeviceFunction void calcWallPhaseIter(){ PhaseF = (myA - sqrt( myA*myA - 4.0 * (myA-1.0)*myPhase))/(myA-1.0) - myPhase; } } - if (isnan(PhaseF)) printf("NaN at (%.1f %.1f)\n", X, Y); } #ifdef OPTIONS_CM From f05a11b8f9e14bf15303f8a1bef6438c6282846d Mon Sep 17 00:00:00 2001 From: TravisMitchell Date: Wed, 20 Dec 2023 16:28:09 +1000 Subject: [PATCH 4/4] + adjust calc grad phi for efficiency, maintainability, and consistency --- .../multiphase/d2q9_pf_velocity/Dynamics.c.Rt | 87 +++++-------------- 1 file changed, 24 insertions(+), 63 deletions(-) diff --git a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt index 0b7ab7929..d53213371 100755 --- a/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt +++ b/models/multiphase/d2q9_pf_velocity/Dynamics.c.Rt @@ -40,6 +40,7 @@ #include #define PI 3.1415926535897 +#define cs2 0.33333333