
[Sponsors] 
November 4, 2010, 03:29 
Nonconstant Diffusion Coefficient

#1 
New Member
Join Date: Aug 2010
Posts: 14
Rep Power: 6 
Hello,
I'm currently solving the transport equation as follows: Code:
solve ( fvm::ddt(C) + fvm::div(phi, C)  fvm::laplacian(Diff, C) ); Code:
Diff = max( upper_value*(1C) , lower_value ) 

November 4, 2010, 09:16 

#2 
Senior Member
Stefan Herbert
Join Date: Dec 2009
Location: Darmstadt, Germany
Posts: 129
Rep Power: 8 
You can hardcode it. All you have to do it to create a volScalarField with a diffusivityvalue 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.0C[cellI]),lowerValue); } solve ( fvm::ddt(C) + fvm::div(phi, C)  fvm::laplacian(Diff, C) ); Stefan P.S. For speed reasons, the declaration of Diff can be done outside the time loop. 

November 4, 2010, 11:02 

#3 
New Member
Join Date: Aug 2010
Posts: 14
Rep Power: 6 
Perfect.. it works.. thank you!


November 5, 2010, 09:52 

#4 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
I think the loop can be avoided, Did you try
Code:
Diff = max(upperValue*(1.0C),lowerValue); Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

January 17, 2011, 09:30 

#5 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17 
Santiago, do you know if there's any difference in using the openfoam operations on fields instead of a forAll loop other than readability?


January 17, 2011, 22:37 

#6 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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++) Hope it might be useful. Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

January 20, 2011, 05:36 

#7 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17 
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.


January 28, 2011, 10:38 

#8 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17 
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 straightforward. 

February 1, 2011, 17:14 

#9 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14 
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. Postdoctoral Fellow Research Center for Computational Mechanics (CIMEC)  CONICET/FICHUNL T.E.: 543424511594 Ext. 1005 Güemes 3450  (3000) Santa Fe Santa Fe  Argentina http://www.cimec.org.ar 

July 14, 2013, 14:36 

#10 
Member
Join Date: Jun 2011
Posts: 64
Rep Power: 6 
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! 

July 14, 2013, 15:07 

#11 
Member
Join Date: Jun 2011
Posts: 64
Rep Power: 6 
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! 

Thread Tools  
Display Modes  


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 
TwoPhase Buoyant Flow Issue  Miguel Baritto  CFX  4  August 31, 2006 12:02 
Species diffusion coefficient  iceabc  FLUENT  1  June 10, 2004 10:04 