[NewFoamer]

Hello,

I'm currently solving the transport equation as follows:

solve
(
fvm::ddt(C)
+ fvm::div(phi, C)
- fvm::laplacian(Diff, C)
);

That is, assuming a constant diffusion coefficient. In order to improve the accuracy of my simulation, I however have to let the diffusion coefficient vary with the concentration according to following:

Diff = max( upper_value*(1-C) , lower_value )

I'm not quite sure how to go about implementing this. Any help would be appreciated!

[herbert]

You can hardcode it. All you have to do it to create a volScalarField with a diffusivity-value for each cell.

volScalarField Diff
(
    IOobject
    (
        "Diff",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("Diff",dimensionSet(0,2,-1,0,0),0.0)
);

forAll(Diff, cellI)
{
    Diff[cellI] = max(upperValue*(1.0-C[cellI]),lowerValue);
}

solve
(
    fvm::ddt(C)
  + fvm::div(phi, C)
  - fvm::laplacian(Diff, C)
);

Regards,
Stefan

P.S. For speed reasons, the declaration of Diff can be done outside the time loop.

[NewFoamer]

Perfect.. it works.. thank you!

[santiagomarquezd]

I think the loop can be avoided, Did you try

Diff = max(upperValue*(1.0-C),lowerValue);

?

Regards.

[akidess]

Santiago, do you know if there's any difference in using the openfoam operations on fields instead of a forAll loop other than readability?

[santiagomarquezd]

Anton, good question! Fields operations are quite hard to hack, but I think I got it, max in this case:

GeometricFieldFunctions.C

00551 BINARY_TYPE_FUNCTION(Type, Type, Type, max)

FieldFunctions.C

00535 BINARY_TYPE_FUNCTION(Type, Type, Type, max)

DimensionedFieldFunctions.C

00633     Foam::OpFunc(tRes().field(), df1.field(), dt2.value());                   \

FieldM.H

00188 // member function : this f1 OP fUNC f2, s
00189 
00190 #define TFOR_ALL_F_OP_FUNC_F_S(typeF1, f1, OP, FUNC, typeF2, f2, typeS, s)  \
00191                                                                             \
00192     /* check the two fields have same Field<Type> mesh */                   \
00193     checkFields(f1, f2, "f1 " #OP " " #FUNC "(f2, s)");                     \
00194                                                                             \
00195     /* set access to f1, f2 and f3 at end of each field */                  \
00196     List_ACCESS(typeF1, f1, f1P);                                           \
00197     List_CONST_ACCESS(typeF2, f2, f2P);                                     \
00198                                                                             \
00199     /* loop through fields performing f1 OP1 f2 OP2 f3 */                   \
00200     List_FOR_ALL(f1, i)                                                     \
00201         List_ELEM(f1, f1P, i) OP FUNC(List_ELEM(f2, f2P, i), (s));          \
00202     List_END_FOR_ALL


ListLoopM.H

00049 #define List_ACCESS(type, f, fp) \
00050     type* const __restrict__ fp = (f).begin()


00052 #define List_CONST_ACCESS(type, f, fp) \
00053     const type* const __restrict__ fp = (f).begin()
00054 
00055 #else


00059 #define List_FOR_ALL(f, i)                      \
00060         register label i = (f).size();          \
00061         while (i--)                             \
00062         {                                       \
00063 
00064 #define List_END_FOR_ALL  }

00066 #define List_ELEM(f, fp, i)  (*fp++)

there is a lot of macros inside and it's difficult to walk into them with gdb. I think the main advantage of using these ops is that they are programmed to copy from memory to memory in a quasi raw form, i.e. copying from list to list and due the macros code is practically hardcoded at compilation time. If you use the forAll loop it is necessary to access to classes dynamically, do the dereferences, etc. I think it's slower.

Hope it might be useful.

Regards.

[akidess]

Santiago, thanks a lot for your effort, I appreciate it. I'll try to do a profiling run to see if I can confirm your findings.

[akidess]

Santiago, your assumption was correct! I tested this on my own solver and the result was astounding. By replacing two forAll loops in my code by the top level notation, the computation got 3x faster.

My initial motivation to use forAll was higher optimization - following code should result in a loop to calculate A, and a second one to calculate D:
A = B * C;
D = (E + A) * (F + A);

With a forAll loop you can do the computation by looping through the fields only once, yet it turns out now it's slower anyway.

The OpFunc function you posted makes me think of expression templates, but to be honest I don't fully understand it yet. As you mentioned hacking through all the macros is not straight-forward.

[santiagomarquezd]

Anton, thanks for sharing your experience. These procedures with macros, etc. are well optimized for up two ops at the same time. In your case you have four ops and they are nested, preprocessing requires to separate the loops and I think it cannot be inlined at all, nevertheless, as you showed, field ops are still faster in this case. I think it isn't the rule ever and it is necessary to do some profiling to assure it.

Regards.

[maalan]

Hi there!!

I am trying to solve the laplace's equation for a scalar phi over an arbitrary domain!! do you know how I should modify the laplacianFoam application?? I have noticed in the .C file the laplacian application has 2 inputs and I just need one of them: laplacian(phi). By the moment, the steps I have followed have been:

1) create my phi volScalarField in the same manner as T in the original laplacianFoam, and comment the rest

2) leave the original libraries in the top of the laplacianFoam.C file and modify the diffusion equation by

solve
(
fvm::laplacian(phi);
);

It should be very easy, but I don't see the mistake!
Thanks a lot!

[maalan]

Hi there!!

I am trying to solve the laplace's equation for a scalar phi over an arbitrary domain!! do you know how I should modify the laplacianFoam application?? I have noticed in the .C file the laplacian application has 2 inputs and I just need one of them: laplacian(phi). By the moment, the steps I have followed have been:

1) create my phi volScalarField in the same manner as T in the original laplacianFoam, and comment the rest

2) leave the original libraries in the top of the laplacianFoam.C file and modify the diffusion equation by

solve
(
fvm::laplacian(phi);
);

It should be very easy, but I don't see the mistake!
Thanks a lot!