Snippet for redefining fixedGradient boundary condition for a patch inside solver

babakflame · March 3, 2017, 18:36

Dear Fellows

I want to redifine fixedGradient boundary condition for a specified patch inside my source code. i.e. whenever the code sees "fixedGradient" in 0/cPlus file (cPlus is a volScalarField and a dependent variable)

My snippet now looks like:

Code:

const polyPatchList& patches = mesh.boundaryMesh();
   forAll (patches,patchI){
     
      if(cPlus.boundaryField()[patchI].type() == 'fixedGradient'){


           fixedGradientPatchField& fgPatchField
         (
          refCast<fixedGradientPatchField>(cPlus.boundaryField()[patchI])
         );
          fgPatchField.gradient() = fvc::snGrad(cPlus) + cPlus * fvc::snGrad(Ue);
       }
    }

In above snippet, I want to tell the code that whenever it saw "fixedGradient" for a patch inside 0/cPlus file calculate the following relation at the boundary surface and assign it as boundary condition:

$\frac{d c^+}{d n} + c^+ \frac{d \varphi}{d n}$ in which

& $\varphi$ are dependent scalar variables that should be calculated at the boundary patch from previous time step.

However the code compilation returns this error:

Code:

createFields.H:501:3: error: ‘fgPatchField’ was not declared in this scope
In file included from interFoamEHDIon.C:79:0:

My first question is How should I declare fgPatchField ? since I tried volScalarField or volVectorField but did not work. Moreover it does not accept fvc::snGrad. Should I try fvc::reconstruct(snGrad(cPlus))?

The second error is:

Code:

createFields.H:532:12: error: ‘fixedGradientFvPatchVectorField’ was not declared in this scope

createFields.H:533:74: error: no matching function for call to ‘refCast(Foam::fvPatchField<Foam::Vector<double> >&)’

By the way, My OpenFoam version is 2.1.x

babakflame · March 4, 2017, 18:09

Dear Fellows

I made my snippet working. This snippet applies above formulation to a specific patch that has been predefined in 0/cPlus file with type fixedGradient. I haven't tested it though.

You need to add this header to your main.C file as well:

Code:

#include "fixedGradientFvPatchFields.H"

The snippet now looks like:

Code:

const polyPatchList& patches = mesh.boundaryMesh();
       forAll (patches,patchI){      

    if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){
        fixedGradientFvPatchScalarField& gradcPlusPatch =
            refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]);
        scalarField& gradcPlusField = gradcPlusPatch.gradient();
        gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();
     }
       }

I will keep ya updated with the results.

Keep Foaming Fellows.

Zeppo · March 4, 2017, 18:27

Code:

fixedGradientPatchField& fgPatchField
(
    refCast<fixedGradientPatchField>(cPlus.boundaryField()[patchI])
);
fgPatchField.gradient() = fvc::snGrad(cPlus) + cPlus * fvc::snGrad(Ue);

->

Code:

fixedGradientFvPatchField<scalar>& fgPatchField
(
    refCast<fixedGradientFvPatchField<scalar> >(cPlus.boundaryField()[patchI])
);
fgPatchField.gradient() = 
    cPlus.boundaryField()[patchI].fvPatchField<scalar>::snGrad()
  + cPlus.boundaryField[patchI].snGrad() * Ue.boundaryField[patchI].snGrad();

Should be compiled but I am not sure the code corresponds to your formula. I think you'll have to develop a custom boundary condition.

babakflame · March 4, 2017, 18:41

Thanks Sergei

Actually, I have doubts on my final line in the snippet addressing formulation i.e.:

Code:

gradcPlusField = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();

I am gonna try the formulation of last line and give the feedbacks.

Regards

bigphil · March 6, 2017, 03:34

Quote:

Originally Posted by babakflame

Thanks Sergei

Actually, I have doubts on my final line in the snippet addressing formulation i.e.:

Code:

gradcPlusField = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();

I am gonna try the formulation of last line and give the feedbacks.

Regards

Hi babakflame,

What equation are you trying to implement as the boundary condition?

I don't think the following makes sense:

Code:

cPlus.boundaryField()[patchI].snGrad() = 
       cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad();

which is essentially what you are doing.

This would only converge if "cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad()" is zero.

Maybe you are trying to set something like the following?

Code:

cPlus.boundaryField()[patchI].snGrad() = 
       - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();

Philip

babakflame · March 6, 2017, 17:08

