Updating Boundary Conditions Each Iteration

Kombinator · July 25, 2018, 03:31

Dear Philip,

Sorry to bother you. It seems that I stuck again but right now with tangential laplacian. So tangential laplacian is:

$\nabla^{2}_{t}$ $\varphi$ = $\nabla$ $\cdot$ ((

- n $\cdot$ n) $\cdot$ (

- n $\cdot$ n) $\nabla$ $\varphi$ )

That means that we want to take a divergence operation, from this line

Code:

/ Get tangential component of gradPotE at the boundary
vectorField tangGradPotE =
    ((I - sqr(nf)) & (I - sqr(nf)) & gradPotE.boundaryFieldRef()[patchID].patchInternalField());

But it is not possible due to we store our variable tangGradPotE in a vectorField and for fvc::div(tangGradPotE) we need tangGradPotE to be a surfaceScalarField. Moreover, I understand that fvc::div(tangGradPotE) is not a correct operation at all in my case, because I need values only on the boundary.
So considering the above it seems that:

I cannot take laplacian before the projection like fvc:: laplacian(PotE) (because firstly I need to project my vector).
I cannot take laplacian afterwards, because I store my variables in a vectorField.

That's why I don't understand what to do in this case. Any help will be appreciated. Thanks in advance.

Best regards,
Tar

Kombinator · July 25, 2018, 09:39

Small update.
I tried to go the second way and take laplacian after the projection. I am not sure that this was the right way, but it was worth to try. Right now it looks something like this:

Code:

// Calculate project tensor
surfaceTensorField pp = (I - sqr(mesh.Sf())) & (I - sqr(mesh.Sf()));
// Calculate lpalacian with projection
volScalarField lp = fvc::laplacian(pp, PotE);

fixedGradientFvPatchScalarField& PotEpatch =
    refCast<fixedGradientFvPatchScalarField>
    (
        PotE.boundaryFieldRef()[patchID]
    );

// Choose laplacian at the boundary   
scalarField lps = lp.boundaryFieldRef()[patchID].patchInternalField(); 

// Set the normal gradient equal to the lpalacian
forAll(PotEpatch, faceI)
{
    PotEpatch.gradient()[faceI] = lps[faceI];
}

It complies but when I run solver it shows the error due to this part:

Code:

sqr(mesh.Sf())

It seems that "sqr" is compatible only with patch.nf(). Maybe anybody could advice something how to overcome this. Thanks in advance.

Best regards,
Tar

bigphil · July 25, 2018, 09:45

Hi Tar,

