[bigphil]

Try the following:

const polyPatchList& patches = mesh.boundaryMesh();
forAll(patches,patchI)
{
    if (patches[patchI].name() == "midup")
    {
        fixedGradientFvPatchScalarField& gradcMinusPatch =
            refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]);
        scalarField& gradcMinusField = gradcMinusPatch.gradient();

        forAll(gradcMinusField, faceI)
        {                                           
            if (cMinus.boundaryField()[patchI][faceI] < 0)
            {
                // Set gradient for this face to zero
                gradcMinusField[faceI] = 0;
            }
            else
            {
                // Set gradient for this face to some calculated value; I assume eKbT, Dinv, and currentj are all scalars (not scalarFields)
                gradcMinusField[faceI] = eKbT*cMinus.boundaryField()[patchI][faceI]*Ue.boundaryField()[patchI].snGrad()[faceI] + Dinv*currentj;
            }
        }        
    }    
}

Also, I suggest reading the OpenFOAM coding style: following it will make your code more readable for others

Philip

[babakflame]

Many Thanks Philip

I am working on your proposed snippet and appologize for the caused confusions.

I had read that page, but maybe I need to re-read again

Regards

[babakflame]

Philip,

thanks for your reply. I think the second forAll loop

forAll(patchI, faceI){
                 if(cMinus.boundaryField()[patchI][faceI] < 0){     
                    gradcMinusField [faceI] = 0;
                  }
                  else{
                  gradcMinusField [faceI] = eKbT * cMinus.boundaryField()[patchI][faceI] * Ue.boundaryField()    [patchI][faceI].snGrad() + Dinv * currentj; //3.13e13;   //midup
               }
         }

brings this error:

Boundary.H: In function ‘int main(int, char**)’:
/home/babak/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/UList.H:394:36: error: request for member ‘size’ in ‘patchI’, which is of non-class type ‘Foam::label {aka int}’
     for (Foam::label i=0; i<(list).size(); i++)
                                    ^
Boundary.H:18:4: note: in expansion of macro ‘forAll’
    forAll(patchI, faceI){
    ^
In file included from interFoamEHDIon5New.C:97:0:
Boundary.H:16:17: warning: unused variable ‘gradcMinusField’ [-Wunused-variable]
    scalarField& gradcMinusField = gradcMinusPatch.gradient();
                 ^
In file included from interFoamEHDIon5New.C:77:0:
/home/babak/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/lnInclude/readTimeControls.H:38:8: warning: unused variable ‘maxDeltaT’ [-Wunused-variable]
 scalar maxDeltaT =
        ^
make: *** [Make/linux64GccDPOpt/interFoamEHDIon5New.o] Error 1

Regards

[bigphil]

Sorry, "forAll(patchI, faceI)" should have been "forAll(gradcMinusField, faceI)": I have editted this in my previous post.

Also, for your second error, please be careful, you have written:

else
{
    gradcMinusField [faceI] =
        eKbT*cMinus.boundaryField()[patchI][faceI]*Ue.boundaryField()[patchI][faceI].snGrad() + Dinv * currentj;
}

whereas I have proposed:

else
{
    gradcMinusField[faceI] =
        eKbT*cMinus.boundaryField()[patchI][faceI]*Ue.boundaryField()[patchI].snGrad()[faceI] + Dinv*currentj;
}

Notice the differences.

[bigphil]

Actually, as "snGrad()" constructs and returns the full patch field, it will be more efficient to store the snGrad as follows:

const polyPatchList& patches = mesh.boundaryMesh();
forAll(patches,patchI)
{
    if (patches[patchI].name() == "midup")
    {
        fixedGradientFvPatchScalarField& gradcMinusPatch =
            refCast<fixedGradientFvPatchScalarField>(cMinus.boundaryField()[patchI]);
        scalarField& gradcMinusField = gradcMinusPatch.gradient();
        const scalarField UeSnGrad = Ue.boundaryField()[patchI].snGrad();

        forAll(gradcMinusField, faceI)
        {                                           
            if (cMinus.boundaryField()[patchI][faceI] < 0)
            {
                // Set gradient for this face to zero
                gradcMinusField[faceI] = 0;
            }
            else
            {
                // Set gradient for this face to some calculated value; I assume eKbT, Dinv, and currentj are all scalars (not scalarFields)
                gradcMinusField[faceI] = eKbT*cMinus.boundaryField()[patchI][faceI]*UeSnGrad[faceI] + Dinv*currentj;
            }
        }        
    }    
}

This will avoid the full Ue patch snGrad field being created for every time a face value is set.

[babakflame]

Thank you philip for the comments.

I understood the points you raised in the snippet. The code compiles successfully.



I am running my cases in parallel. Do I need to add reduce function after the last line of each section (after gradcMinusField=.. and gradcPlusField=.. lines) to make sure that the master processor gets the correct value??

sth like this:

reduce(gradcMinusField = maxOp<scalar>());

Regards

[bigphil]

Without seeing the rest of your code, I cannot say whether you need to explicitly perform a parallel sync; however, from looking at your code snippet, the only terms that might need it are those you have defined before the code snippet e.g. eKbT, Dinv, currentj.

How do you calculate these terms?

[babakflame]

Philip,

These terms are all constant scalars (eKbT, Dinv) or constant outward vector (currentj) with known values given in the case folder (defined this way in createFields.H file)

Regards