CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Non-constant Diffusion Coefficient

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes
  • 1 Post By akidess
  • 2 Post By santiagomarquezd

Reply
 
LinkBack Thread Tools Display Modes
Old   November 4, 2010, 03:29
Default Non-constant Diffusion Coefficient
  #1
New Member
 
Join Date: Aug 2010
Posts: 14
Rep Power: 6
NewFoamer is on a distinguished road
Hello,

I'm currently solving the transport equation as follows:

Code:
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:

Code:
Diff = max( upper_value*(1-C) , lower_value )
I'm not quite sure how to go about implementing this. Any help would be appreciated!
NewFoamer is offline   Reply With Quote

Old   November 4, 2010, 09:16
Default
  #2
Senior Member
 
Stefan Herbert
Join Date: Dec 2009
Location: Darmstadt, Germany
Posts: 129
Rep Power: 8
herbert is on a distinguished road
You can hardcode it. All you have to do it to create a volScalarField with a diffusivity-value for each cell.

Code:
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.
herbert is offline   Reply With Quote

Old   November 4, 2010, 11:02
Default
  #3
New Member
 
Join Date: Aug 2010
Posts: 14
Rep Power: 6
NewFoamer is on a distinguished road
Perfect.. it works.. thank you!
NewFoamer is offline   Reply With Quote

Old   November 5, 2010, 09:52
Default
  #4
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
I think the loop can be avoided, Did you try

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

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 17, 2011, 09:30
Default
  #5
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
Santiago, do you know if there's any difference in using the openfoam operations on fields instead of a forAll loop other than readability?
ahmmedshakil likes this.
akidess is offline   Reply With Quote

Old   January 17, 2011, 22:37
Default
  #6
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Anton, good question! Fields operations are quite hard to hack, but I think I got it, max in this case:

Code:
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.
hua1015 and ahmmedshakil like this.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   January 20, 2011, 05:36
Default
  #7
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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 is offline   Reply With Quote

Old   January 28, 2011, 10:38
Default
  #8
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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.
akidess is offline   Reply With Quote

Old   February 1, 2011, 17:14
Default
  #9
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   July 14, 2013, 14:36
Default
  #10
Member
 
Join Date: Jun 2011
Posts: 64
Rep Power: 6
maalan is on a distinguished road
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 is offline   Reply With Quote

Old   July 14, 2013, 15:07
Default
  #11
Member
 
Join Date: Jun 2011
Posts: 64
Rep Power: 6
maalan is on a distinguished road
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 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20
mass flow in is not equal to mass flow out saii CFX 2 September 18, 2009 08:07
Automotive test case vinz OpenFOAM Running, Solving & CFD 98 October 27, 2008 09:43
Two-Phase Buoyant Flow Issue Miguel Baritto CFX 4 August 31, 2006 12:02
Species diffusion coefficient iceabc FLUENT 1 June 10, 2004 10:04


All times are GMT -4. The time now is 11:02.