Hi BigPhil

Many thanks for your reply.

In my domain, I have upper and lower walls which fluid is passing inbetween.

Actually I am trying to set a fixed flux (at lower wall) and zero flux (at upper wall) boundary conditions for scalar cPlus and zero flux (at lower wall) and fixed flux (at upper wall) boundary conditions for scalar cMinus.

My final boundary condition relations are as follows: (

and

are constants.) Actually,

is the value of flux. However, these relations are the boundary conditions I have for cPlus and cMinus.
$\left\{ \begin{array}{ll} \frac{d c^+}{dn} = - A_1 c^+ \frac{d \varphi}{dn}\\ \frac{d c^-}{dn} = - B_1 + A_1 c^- \frac{d \varphi}{dn} \end{array} \right.$ for the upper wall

and the following relations:
$\frac{d c^+}{dn} = - B_1 - A_1 c^+ \frac{d \varphi}{dn}$

$\frac{d c^-}{dn} = - A_1 c^- \frac{d \varphi}{dn}$ for the lower wall.
n is the direction normal to the wall and outwards.

At the moment, I am trying this snippet in my main.C file:
( $\varphi$ is denoted by Ue in my code which is another scalar (dependent variable)).

Code:

const label& ID= mesh.boundaryMesh().findPatchID("midup");
if(cPlus.boundaryField()[ID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& gradcPlusPatch =
		refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[ID]);
	scalarField& gradcPlusField = gradcPlusPatch.gradient();
	gradcPlusField = - eKbT * gradcPlusPatch * Ue.boundaryField()[ID].snGrad();
}
else if(cMinus.boundaryField()[ID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& gradcMinusPatch =
		refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[ID]);
	scalarField& gradcMinusField = gradcMinusPatch.gradient();
	gradcMinusField =  Dinv * 2.8e-11 + eKbT * gradcMinusPatch * Ue.boundaryField()[ID].snGrad();
}
       
const label& patchID = mesh.boundaryMesh().findPatchID("middown");
if(cPlus.boundaryField()[patchID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& cPlusPatch =
		refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchID]);
	scalarField& gradcPlusField = cPlusPatch.gradient();
	gradcPlusField =  Dinv * 2.8e-11 + eKbT * cPlusPatch * Ue.boundaryField()[patchID].snGrad();
	 }

else if(cMinus.boundaryField()[patchID].type()=="fixedGradient"){
	fixedGradientFvPatchScalarField& cMinusPatch =
		refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchID]);
	scalarField& gradcMinusField = cMinusPatch.gradient();
		gradcMinusField = - eKbT * cMinusPatch * Ue.boundaryField()[patchID].snGrad();
}

Attached, you can see a picture of my domain with coordinates and n direction.

Another quick question:

I have defined the type of boundary condition as fixedGradient in 0/cPlus and 0/cMinus for upper and lower walls. However, the value of this gradient in 0 folder (I need to assign a value there like 0) still affects my results. When we use this method, the code should use my snippet as a custom boundary condition not the initial value in 0 folder, Isn't it?

Regards

Zeppo · March 7, 2017, 11:48

Quote:

Originally Posted by babakflame

