# Non-constant Diffusion Coefficient

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

 November 4, 2010, 03:29 Non-constant Diffusion Coefficient #1 New Member   Join Date: Aug 2010 Posts: 14 Rep Power: 9 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!

 November 4, 2010, 09:16 #2 Senior Member   Stefan Herbert Join Date: Dec 2009 Location: Darmstadt, Germany Posts: 129 Rep Power: 10 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.

 November 4, 2010, 11:02 #3 New Member   Join Date: Aug 2010 Posts: 14 Rep Power: 9 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: 438 Rep Power: 17 I think the loop can be avoided, Did you try Code: `Diff = max(upperValue*(1.0-C),lowerValue);` ? Regards. arvindpj likes this. __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC) - CONICET/UNL Tel: 54-342-4511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe - Argentina. http://www.cimec.org.ar

 January 17, 2011, 09:30 #5 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,264 Rep Power: 23 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.

 January 17, 2011, 22:37 #6 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 438 Rep Power: 17 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 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. Research Scientist Research Center for Computational Methods (CIMEC) - CONICET/UNL Tel: 54-342-4511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe - Argentina. http://www.cimec.org.ar

 January 20, 2011, 05:36 #7 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,264 Rep Power: 23 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: Germany Posts: 1,264 Rep Power: 23 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.

 February 1, 2011, 17:14 #9 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 438 Rep Power: 17 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. Research Scientist Research Center for Computational Methods (CIMEC) - CONICET/UNL Tel: 54-342-4511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe - Argentina. http://www.cimec.org.ar

 July 14, 2013, 14:36 #10 Member   Join Date: Jun 2011 Posts: 79 Rep Power: 8 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: 79 Rep Power: 8 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 Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20 saii CFX 2 September 18, 2009 08:07 vinz OpenFOAM Running, Solving & CFD 98 October 27, 2008 09:43 Miguel Baritto CFX 4 August 31, 2006 12:02 iceabc FLUENT 1 June 10, 2004 10:04

All times are GMT -4. The time now is 04:28.