I am not quite sure what a tangential Laplacian operator is (maybe it is like div(grad(PotE).T() ...?), but your proposed formulae/code look very strange to me.

I suggest you examine how others in literature have performed this discretisation, to be clear on what coefficients you expect etc.

Philip

Kombinator · July 25, 2018, 11:46

Quote:

Originally Posted by bigphil

Hi Tar,

I am not quite sure what a tangential Laplacian operator is (maybe it is like div(grad(PotE).T() ...?), but your proposed formulae/code look very strange to me.

I suggest you examine how others in literature have performed this discretisation, to be clear on what coefficients you expect etc.

Philip

Dear Philip, thank you for your answer.
The tangential laplacian means that we use only the tangential part of the gradient (along the wall) when we calculate the laplacian:

$\nabla^{2}_{t}$ $\varphi$ = $\nabla^{}_{t}$ $\cdot$ ( $\nabla^{}_{t}$ $\varphi$ )

Unfortunately I couldn't find any information about tangential laplacian so let me explain my formula by myself:
A vector has a normal and a tangential component.A tangential component of any vector is a full vector minus the normal one:

$\vec{X}^{}_{t}$ = $\vec{X}$ - $\vec{X}^{}_{n}$ = $\vec{X}$ - ( $\vec{X}$ $\cdot$ $\vec{n}$ ) $\vec{n}$ = (I - $\vec{n}$ $\otimes$ $\vec{n}$ ) $\vec{X}$

That means that (I - $\vec{n}$ $\otimes$ $\vec{n}$ ) is a projection tensor.

Laplacian is a dot product of two vectors:

$\nabla^{2}_{}$ $\varphi$ = $\nabla$ $\cdot$ $\nabla$ $\varphi$

So, we want to use only the tangential part. That means that we have to project our gradients before taking the laplacian.

$\nabla^{2}_{t}$ $\varphi$ = ((I - $\vec{n}$ $\otimes$ $\vec{n}$ ) $\nabla$ ) $\cdot$ ((I - $\vec{n}$ $\otimes$ $\vec{n}$ ) $\nabla$ $\varphi$ ) = $\nabla$ $\cdot$ ((I - $\vec{n}$ $\otimes$ $\vec{n}$ )(I - $\vec{n}$ $\otimes$ $\vec{n}$ ) $\nabla$ $\varphi$ )

I understand how to calculate the whole right part besides a divergence: ((I - $\vec{n}$ $\otimes$ $\vec{n}$ )(I - $\vec{n}$ $\otimes$ $\vec{n}$ ) $\nabla$ $\varphi$ )
I do not understand how to take a divergence of this term, because we store this part in vectorField:

Code:

vectorField tangGradPotE =
    ((I - sqr(nf))&(I – sqr(nf)) & gradPotE.boundaryFieldRef()[patchID].patchInternalField());

And it is not allowed to take a divergence of vectorField. This is the main question. Hope this post is not very messy.

Best regards,
Tar

bigphil · July 26, 2018, 04:33

Hi Tar,

OK, I understand a bit better now.

Quote:

And it is not allowed to take a divergence of vectorField.

Why not? The divergence of a vector field is a scalar field.

You can explicitly calculate the term $\nabla \cdot [(I - nn) \cdot (I - nn) \cdot (\nabla \phi) ]$ as:

Code:

volScalarField divTangGrad =
    fvc::div
    (
        Sf
     &  (
             (I - sqr(n)) & (I - sqr(n))
           & fvc::interpolate(fvc::grad(phi))
        )
    );

where 'phi' is a volScalarField and the area vectors and unit area vectors are set as:

Code:

const surfaceVectorField& Sf = mesh.Sf();
const surfaceVectorField n = mesh.Sf()/mesh.magSf();

Note that you have to interpolate grad(phi) to the cell faces because we can't multiply a volField (defined at cell centres) by a surface field (defined at face centres). Also, in OpenFOAM the differential operators are actually evaluated in integral form (i.e. they multiple my areas and volumes) so when we pass a surfaceField to the fvc::div operator, it expects us to have multiplied by the area vector, and then it in facts just performs a surface integration.
I think this makes sense but please double check it.

Philip

Kombinator · July 26, 2018, 10:22

Quote:

Originally Posted by bigphil

Hi Tar,

I think this makes sense but please double check it.

Philip

Dear Philip,

Thank you very much for your answer. I have checked it and concluded that this part is not right:

Code:

 Sf &  ((I - sqr(n)) & (I - sqr(n))
           & fvc::interpolate(fvc::grad(phi))

because it is equal to zero in any case. The point is that we are multiplying our projected vector on the plane surface by the normal Sf to the same surface.

(I - sqr(n)) & (I - sqr(n)) & fvc::interpolate(fvc::grad(phi)) - projected vector
Sf - normal vector

That means that we have a scalar multiplication between two orthogonal vectors and it equals to zero. I also tried this in openFoam and got zero everywhere.

Can we do some tricks and change area normal with something else?

bigphil · July 26, 2018, 10:33

Ah yes, that makes sense. The projection removes the normal component but then the divergence takes only the normal component.

Hmnn actually, we can note the following equivalence: $[(I - nn) \cdot (I - nn)] = (I - nn)$ .

And then also we can note that the divergence of the tangential component of the gradient of a scalar is always zero: $\nabla \cdot [(I - nn) \cdot \nabla \phi] = 0$ .

So now I realise that I once again don't understand what exactly this tangential Laplacian operator is meant to be :O

I guess if you know what you want to do for a simple structured grid, then generalising it to unstructured grids in OpenFOAM should not be an issue if you can explain the behaviour you want for the simple case.

Philip

Kombinator · July 27, 2018, 05:44

Dear Philip,

Thank you very much for your effort. Anyway, I learned a lot from this topic. I will try to do something and if I succeed I will share my ideas in this topic.

Best regards,
Tar

pvpnrao · September 9, 2018, 13:29

Quote:

Originally Posted by marupio

I'm not infront of my FOAM machine right now, but I seem to recall an issue where the boundary condition would update only once per timestep. I was doing multiple iterations per timestep, and I had to call some function to force it to update. I forget the details... but if you encounter the situation where it only seems to update once, then let me know, I'll dig for the answer.

-Dave

David
Can you let me know the function which can be called to force a boundary update during the iteration witina time loop

shiyu · September 7, 2020, 09:50

Hi Philip,

Hope you are well.
I found that you have provided several answers about updating the boundary conditions inside the OpenFOAM solver. I have also tried the code in your replies, which is working quite well. Thanks for your help.

However, I have a question related to the code (I quoted) below.
Does the value we extract from the boundary cell or face is based on the last timestep?
I am wondering if you know how to update them at the same iteration (implicitly), like the fixedGradient B.C. ?
What I mean is as below,
U_f = a*U_c + b
where the U_f and Uc are the boundary face values and boundary cell values, and most importantly, they are required to updated at the same time/iteration. The a and b are constant coefficients. That is to say, we need to implement this updating implicitly, rather than explicitly(if I am right).

Many thanks,
Shiyu

Quote:

Originally Posted by bigphil

Hi Wei,

You can set non-uniform boundary values based on anything you want:

Code:

// find ID of patch
label patchID = mesh.boundaryMesh().findPatchID("patch_of_interest");
// check patch has been found
if(patchID == -1)
{
  FatalError << "patch not found!" << exit(FatalError);
}
 
// set value on patch
forAll(U.boundaryField()[patchID], facei)
{
   vector faceOldVel = U.boundaryField()[patchID][facei];
   vector faceCellOldVel =
         U.internalField()[mesh.boundaryMesh()[patchID].faceCells()[facei]];
   vector faceSnGrad = U.boundaryField()[patchID].snGrad()[facei];
   
   // set whatever you want here 
   U.boundaryField()[patchID][facei] == faceOldVel + faceCellOldVel + faceSnGrad;
}

Best regards,
Philip

bigphil · September 7, 2020, 10:18

Hi Shiyu,

In the code you have shown, the boundary value is updated explicitly using the latest value of the internal field. So when the equation is solved after this, the boundary value does not change and so this approach would be explicit rather than implicit.

To make it implicit, you could do one of two things:

Add a loop around this code and solving the equation; once this loop converges then this would become implicit in time
Or, re-formulate this condition as a "fixedGradient"-type boundary and set the appropriate surface normal gradient i.e. re-arranging U_f = a*U_c + b into the form (U_f - U_c)/delta = gradientYouWant, where delta is the distance from the cell-centre to the boundary; you can then update this gradient as required at the start of each time-step, for example.

Philip

mostanad · December 13, 2021, 03:55

Quote:

Originally Posted by bigphil

Hi Shiyu,

In the code you have shown, the boundary value is updated explicitly using the latest value of the internal field. So when the equation is solved after this, the boundary value does not change and so this approach would be explicit rather than implicit.

To make it implicit, you could do one of two things:

Add a loop around this code and solving the equation; once this loop converges then this would become implicit in time
Or, re-formulate this condition as a "fixedGradient"-type boundary and set the appropriate surface normal gradient i.e. re-arranging U_f = a*U_c + b into the form (U_f - U_c)/delta = gradientYouWant, where delta is the distance from the cell-centre to the boundary; you can then update this gradient as required at the start of each time-step, for example.

Philip

Hi Philip,
I know it is a while that you are answering to the questions. But for me, this kind of BC update, every iteration for the fixedValue velocity, does not work. So, I hope you can help me to figure out what is going on here:

Code:

        // find ID of patch
        label patchID = mesh.boundaryMesh().findPatchID("sphere");
        // check patch has been found
        if(patchID == -1)
        {
            FatalError << "patch not found!" << exit(FatalError);
        }        
        Info<<"patchID: "<<patchID<<U.boundaryField()[patchID].type()<<endl;

        vectorField newPatchValues(U.boundaryField()[patchID].size(), vector::zero);

        forAll(U.boundaryField()[patchID],i)
        {
            newPatchValues[i] = vector(0, -0.1, 0);
        }

        U.boundaryField()[patchID] == newPatchValues;

As I see, the patch velocity does not change through the iterations and the 0/fixedValue remains unchanged.

Cheers,
Mohammad

bigphil · December 13, 2021, 08:06

What U boundary condition are you using for the "sphere" patch?
Where in your solver do you call this code?
You can check the current boundary condition and its values in the code, by printing it to the standard output, e.g.

Code:

Info<< " U.boundaryField()[patchID] = " <<  U.boundaryField()[patchID] << endl;

mostanad · December 13, 2021, 17:49

Quote:

Originally Posted by bigphil

What U boundary condition are you using for the "sphere" patch?
Where in your solver do you call this code?
You can check the current boundary condition and its values in the code, by printing it to the standard output, e.g.

Code:

Info<< " U.boundaryField()[patchID] = " <<  U.boundaryField()[patchID] << endl;

By this printing output, it gives me the following lines:

Code:

U.boundaryField()[patchID] = type            fixedValue;
value           uniform (0 0 0);

Also, I am using the code before the momentum predictor step, like this in pisoFOAM:

Code:

while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"


        label patchID = mesh.boundaryMesh().findPatchID("sphere");
        if(patchID == -1)
        {
            FatalError << "patch not found!" << exit(FatalError);
        }        
        vectorField newPatchValues(U.boundaryField()[patchID].size(), vector::zero);
        forAll(U.boundaryField()[patchID],i)
        {
            newPatchValues[i] = vector(0, -0.1, 0);
        }
        U.boundaryField()[patchID] == newPatchValues;
        Info<< " U.boundaryField()[patchID] = " <<  U.boundaryField()[patchID] << endl;


        // Pressure-velocity PISO corrector
        {
            #include "UEqn.H"

            // --- PISO loop
            while (piso.correct())
            {
                #include "pEqn.H"
            }
        }

        laminarTransport.correct();
        turbulence->correct();

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

bigphil · December 14, 2021, 06:09

Which version of OpenFOAM are you using?

Maybe the "==" operator has been redefined, and the values need to be set using another approach.

mostanad · December 14, 2021, 06:48

Quote:

Originally Posted by bigphil

Which version of OpenFOAM are you using?

Maybe the "==" operator has been redefined, and the values need to be set using another approach.

I am using OF 5.x. GeometricField.C shows the following definitions. Should I change it?

Code:

template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
(
    const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
)
{
    const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();

    checkField(*this, gf, "==");

    // Only assign field contents not ID

    ref() = gf();
    boundaryFieldRef() == gf.boundaryField();

    tgf.clear();
}


template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
(
    const dimensioned<Type>& dt
)
{
    ref() = dt;
    boundaryFieldRef() == dt.value();
}

bigphil · December 14, 2021, 07:45

Both of those functions refer to GeometricFields, however, in your case you are calling the "==" operator for a fvPatchField. From a quick check of the OF 5 Github code (https://github.com/OpenFOAM/OpenFOAM...fvPatchField.C) this should work; I am not sure what the problem is; however, instead you could try set the values face by face:

Code:

        forAll(U.boundaryField()[patchID], faceI)
        {
            U.boundaryField()[patchID][faceI] = vector(0, -0.1, 0);
        }
        Info<< "U.boundaryField()[patchID] = " << U.boundaryField()[patchID] << endl;

mostanad · December 14, 2021, 08:09

Quote:

Originally Posted by bigphil

Both of those functions refer to GeometricFields, however, in your case you are calling the "==" operator for a fvPatchField. From a quick check of the OF 5 Github code (https://github.com/OpenFOAM/OpenFOAM...fvPatchField.C) this should work; I am not sure what the problem is; however, instead you could try set the values face by face:

Code:

        forAll(U.boundaryField()[patchID], faceI)
        {
            U.boundaryField()[patchID][faceI] = vector(0, -0.1, 0);
        }
        Info<< "U.boundaryField()[patchID] = " << U.boundaryField()[patchID] << endl;

Dear Philip,
I have tried this as well. But nothing changed. Can you please just try it on a titorial case and see what happens?
I really need to have this changing BC option. Thank you for your help.

bigphil · December 14, 2021, 08:28

The following works in OpenFOAM-7 (I don't have OF-5), where I edited icoFoam.C:

Code:

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        const label patchID = mesh.boundaryMesh().findPatchID("movingWall");
        if (patchID == -1)
        {
            FatalError
                << "patch not found!" << exit(FatalError);
        }
        U.boundaryFieldRef()[patchID] == vector(0.001, 0, 0);
        Info<< "U.boundaryField()[patchID] = " << U.boundaryField()[patchID] << endl;

	// Momentum predictor                                                                                                                                                                                

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

...

where notice that you need to use "boundaryFieldRef()" as opposed to "boundaryField()" to get non-const access to the patch (I guess this may have been your problem). In fact, the code will not compile if I try assign "boundaryField()" with "==" as this is read only.

The output on the cavity case for the last time-step is:

Code:

Time = 0.5

Courant Number mean: 0.000222305 max: 0.000852393
U.boundaryField()[patchID] = type            fixedValue;
value           uniform (0.001 0 0);

smoothSolver:  Solving for Ux, Initial residual = 2.41821e-07, Final residual = 2.41821e-07, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 5.28524e-07, Final residual = 5.28524e-07, No Iterations 0
DICPCG:  Solving for p, Initial residual = 9.68699e-07, Final residual = 9.68699e-07, No Iterations 0
time step continuity errors : sum local = 8.74098e-12, global = -3.24571e-22, cumulative = 1.37819e-21
DICPCG:  Solving for p, Initial residual = 1.02537e-06, Final residual = 8.70812e-07, No Iterations 1
time step continuity errors : sum local = 7.96321e-12, global = -7.50117e-22, cumulative = 6.28071e-22
ExecutionTime = 0.08 s  ClockTime = 0 s

End

mostanad · December 14, 2021, 17:23

Quote:

Originally Posted by bigphil

where notice that you need to use "boundaryFieldRef()" as opposed to "boundaryField()" to get non-const access to the patch (I guess this may have been your problem). In fact, the code will not compile if I try assign "boundaryField()" with "==" as this is read only.

Thanks Philip,
Your are the master of OpenFOAM coding. Now it is working for OF5.x with U.boundaryFieldRef. I just insert the code to be in access for everybody working with this version:

Code:

label patchID = mesh.boundaryMesh().findPatchID("sphere");
// check patch has been found
if(patchID == -1)
{
     FatalError << "patch not found!" << exit(FatalError);
}
vectorField newPatchValues(U.boundaryField()[patchID].size(), vector::zero);
forAll(U.boundaryField()[patchID],i)
{
     newPatchValues[i] = vector(0, -0.1, 0);
}
U.boundaryFieldRef()[patchID] == newPatchValues

July 25, 2018, 03:31		#61
Kombinator New Member Artem Join Date: Apr 2014 Posts: 29 Rep Power: 12	Dear Philip, Sorry to bother you. It seems that I stuck again but right now with tangential laplacian. So tangential laplacian is: $\nabla^{2}_{t}$ $\varphi$ = $\nabla$ $\cdot$ (( - n $\cdot$ n) $\cdot$ ( - n $\cdot$ n) $\nabla$ $\varphi$ ) That means that we want to take a divergence operation, from this line Code: / Get tangential component of gradPotE at the boundary vectorField tangGradPotE = ((I - sqr(nf)) & (I - sqr(nf)) & gradPotE.boundaryFieldRef()[patchID].patchInternalField()); But it is not possible due to we store our variable tangGradPotE in a vectorField and for fvc::div(tangGradPotE) we need tangGradPotE to be a surfaceScalarField. Moreover, I understand that fvc::div(tangGradPotE) is not a correct operation at all in my case, because I need values only on the boundary. So considering the above it seems that: I cannot take laplacian before the projection like fvc:: laplacian(PotE) (because firstly I need to project my vector). I cannot take laplacian afterwards, because I store my variables in a vectorField. That's why I don't understand what to do in this case. Any help will be appreciated. Thanks in advance. Best regards, Tar

September 7, 2020, 10:18		#71
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Hi Shiyu, In the code you have shown, the boundary value is updated explicitly using the latest value of the internal field. So when the equation is solved after this, the boundary value does not change and so this approach would be explicit rather than implicit. To make it implicit, you could do one of two things: Add a loop around this code and solving the equation; once this loop converges then this would become implicit in time Or, re-formulate this condition as a "fixedGradient"-type boundary and set the appropriate surface normal gradient i.e. re-arranging U_f = a*U_c + b into the form (U_f - U_c)/delta = gradientYouWant, where delta is the distance from the cell-centre to the boundary; you can then update this gradient as required at the start of each time-step, for example. Philip shiyu likes this.

December 13, 2021, 08:06		#73
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	What U boundary condition are you using for the "sphere" patch? Where in your solver do you call this code? You can check the current boundary condition and its values in the code, by printing it to the standard output, e.g. Code: Info<< " U.boundaryField()[patchID] = " << U.boundaryField()[patchID] << endl;

December 14, 2021, 07:45		#77
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Both of those functions refer to GeometricFields, however, in your case you are calling the "==" operator for a fvPatchField. From a quick check of the OF 5 Github code (https://github.com/OpenFOAM/OpenFOAM...fvPatchField.C) this should work; I am not sure what the problem is; however, instead you could try set the values face by face: Code: forAll(U.boundaryField()[patchID], faceI) { U.boundaryField()[patchID][faceI] = vector(0, -0.1, 0); } Info<< "U.boundaryField()[patchID] = " << U.boundaryField()[patchID] << endl;

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Impinging Jet Boundary Conditions	Anindya	Main CFD Forum	25	February 27, 2016 12:58
OpenFOAM Variable Velocity Boundary Conditions	NickolasPl	OpenFOAM Programming & Development	2	May 19, 2011 05:37
Concentric tube heat exchanger (Air-Water)	Young	CFX	5	October 6, 2008 23:17
[Commercial meshers] Trimmed cell and embedded refinement mesh conversion issues	michele	OpenFOAM Meshing & Mesh Conversion	2	July 15, 2005 04:15
A problem about setting boundary conditions	lyang	Main CFD Forum	0	September 19, 1999 18:29

July 25, 2018, 09:45		#63
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Hi Tar, I am not quite sure what a tangential Laplacian operator is (maybe it is like div(grad(PotE).T() ...?), but your proposed formulae/code look very strange to me. I suggest you examine how others in literature have performed this discretisation, to be clear on what coefficients you expect etc. Philip

July 26, 2018, 10:33		#67
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Ah yes, that makes sense. The projection removes the normal component but then the divergence takes only the normal component. Hmnn actually, we can note the following equivalence: $[(I - nn) \cdot (I - nn)] = (I - nn)$ . And then also we can note that the divergence of the tangential component of the gradient of a scalar is always zero: $\nabla \cdot [(I - nn) \cdot \nabla \phi] = 0$ . So now I realise that I once again don't understand what exactly this tangential Laplacian operator is meant to be :O I guess if you know what you want to do for a simple structured grid, then generalising it to unstructured grids in OpenFOAM should not be an issue if you can explain the behaviour you want for the simple case. Philip

July 27, 2018, 05:44		#68
Kombinator New Member Artem Join Date: Apr 2014 Posts: 29 Rep Power: 12	Dear Philip, Thank you very much for your effort. Anyway, I learned a lot from this topic. I will try to do something and if I succeed I will share my ideas in this topic. Best regards, Tar

December 14, 2021, 06:09		#75
bigphil Super Moderator Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,089 Rep Power: 34	Which version of OpenFOAM are you using? Maybe the "==" operator has been redefined, and the values need to be set using another approach.