$\left\{ \begin{array}{ll} \frac{d c^+}{dn} = - A_1 c^+ \frac{d \varphi}{dn}\\ \frac{d c^-}{dn} = - B_1 + A_1 c^- \frac{d \varphi}{dn} \end{array} \right.$

Now as we know what Bobi really wanted to get Philip's suggestion looks rather feasible

Quote:

Originally Posted by bigphil

Maybe you are trying to set something like the following?

Code:

cPlus.boundaryField()[patchI].snGrad() = 
       - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();

bigphil · March 7, 2017, 15:47

Hi babakflame,

A couple of comments/questions:

in the picture of your domain, the normal vectors "n" point inwards (i.e. in to the domain); however, OpenFOAM convention assumes the normals point outwards;
where do you place this code snippet? If you place it in createFields.H then the boundary gradient will only be set once at the start of the simulation e.g. dcPlusdn will be set based on the initial value of cPlus on the boundary, so if cPlus is zero at the start then dcPlusdn will also be set to zero. So maybe you want to put this code at the start of the time-loop so it is called once at the start of every time-step: this would set dcPlusdn based on the previous time-step value of cPlus (an explicit method), or if you need it more implicit then you would need to call this code multiple times per time-step (e.g. put a outer loop around your cPlus equation).

Philip

babakflame · March 7, 2017, 23:17

Dear Philip

Many thanks for your reply.

As a respond to your comments:

I wasn't aware of the direction of n in OpenFoam in the way that you had mentioned, but by resimulating my case couple of times, I figured out that I need to multiply the lower wall boundary conditions by a negative. (This with your description now makes sense to me).
I put this snippet inside my time-loop before the pimple algorithm inner-loop. So, I think it updates the boundary conditions once at each time step.

I have a doubt here:
Although I have a gradient of $\varphi$ inside my domain, I think the term $\frac{d \varphi}{d n}$ within the boundary condition relations is not calculated (The code considers it as zero). Basically, its just the value of flux

in my above boundary conditions that affects the gradient of

or

.

Just for clarification, this is my final snippet:

Code:

const polyPatchList& patches = mesh.boundaryMesh(); 
    forAll (patches,patchI){   
 
    if(patches[patchI].name()=="midup"){     
 
         if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcPlusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
             scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
             gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();  // midup
         } 
 
         if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
              fixedGradientFvPatchScalarField& gradcMinusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
              scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
              gradcMinusField =  cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + Dinv * 2.8e-11;  //midup 
          } 
      }      
      else if(patches[patchI].name()=="middown"){ 
          if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcPlusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
             scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
             gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + Dinv * 2.8e-11;  // middown 
 
          } 
 
          if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
             fixedGradientFvPatchScalarField& gradcMinusPatch = 
                 refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
              scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
              gradcMinusField =  cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() ;  //middown 
 
            } 
        } 
   }

A question: Is the way I have implemented $c^+ \frac{d \varphi}{d n}$ or $c^+ \frac{d \varphi}{d n}$ in above snippet correct or should I just use :

Code:

fvc::snGrad(Ue)

for $\frac{d \varphi}{d n}$ term?

Regards

bigphil · March 8, 2017, 03:35

The snGrad at a patch is correctly calculated from:

Code:

Ue.boundaryField()[patchI].snGrad()

If you were to use:

Code:

fvc::snGrad(Ue)

this would give you the snGrad on every internal face and boundary face; you could then take the patch values; of course, in your case you only want the patch values so the first makes sense (note: they will both give the same values).

So, I am not sure if your method is working for you now but the general implementation seems OK; if the values are still not as expected, I would suggest printing them to the screen and figuring out whats wrong.

Philip

babakflame · March 8, 2017, 10:21

Dear Philip

Many thanks. I am gonna investigate term by term (in the way that you mentioned) to see exactly whats going on.

Regards

babakflame · May 19, 2017, 14:22

Dear Fellows

As abovementioned, I have successfully implemented fixed flux boundary condition. However, after non-dimensionalizing my solver and boundary conditions, some scales are introduced and boundary conditions (for example for lower wall) becomes as follow

$\frac {d c^-}{d n}= - \frac {e \varphi_0}{k_B T} c^- \frac{d \varphi}{d n}$
Here, the coefficient $\frac {e \varphi_0}{k_B T}$ has the value 0f

which causes my solver crashing.

If I ran the non-dimensional code with the non-dimensional boundary condition without this coefficient every thing is fine. I will post the error soon. However, is there any way to improve the solver performance?

What happens is that my courant number of non-dimensional code starts oscilating before divergence in initial time steps.

I am aware of some approaches such as:

grid refinement

increasing relTol to 0.1

to improve the performance of the solver, However, none of these worked for me.

I tried to increase that coefficient gradually, however, did not work. here is the error:

Code:

Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#4  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#5  Foam::fvMatrix<double>::solve() at ??:?
#6  
 at ??:?
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"

Any suggestion or previous experience?

babakflame · May 19, 2017, 21:52

Even underrelaxing the related equations as follows did not work.

Code:

relaxationFactors
  {
     fields
      {

      }
     equations
     {
          Ue               0.1;
         "(cPlus|cMinus)*" 0.1;
    }
  }

bigphil · May 22, 2017, 09:27

Quote:

Originally Posted by babakflame

Even underrelaxing the related equations as follows did not work.

Code:

relaxationFactors
  {
     fields
      {

      }
     equations
     {
          Ue               0.1;
         "(cPlus|cMinus)*" 0.1;
    }
  }

Hi babakflame,

Can you post the code you are using?
This might help give an idea of how the stability of your method could be improved (also it will help others who might want to try something similar).

