try to fix fluxBC for greyscaleSC model but didn't work; save the work

models/GreyscaleSCModel.cpp

@@ -15,7 +15,7 @@ void DeleteArray( const TYPE *p )
 
 ScaLBL_GreyscaleSCModel::ScaLBL_GreyscaleSCModel(int RANK, int NP, MPI_Comm COMM):
 rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauA_eff(0),tauB_eff(0),Gsc(0),
-rhoA(0),rhoB(0),rhoA_minor(0),rhoB_minor(0),Fx(0),Fy(0),Fz(0),flux(0),dinA(0),doutA(0),dinB(0),doutB(0),GreyPorosity(0),
+rhoA(0),rhoB(0),rhoA_minor(0),rhoB_minor(0),Fx(0),Fy(0),Fz(0),fluxA(0),fluxB(0),dinA(0),doutA(0),dinB(0),doutB(0),GreyPorosity(0),
 Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
 	SignDist.resize(Nx,Ny,Nz);           
@@ -50,7 +50,7 @@ void ScaLBL_GreyscaleSCModel::ReadParams(string filename){
     dinB=rhoB_minor;//inlet density for fluid B
     doutA=rhoA_minor;//outlet denisty for fluid A
     doutB=rhoB;//outlet density for fluid B
-	flux=0.0;
+	fluxA=fluxB=0.0;
 	
 	// ---------------------- Greyscale Model parameters -----------------------//
 	if (greyscaleSC_db->keyExists( "timestepMax" )){
@@ -83,8 +83,11 @@ void ScaLBL_GreyscaleSCModel::ReadParams(string filename){
 	if (greyscaleSC_db->keyExists( "Restart" )){
 		Restart = greyscaleSC_db->getScalar<bool>( "Restart" );
 	}
-	if (greyscaleSC_db->keyExists( "flux" )){
+	if (greyscaleSC_db->keyExists( "fluxA" )){
-		flux = greyscaleSC_db->getScalar<double>( "flux" );
+		fluxA = greyscaleSC_db->getScalar<double>( "fluxA" );
+	}
+	if (greyscaleSC_db->keyExists( "fluxB" )){
+		fluxB = greyscaleSC_db->getScalar<double>( "fluxB" );
 	}
 	if (greyscaleSC_db->keyExists( "tolerance" )){
 		tolerance = greyscaleSC_db->getScalar<double>( "tolerance" );
@@ -188,6 +191,9 @@ void ScaLBL_GreyscaleSCModel::ReadInput(){
         case 3:
             if (rank==0) printf("BoundaryCondition=%i: Constant pressure boundary condition\n",BoundaryCondition);
             break;
+        case 4:
+            if (rank==0) printf("BoundaryCondition=%i: Constant flux boundary condition\n",BoundaryCondition);
+            break;
         default:
             if (rank==0) printf("BoundaryCondition=%i: is currently not supported! Periodic boundary condition is used.\n",BoundaryCondition);
             BoundaryCondition=0;
@@ -912,6 +918,11 @@ void ScaLBL_GreyscaleSCModel::Run(){
 			ScaLBL_Comm->GreyscaleSC_Pressure_BC_z(NeighborList, fqA, fqB,  dinA,  dinB, timestep);
 			ScaLBL_Comm->GreyscaleSC_Pressure_BC_Z(NeighborList, fqA, fqB, doutA, doutB, timestep);
 		}
+		if (BoundaryCondition == 4){
+			dinA = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fqA, fluxA, timestep);
+			dinB = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fqB, fluxB, timestep);
+			ScaLBL_Comm->GreyscaleSC_Pressure_BC_Z(NeighborList, fqA, fqB, doutA, doutB, timestep);
+		}
 		ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(NeighborList, dvcMap, fqA, fqB, DenA, DenB, 0, ScaLBL_Comm->LastExterior(), Np);
 
 		//if (BoundaryCondition > 0){
@@ -925,7 +936,7 @@ void ScaLBL_GreyscaleSCModel::Run(){
 		ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(DenA);
 		ScaLBL_DeviceBarrier();
-        if (BoundaryCondition >0 ){
+        if (BoundaryCondition ==3 || BoundaryCondition ==4){//not necessarily applied to velBC (BC=2)
             if (Dm->kproc()==0){
                 ScaLBL_SetSlice_z(DenA,dinA,Nx,Ny,Nz,0);
             }
@@ -938,7 +949,7 @@ void ScaLBL_GreyscaleSCModel::Run(){
 		ScaLBL_Comm_Regular->SendHalo(DenB);
 		ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(DenB);
-        if (BoundaryCondition >0 ){
+        if (BoundaryCondition ==3 || BoundaryCondition ==4){//not necessarily applied to velBC (BC=2)
             if (Dm->kproc()==0){
                 ScaLBL_SetSlice_z(DenB,dinB,Nx,Ny,Nz,0);
             }
@@ -973,6 +984,11 @@ void ScaLBL_GreyscaleSCModel::Run(){
 			ScaLBL_Comm->GreyscaleSC_Pressure_BC_z(NeighborList, fqA, fqB,  dinA,  dinB, timestep);
 			ScaLBL_Comm->GreyscaleSC_Pressure_BC_Z(NeighborList, fqA, fqB, doutA, doutB, timestep);
 		}
+		if (BoundaryCondition == 4){
+			dinA = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fqA, fluxA, timestep);
+			dinB = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fqB, fluxB, timestep);
+			ScaLBL_Comm->GreyscaleSC_Pressure_BC_Z(NeighborList, fqA, fqB, doutA, doutB, timestep);
+		}
 		ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(dvcMap, fqA, fqB, DenA, DenB, 0, ScaLBL_Comm->LastExterior(), Np);
 
 		//if (BoundaryCondition > 0){
@@ -986,7 +1002,7 @@ void ScaLBL_GreyscaleSCModel::Run(){
 		ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(DenA);
 		ScaLBL_DeviceBarrier();
-        if (BoundaryCondition >0 ){
+        if (BoundaryCondition ==3 || BoundaryCondition ==4){//not necessarily applied to velBC (BC=2)
             if (Dm->kproc()==0){
                 ScaLBL_SetSlice_z(DenA,dinA,Nx,Ny,Nz,0);
             }
@@ -1000,7 +1016,7 @@ void ScaLBL_GreyscaleSCModel::Run(){
 		ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(DenB);
 		ScaLBL_DeviceBarrier();
-        if (BoundaryCondition >0 ){
+        if (BoundaryCondition ==3 || BoundaryCondition ==4){//not necessarily applied to velBC (BC=2)
             if (Dm->kproc()==0){
                 ScaLBL_SetSlice_z(DenB,dinB,Nx,Ny,Nz,0);
             }
@@ -1149,6 +1165,14 @@ void ScaLBL_GreyscaleSCModel::Run(){
 //                fclose(log_file);
 //            }
 //		}
+        if (timestep==2&&BoundaryCondition==4){
+            if (rank==0) printf("    Timestep    dinA    dinB    doutA    doutB\n");
+        } 
+		if (timestep%analysis_interval==0){
+            if (BoundaryCondition==4){
+                if (rank==0) printf("    %i    %.3g    %.3g    %.3g    %.3g\n",timestep,dinA,dinB,doutA,doutB); 
+            }
+        }
 
 		if (timestep%visualization_interval==0){
             WriteOutput();

models/GreyscaleSCModel.h

@@ -42,7 +42,8 @@ public:
     double rhoA,rhoB;
     double rhoA_minor,rhoB_minor;//dissolved density
 	double tolerance;
-	double Fx,Fy,Fz,flux;
+	double Fx,Fy,Fz;
+    double fluxA,fluxB;
 	double dinA,doutA;
 	double dinB,doutB;
     double GreyPorosity;