diff --git a/solvers_post/LagrangianCMCFoam/Peta.H b/solvers_post/LagrangianCMCFoam/Peta.H index 891f671..3f4fbc3 100644 --- a/solvers_post/LagrangianCMCFoam/Peta.H +++ b/solvers_post/LagrangianCMCFoam/Peta.H @@ -4,141 +4,136 @@ // eta = (1/etamax)*etacount { - scalar betaftn=0.0, gammaCMC=0.0, alphaCMC=0.0, betaCMC=0.0, PTOTAL=0.0; - scalar loop_counter = 0, mf_diff=0.0, mfnew=0.0; - forAll(mf, celli) //cell number loop + forAll(mf, celli) //cell number loop + { + scalar mfi = mf[celli]; //karam, test + + Switch pdfCalculated = false; + + scalarField pdf(etamax+1,0.0); + + label loop_counter = 0; + do { - mf_temp[celli] = mf[celli]; //karam, test + //Adjust mixture fraction and its variance + if(mfVar[celli] < 1.0e-30) + { + mfVar[celli] = 1.0e-30; + } -Re_calculate_pdf:; + if(mfi < 0.0) + { + mfi = 0.0; + } - //Adjust mixture fraction and its variance - if(mfVar[celli] < 1.0e-30) - { - mfVar[celli] = 1.0e-30; - } + //calculate alpha and beta for PDF function + scalar gammaCMC = mfi*(1.0-mfi) / mfVar[celli] - 1.0; + scalar betaCMC = (1.0-mfi) * gammaCMC; + scalar alphaCMC = mfi * gammaCMC; - if(mf_temp[celli] < 0.0) - { - mf_temp[celli] = 0.0; - } - - if(mf_temp[celli] < 5.0e-4) //karam, 5e-4 -> 1e-4 ) - { - pdf = 0.0; - pdf[0] = 1.0 / deta[0];//1; - deltaftn[celli] = 0; - goto Peta_NextCell; - } - - gammaCMC = mf_temp[celli]*(1.0-mf_temp[celli]) / mfVar[celli] - 1.0; - - if(gammaCMC <= 0.0) - { + if(mfi < 5.0e-4) //karam, 5e-4 -> 1e-4 ) + { + pdf = 0.0; + pdf[0] = 1.0 / deta[0];//1; + deltaftn[celli] = 0; + pdfCalculated = true; + } + else if(gammaCMC <= 0.0) + { // Info<<"GammaCMC is lower than zero -- two delta function"; - pdf = 0.0; - pdf[0] = 0.5 * 1.0/deta[0]; //0.5; - pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5; - deltaftn[celli] = 1; - goto Peta_NextCell; - } - - //calculate alpha and beta for PDF function - betaCMC = (1.0-mf_temp[celli]) * gammaCMC; - alphaCMC = mf_temp[celli] * gammaCMC; - - if(alphaCMC <= 1.0 && betaCMC <= 1.0) - { + pdf = 0.0; + pdf[0] = 0.5 * 1.0/deta[0]; //0.5; + pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5; + deltaftn[celli] = 1; + pdfCalculated = true; + } + else if(alphaCMC <= 1.0 && betaCMC <= 1.0) + { // Info<<"AlphaCMC and betaCMC are lower than one -- two delta function"; - pdf = 0.0; - pdf[0] = 0.5 * 1.0/deta[0]; //0.5; - pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5; - deltaftn[celli] = 1; - goto Peta_NextCell; - } + pdf = 0.0; + pdf[0] = 0.5 * 1.0/deta[0]; //0.5; + pdf[etamax] = 0.5 * 1.0/deta[etamax]; //0.5; + deltaftn[celli] = 1; + pdfCalculated = true; + } + else if(alphaCMC < 1.0 && betaCMC >= 1.0) + { + pdf = 0.0; + pdf[0] = 1.0/deta[0]; //1; + deltaftn[celli] = 0; + pdfCalculated = true; + } + else if(alphaCMC >= 1.0 && betaCMC < 1.0) + { + pdf = 0.0; + pdf[etamax] = 1.0/deta[etamax]; //1; + deltaftn[celli] = 0; + pdfCalculated = true; + } + else + { + scalar betaftn = std::exp( gammaln(alphaCMC) + gammaln(betaCMC) - gammaln(alphaCMC + betaCMC) ); //Numerial recipe 2nd Edition 206page - else if(alphaCMC < 1.0 && betaCMC >= 1.0) + if (betaftn > 1.0e-50) { - pdf = 0.0; - pdf[0] = 1.0/deta[0]; //1; - deltaftn[celli] = 0; - goto Peta_NextCell; - } + pdf[0] = 0.0; + pdf[etamax] = 0.0; + for(label etacount = 1 ; etacount < etamax ; etacount++) + { + pdf[etacount] = Foam::pow( etaValue[etacount] , alphaCMC - 1.0 ) + * Foam::pow ( 1.0 - etaValue[etacount] , betaCMC - 1.0 ) / betaftn ; + } - else if(alphaCMC >= 1.0 && betaCMC < 1.0) + deltaftn[celli] = 2; + + scalarField f(etamax+1, 1.0); + scalar pTotal = integration(deltaftn[celli], MFcut, Neta, pdf, f); + + for(label etacount = 0; etacount <= etamax ; etacount++) + { + pdf[etacount] = pdf[etacount]/pTotal; + } + + scalar mfnew = integration(deltaftn[celli], MFcut, Neta, pdf, etaValue); + scalar mf_diff = (mf[celli] - mfnew)/mf[celli]; + + if(std::fabs(mf_diff) < 0.00001) // 1e-5) + { + pdfCalculated = true; + } + else + { + mfi = mfi + 0.5*mf_diff*mf[celli]; + } + } + else { - pdf = 0.0; - pdf[etamax] = 1.0/deta[etamax]; //1; - deltaftn[celli] = 0; - goto Peta_NextCell; + loop_counter = 1000; } + } - betaftn = std::exp( gammaln(alphaCMC) + gammaln(betaCMC) - gammaln(alphaCMC + betaCMC) ); //Numerial recipe 2nd Edition 206page + loop_counter++; - if(betaftn < 1.0e-50) + } while ((!pdfCalculated) && (loop_counter < 100)); + + if(!pdfCalculated) + { + pdf = 0.0; + for(label etacount = 0; etacount<=etamax ; etacount++) //set one delta function + { + if( etaValue[etacount] >= mfi) { -set_onedelta:; - pdf = 0.0; - for(label etacount = 0; etacount<=etamax ; etacount++) //set one delta function - { - if( etaValue[etacount] >= mf_temp[celli]) - { - pdf[etacount] = 1.0/deta[etacount]; //1; - deltaftn[celli] = 0; - goto Peta_NextCell; - } - } - + pdf[etacount] = 1.0/deta[etacount]; //1; + deltaftn[celli] = 0; + break; } - - for(label etacount=0 ; etacount <= etamax ; etacount++) - { - if(etacount == 0 || etacount == etamax) - { - pdf[etacount] = 0.0; - } - else - { - pdf[etacount] = Foam::pow( etaValue[etacount] , alphaCMC - 1.0 ) * Foam::pow ( 1.0 - etaValue[etacount] , betaCMC -1.0 ) / betaftn ; - } - } - deltaftn[celli] = 2; - - PTOTAL = 0.0; //PTOTAL initialization in each cell - - f = 1.0; - - PTOTAL = integration(deltaftn[celli], MFcut, Neta, pdf, f); - - for(label etacount=0; etacount<=etamax ; etacount++) - { - pdf[etacount] = pdf[etacount]/PTOTAL; - } - - mfnew = integration(deltaftn[celli], MFcut, Neta, pdf, etaValue); - - mf_diff = (mf[celli] - mfnew)/mf[celli]; - - if(std::fabs(mf_diff) > 0.00001)//1e-5) - { - mf_temp[celli] = mf_temp[celli] + 0.5*mf_diff*mf[celli]; - if(loop_counter > 100) - { - Info<<"Force set pdf to one delta"<