Philip

babakflame · May 25, 2017, 14:19

Dear Philip

Many thanks for your help.

Here is the main part of code (.C) file:

Code:

int main(int argc, char *argv[]) 
{ 
    #include "setRootCase.H" 
    #include "createTime.H" 
    #include "createMesh.H" 
 
     pimpleControl pimple(mesh); 
 
    #include "initContinuityErrs.H" 
    #include "createFields.H" 
    #include "readTimeControls.H" 
    #include "correctPhi.H" 
    #include "CourantNo.H" 
    #include "setInitialDeltaT.H" 
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
 
    Info<< "\nStarting time loop\n" << endl; 
 
    while (runTime.run()) 
    { 
        #include "readTimeControls.H" 
        #include "CourantNo.H" 
        #include "alphaCourantNo.H" 
        #include "setDeltaT.H" 
 
        runTime++; 
 
        Info<< "Time = " << runTime.timeName() << nl << endl; 
 
        twoPhaseProperties.correct(); 
 
    const polyPatchList& patches = mesh.boundaryMesh(); 
           forAll (patches,patchI){   
 
        if(patches[patchI].name()=="midup"){     
 
            if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcPlusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
                scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
                gradcPlusField = - cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad();  // midup 
 
             } 
 
                  if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcMinusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
                scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
                gradcMinusField = cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;  // 1D case 
 
             } 
              }      
        else if(patches[patchI].name()=="middown"){ 
            if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcPlusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); 
                scalarField& gradcPlusField = gradcPlusPatch.gradient(); 
                gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + DcinfLinv * 3.68e-14;  //channelBB case 
                                 
             } 
 
                  if(cMinus.boundaryField()[patchI].type()=="fixedGradient"){ 
                fixedGradientFvPatchScalarField& gradcMinusPatch = 
                    refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]); 
                scalarField& gradcMinusField = gradcMinusPatch.gradient(); 
                gradcMinusField = - cMinus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() ;  //middown  channelBB case 
                               
             } 
        } 
    } 
 
 
        #include "alphaEqnSubCycle.H" 
        #include "ElectricEqn.H" 
         //Info<< "ElectricEqn over " << nl << endl; 
        // --- Pressure-velocity PIMPLE corrector loop 
        while (pimple.loop()) 
        { 
         
            #include "UEqn.H" 
 
            // --- Pressure corrector loop 
            while (pimple.correct()) 
            { 
                #include "pEqn.H" 
            } 
 
            if (pimple.turbCorr()) 
            { 
                turbulence->correct(); 
            } 
        } 
 
        runTime.write(); 
 
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" 
            << "  ClockTime = " << runTime.elapsedClockTime() << " s" 
            << nl << endl; 
    } 
 
    Info<< "End\n" << endl; 
 
    return 0; 
}

And the ElectricEqn.H is as follows:

Code:

        surfaceScalarField cPlusFlux  = -DLU * fvc::interpolate(cPlus)*mesh.magSf()*fvc::snGrad(Ue);

        surfaceScalarField cMinusFlux =  DLU * fvc::interpolate(cMinus)*mesh.magSf()*fvc::snGrad(Ue);

 //           surfaceScalarField cPlusFlux  = -fvc::interpolate(sgm)*mesh.magSf()*fvc::snGrad(Ue);
            fvScalarMatrix cPlusEqn
            (
                fvm::ddt(cPlus)
          + fvm::div(phi, cPlus)
              + fvc::div(cPlusFlux)
           - fvm::laplacian(kappa, cPlus)    
            );
//            cPlusEqn.relax();
            cPlusEqn.solve();


            fvScalarMatrix cMinusEqn
            (
                fvm::ddt(cMinus)
          + fvm::div(phi, cMinus)
              + fvc::div(cMinusFlux)
           - fvm::laplacian(kappa, cMinus)    
            );

//        cMinusEqn.relax();
            cMinusEqn.solve();


        rhoE = (cPlus - cMinus);
            
              // solve
            fvScalarMatrix UeEqn
            (
                fvm::laplacian(eps,Ue) == -rhoE * gamma
            );
//            UeEqn.relax();
            UeEqn.solve();

As you can see above, there are three variables : cPlus, cMinus and Ue that are solved once in each pimple loop. There are two transport equations for cPlus and cMinus and one Poisson equation for Ue.
gamma, DLU, Kappa and eta are the non-dimensional groups (constants) appear as non-dimensionalization.
The UEqn.H is as follow:

Code:

 surfaceScalarField muEff
    (
        "muEff",
        twoPhaseProperties.muf()
      + fvc::interpolate(rho*turbulence->nut())
    );




    tmp1 = rhoE * E;

    tmp2 = rhoE * Eext;

       Ux = U.component(vector::X);

       rhoEUx = rhoE * Ux;

 
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muEff, U)
      - (fvc::grad(U) & fvc::grad(muEff))
      - tmp2 * eta
      - tmp1 * eta
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)

                ) * mesh.magSf()
            )

        );
    }

The settings I am using for fvSolutions and fvSchemes are as follows:

Code:

solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       0.001;
            relTol          0;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       0.0001;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.05;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration false;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-08;
            relTol          0;
            nVcycles        2;
            smoother        GaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-08;
        relTol          0;
        maxIter         20;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-09;
        relTol          0;
        nSweeps         1;
    }

    "(Ue|UeExt)"
    {
//        solver         PCG;
//        preconditioner   DIC;
//        tolerance        1e-14;
//        relTol           0.2;
        solver smoothSolver;
        smoother GaussSeidel;
        preconditioner DIC;
        tolerance 1e-14;
        relTol 0.1;
    };

//    rhoE
    "(cPlus|cMinus)"
    {
        solver         PBiCG;
        preconditioner   DILU;
        tolerance        1e-16;
        relTol           0.000;
    };


    "(k|B|nuTilda)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    nAlphaCorr      1;
    nAlphaSubCycles 3;
    cAlpha          1;

    pRefPoint      (0 5e-5 1e-5);
    pRefValue      0;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
/*         phi        0.1;
         phiFinal     1;
         cPlus      0.1;
         cPlusFinal   1;
         cMinus     0.1;
         cMinusFinal  1;
*/     }
}

Code:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default          Gauss linear;
    div(rho*phi,U)   Gauss upwind;
//    div(phi,rhoE)    Gauss upwind;
    div(phi,cMinus)    Gauss upwind;
    div(phi,cPlus)    Gauss upwind;
    div(phi,alpha)   Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
    rhoE;
}

As you can see in the main file (.C) one, for some patches I have defined fixedGradient which will be overwritten by my defined boundary condition.

for the Ue boundary condition on that patches, I want to use Dirichlet boundary condition, with high values, However, the code diverges, just for low values it will continue the simulation.

The code is based on interFoam in OpenFoam version 2.1.x.

Many thanks

babakflame · May 26, 2017, 12:07

As you can see above philip, the Ue variable is found through the Poisson equation (laplace with source term).

My boundary condition for Ue is Dirichlet type. I have tried the followings:

smooth solver (Gauss seidel with relTol = 0.1)
refining grid close to boundaries
relaxing the equation
higher order fv schemes

None of these worked. The initial residual of either cPlus or cMinus or Ue starts to grow and reaching to one which throws a floating point exception error.

bigphil · May 26, 2017, 13:03

Quote:

Originally Posted by babakflame

As you can see above philip, the phi variable is found through the poisson equation (laplacae with source term).

My boundary condition for phi is Dirichlet type. I have tried the followings:

smooth solver (Gauss seidel with relTol = 0.1)
refining grid close to boundaries
relaxing the equation
higher order fv schemes

None of these worked. The initial residual of either cPlus or cMinus or Ue starts to grow and reaching to one which throws a floating point exception error.

Hi babakflame,

Can you also post the output from the solver? What exactly is diverging? the linear solver for one of the equations? or the outer pimple loop? or time-steps?

From a quick look through the code you have posted, it seems that you are solving three coupled PDEs in a staggered manner; of course, convergence in these cases may not be easy to achieve using a staggered method (if the PDEs are strongly coupled).
A fruitful direction might be performing outer iterations over the equations with significant under-relaxation of the key variables and/or under-relaxation of the boundary conditions (I would start with the BCs). At least try figure out what exactly is causing the divergence.

Philip

babakflame · May 26, 2017, 14:00

Dear Philip

many thanks for your help.

First of all, these are boundary files in 0 folder:
cPlus

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      cPlus;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 -3 1 0 0 1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

cMinus:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      cMinus;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 -3 1 0 0 1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedGradient;
        gradient        uniform 0;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      Ue;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 2 -3 0 0 -1 0];

internalField   uniform 0;

boundaryField
{
    leftup
    {
        type            zeroGradient;
    }
    midup
    {
        type            fixedValue;
        value           uniform 20;
    }
    rightup
    {
        type            zeroGradient;
    }
    leftdown
    {
        type            zeroGradient;
    }
    middown
    {
        type            fixedValue;
        value           uniform -20;
    }
    rightdown
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

This is a channel flow case with these additional equations. As you can see in 0/Ue file only I am capable to run the case with values of Ue in order of 10 for midup and middown patches (middle part of channel). Even with this order of values in the middle of simulation, the divergence happens with the following error:

Code:

Courant Number mean: 5.03946621477e+51 max: 7.02462243817e+55
Interface Courant Number mean: 0 max: 0
deltaT = 3.52225240641e-108
Time = 2.16860254237690730505505598558

--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 30 to 31 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 31 to 32 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 32 to 33 to distinguish between timeNames at time 2.16860254238
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
--> FOAM Warning : 
    From function Time::operator++()
    in file db/Time/Time.C at line 1024
    Increased the timePrecision from 33 to 34 to distinguish between timeNames at time 2.16860254238
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 2.62480833735e-11, Final residual = 1.57380692088e-25, No Iterations 1
DILUPBiCG:  Solving for cMinus, Initial residual = 1, Final residual = 1.16439623825e-18, No Iterations 3
smoothSolver:  Solving for Ue, Initial residual = 1, Final residual = 0.0999529656162, No Iterations 33
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  void Foam::dot<Foam::Vector<double>, Foam::Vector<double> >(Foam::Field<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type>&, Foam::UList<Foam::Vector<double> > const&, Foam::UList<Foam::Vector<double> > const&) at ??:?
#4  void Foam::dot<Foam::Vector<double>, Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#5  Foam::tmp<Foam::GeometricField<Foam::innerProduct<Foam::Vector<double>, Foam::Vector<double> >::type, Foam::fvPatchField, Foam::volMesh> > Foam::operator&<Foam::Vector<double>, Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6  
 at ??:?
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8  
 at ??:?
Floating point exception (core dumped)

This error is similar to the one that arises in higher orders of Ue in the very first time steps.

Regards

babakflame · May 26, 2017, 14:08

This is the situation of residuals a few time steps before divergence:

Code:

Courant Number mean: 0.0135133107361 max: 0.0199948706422
Interface Courant Number mean: 0 max: 0
deltaT = 9.99750062485e-05
Time = 2.16858285429

MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 2.31041374006e-06, Final residual = 3.23788533007e-17, No Iterations 3
DILUPBiCG:  Solving for cMinus, Initial residual = 0.961364873856, Final residual = 3.50366344753e-18, No Iterations 5
smoothSolver:  Solving for Ue, Initial residual = 0.946454822981, Final residual = 0.0928075889037, No Iterations 14
GAMG:  Solving for p_rgh, Initial residual = 0.99287837977, Final residual = 0.0320233117181, No Iterations 3
time step continuity errors : sum local = 1.7866809096e-06, global = 4.76494332971e-10, cumulative = 4.7773692495e-10
GAMG:  Solving for p_rgh, Initial residual = 0.347493694957, Final residual = 0.0169072858213, No Iterations 2
time step continuity errors : sum local = 1.99215077851e-06, global = 1.3176994971e-09, cumulative = 1.79543642205e-09
GAMGPCG:  Solving for p_rgh, Initial residual = 0.251746573971, Final residual = 6.03393734406e-09, No Iterations 18
time step continuity errors : sum local = 7.01209128599e-13, global = 2.68332995554e-15, cumulative = 1.79543910538e-09
ExecutionTime = 1580.05 s  ClockTime = 1585 s

Courant Number mean: 0.0135540948435 max: 0.507795570766
Interface Courant Number mean: 0 max: 0
deltaT = 1.96873675724e-05
Time = 2.16860254165

MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
MULES: Solving for alpha1
Liquid phase volume fraction = 0  Min(alpha1) = 0  Max(alpha1) = 0
Two conductive liquids
DILUPBiCG:  Solving for cPlus, Initial residual = 4.54849797817e-07, Final residual = 4.1721661111e-17, No Iterations 2
DILUPBiCG:  Solving for cMinus, Initial residual = 0.991818666534, Final residual = 3.3247608014e-20, No Iterations 5
smoothSolver:  Solving for Ue, Initial residual = 0.993844252023, Final residual = 0.0980642096096, No Iterations 10
GAMG:  Solving for p_rgh, Initial residual = 0.999978321915, Final residual = 0.0355731771938, No Iterations 3
time step continuity errors : sum local = 0.0332311172846, global = 6.40661804589e-06, cumulative = 6.40841348499e-06
GAMG:  Solving for p_rgh, Initial residual = 0.251482888814, Final residual = 0.0110229824325, No Iterations 3
time step continuity errors : sum local = 0.0239228239629, global = 2.60788641251e-06, cumulative = 9.0162998975e-06
GAMGPCG:  Solving for p_rgh, Initial residual = 0.145311847528, Final residual = 7.59825317339e-09, No Iterations 20
time step continuity errors : sum local = 1.75845002939e-08, global = -7.27185444404e-11, cumulative = 9.01622717895e-06
ExecutionTime = 1581.15 s  ClockTime = 1586

As you can see above, the intial residuals of p_rgh, cMinus and Ue start to reach 1 just before crashing.

babakflame · May 26, 2017, 16:35

Quote:

Originally Posted by bigphil

Hi babakflame,

A fruitful direction might be performing outer iterations over the equations with significant under-relaxation of the key variables and/or under-relaxation of the boundary conditions (I would start with the BCs). At least try figure out what exactly is causing the divergence.

Philip

I have applied under-relaxation for the 3 equations of cPlus, cMinus and Ue. But how can I under-relax boundary conditions??

If you meant Fields, I have applied that too with values around 0.1, no results.

Is there any syntax besides fields and equations (in fvSolutions file) that does under-relaxation in openFoam?

March 3, 2017, 18:36	Snippet for redefining fixedGradient boundary condition for a patch inside solver	#1
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Fellows I want to redifine fixedGradient boundary condition for a specified patch inside my source code. i.e. whenever the code sees "fixedGradient" in 0/cPlus file (cPlus is a volScalarField and a dependent variable) My snippet now looks like: Code: const polyPatchList& patches = mesh.boundaryMesh(); forAll (patches,patchI){ if(cPlus.boundaryField()[patchI].type() == 'fixedGradient'){ fixedGradientPatchField& fgPatchField ( refCast<fixedGradientPatchField>(cPlus.boundaryField()[patchI]) ); fgPatchField.gradient() = fvc::snGrad(cPlus) + cPlus * fvc::snGrad(Ue); } } In above snippet, I want to tell the code that whenever it saw "fixedGradient" for a patch inside 0/cPlus file calculate the following relation at the boundary surface and assign it as boundary condition: $\frac{d c^+}{d n} + c^+ \frac{d \varphi}{d n}$ in which & $\varphi$ are dependent scalar variables that should be calculated at the boundary patch from previous time step. However the code compilation returns this error: Code: createFields.H:501:3: error: ‘fgPatchField’ was not declared in this scope In file included from interFoamEHDIon.C:79:0: My first question is How should I declare fgPatchField ? since I tried volScalarField or volVectorField but did not work. Moreover it does not accept fvc::snGrad. Should I try fvc::reconstruct(snGrad(cPlus))? The second error is: Code: createFields.H:532:12: error: ‘fixedGradientFvPatchVectorField’ was not declared in this scope createFields.H:533:74: error: no matching function for call to ‘refCast(Foam::fvPatchField<Foam::Vector<double> >&)’ By the way, My OpenFoam version is 2.1.x Last edited by babakflame; March 4, 2017 at 17:17.

March 8, 2017, 03:35		#10
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	The snGrad at a patch is correctly calculated from: Code: Ue.boundaryField()[patchI].snGrad() If you were to use: Code: fvc::snGrad(Ue) this would give you the snGrad on every internal face and boundary face; you could then take the patch values; of course, in your case you only want the patch values so the first makes sense (note: they will both give the same values). So, I am not sure if your method is working for you now but the general implementation seems OK; if the values are still not as expected, I would suggest printing them to the screen and figuring out whats wrong. Philip babakflame likes this.

May 19, 2017, 14:22		#12
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Fellows As abovementioned, I have successfully implemented fixed flux boundary condition. However, after non-dimensionalizing my solver and boundary conditions, some scales are introduced and boundary conditions (for example for lower wall) becomes as follow $\frac {d c^-}{d n}= - \frac {e \varphi_0}{k_B T} c^- \frac{d \varphi}{d n}$ Here, the coefficient $\frac {e \varphi_0}{k_B T}$ has the value 0f which causes my solver crashing. If I ran the non-dimensional code with the non-dimensional boundary condition without this coefficient every thing is fine. I will post the error soon. However, is there any way to improve the solver performance? What happens is that my courant number of non-dimensional code starts oscilating before divergence in initial time steps. I am aware of some approaches such as: grid refinement increasing relTol to 0.1 to improve the performance of the solver, However, none of these worked for me. I tried to increase that coefficient gradually, however, did not work. here is the error: Code: Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 in "/lib/x86_64-linux-gnu/libc.so.6" #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:? #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:? #5 Foam::fvMatrix<double>::solve() at ??:? #6 at ??:? #7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6" Any suggestion or previous experience? Last edited by babakflame; May 19, 2017 at 21:53.

May 19, 2017, 21:52		#13
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Even underrelaxing the related equations as follows did not work. Code: relaxationFactors { fields { } equations { Ue 0.1; "(cPlus\|cMinus)*" 0.1; } }

May 26, 2017, 12:07		#16
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	As you can see above philip, the Ue variable is found through the Poisson equation (laplace with source term). My boundary condition for Ue is Dirichlet type. I have tried the followings: smooth solver (Gauss seidel with relTol = 0.1) refining grid close to boundaries relaxing the equation higher order fv schemes None of these worked. The initial residual of either cPlus or cMinus or Ue starts to grow and reaching to one which throws a floating point exception error. Last edited by babakflame; May 26, 2017 at 16:14.

March 4, 2017, 18:09		#2
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Fellows I made my snippet working. This snippet applies above formulation to a specific patch that has been predefined in 0/cPlus file with type fixedGradient. I haven't tested it though. You need to add this header to your main.C file as well: Code: #include "fixedGradientFvPatchFields.H" The snippet now looks like: Code: const polyPatchList& patches = mesh.boundaryMesh(); forAll (patches,patchI){ if(cPlus.boundaryField()[patchI].type()=="fixedGradient"){ fixedGradientFvPatchScalarField& gradcPlusPatch = refCast<fixedGradientFvPatchScalarField>(cPlus.boundaryField()[patchI]); scalarField& gradcPlusField = gradcPlusPatch.gradient(); gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad(); } } I will keep ya updated with the results. Keep Foaming Fellows. Luttappy likes this.

March 4, 2017, 18:41		#4
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Thanks Sergei Actually, I have doubts on my final line in the snippet addressing formulation i.e.: Code: gradcPlusField = cPlus.boundaryField()[patchI] * Ue.boundaryField()[patchI].snGrad() + cPlus.boundaryField()[patchI].snGrad(); I am gonna try the formulation of last line and give the feedbacks. Regards Last edited by babakflame; March 4, 2017 at 19:45.

March 7, 2017, 15:47		#8
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Hi babakflame, A couple of comments/questions: in the picture of your domain, the normal vectors "n" point inwards (i.e. in to the domain); however, OpenFOAM convention assumes the normals point outwards; where do you place this code snippet? If you place it in createFields.H then the boundary gradient will only be set once at the start of the simulation e.g. dcPlusdn will be set based on the initial value of cPlus on the boundary, so if cPlus is zero at the start then dcPlusdn will also be set to zero. So maybe you want to put this code at the start of the time-loop so it is called once at the start of every time-step: this would set dcPlusdn based on the previous time-step value of cPlus (an explicit method), or if you need it more implicit then you would need to call this code multiple times per time-step (e.g. put a outer loop around your cPlus equation). Philip cheng1988sjtu, babakflame and piu58 like this.

March 8, 2017, 10:21		#11
babakflame Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15	Dear Philip Many thanks. I am gonna investigate term by term (in the way that you mentioned) to see exactly whats going on. Regards

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Wind turbine simulation	Saturn	CFX	58	July 3, 2020 01:13
Problem with cyclic boundaries in Openfoam 1.5	fs82	OpenFOAM	36	January 7, 2015 00:31
Error finding variable "THERMX"	sunilpatil	CFX	8	April 26, 2013 07:00
[blockMesh] Cyclic BC's: Possible face ordering problem? (Channel flow)	sega	OpenFOAM Meshing & Mesh Conversion	3	September 28, 2010 12:46
[Gmsh] Import gmsh msh to Foam	adorean	OpenFOAM Meshing & Mesh Conversion	24	April 27, 2005 08:19