# Solve poisson equation just add a source term

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 November 3, 2008, 18:32 Hi I am new to OpenFoam, I w #1 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 Hi I am new to OpenFoam, I want to solve a poisson equation of the form laplacian (T) = exp(x)+sin(y) (some source term).I am using the laplacianFoam for this, but I am not able to incorporate the source term in the equation. I have read the programmers guide and it stated that i need to add a volVectorField to carry out this operation. I am not able to comprehend this and could you tell me where exactly I have to add this vectorfield and how to go about coding for the same. Thanks Vishal

 November 6, 2008, 14:19 Hi, I do these stuff in this #2 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 Hi, I do these stuff in this way; first get reference to your x and y position and define S as the source term: volScalarField yPos = mesh.C().component(vector::Y); volScalarField xPos = mesh.C().component(vector::x); volScalarField S ( IOobject ( "S", runTime.timeName(), mesh, IOobject:: NO_READ, IOobject:: NO_WRITE, ), mesh ); Then do whatever you want with your source term for example: forAll(S, counter) { S[counter] = exp(xPos[counter]) + sin(yPos[counter]); } Now solve it: solve(fvm::laplacian(T) == S); This is a little bit messy, someone may do it in a better way. Best, Mahdi

 November 7, 2008, 16:57 Hi, do the following changes #3 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 Hi, do the following changes: volScalarField S ( IOobject ( "S", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); forAll(S,counter) { S[counter] = Foam::sin(M_PI*xPos[counter]) + Foam::sin(M_PI*yPos[counter])); } The first change is to be able to read and write your source term and see if you evaluate it right or wrong. Also to be able to set the dimension. You may also be able to do it in a better way. The second change was to have a source term which is compatible with my boundary conditions. I use zero boudary condition in a 1 by 1 square. Also, I added the "Foam::" because otherwise I cannot compile it. Maybe that's not the case in your version. I use 1.4.1. Also M_PI is pi number. have fun, Mahdi

 November 10, 2008, 17:06 I am also attaching the entire #4 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 I am also attaching the entire code here. Kindly let me know where I am making an error. Regards Vishal

 November 10, 2008, 17:09 #5 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17

 November 10, 2008, 17:13 http://www.cfd-online.com/Ope #6 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 mm.abdollahzadeh likes this.

 November 10, 2008, 18:06 A small change my solver looks #7 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 A small change my solver looks like this.. in the attachment it s wrong,, for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { solve (fvm::poissonian(DT, T)== S); }

 November 10, 2008, 19:19 Hi Mahdi, Thanks a lot for #8 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 Hi Mahdi, Thanks a lot for the statements. 1. You were right I didn't compile them, and hence was getting the same results. 2. Regarding the second comment, shouldn't the object be S instead of T FoamFile { version 2.0; format ascii; class volScalarField; object T;(Shouldn't it be S here) } 3. Regarding the source term factor, then should i have to multiply 1000 to S while solving the equation, that is solve ( fvm::laplacian(DT, T) == 1000 *S (Instead of S) ); Thanks Again for all the help. Regards Vishal

 November 10, 2008, 19:47 Hi, 2. Yeah you are right it #9 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 Hi, 2. Yeah you are right it should be S, but it does not make any difference. 3. You can do it here or in the creatFile.H Enjoy Mahdi

 November 10, 2008, 19:59 Hi Mahdi, Thanks a ton for #10 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 Hi Mahdi, Thanks a ton for all the suggestions. The program now is working perfectly fine. Thanks a lot. Now, I would like to incorporate boundary conditions that are not fixed. Rather they are also functions of x and y. Like my boundary conditions are inlet (AT x =0) - Dirichlet B.C's = -y^3 + sin(y) (something of this type). Similarly (at y=0) (lowerwall) = -x^3 + sin (x).. Neumann B.C's at x =1 (outlet) = -3-2*sin(x)(some arbitary function) Neumann B.C's at y =1 (upperwall) = 4+sin(y) Now, How exactly should I call these in the temperature file (in the 0 folder). I searched the forum (couldn't comprehend what exactly they were saying.) . Can you give your suggestions for these. Regards Vishl

 November 10, 2008, 20:56 Hi, For example you can do t #11 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 Hi, For example you can do this: label inletPatchID = mesh.boundaryMesh().findPatchID("inlet"); fvPatchScalarField& inletT = T.boundaryField()[inletPatchID]; const vectorField& faceCentresInlet = mesh.Cf().boundaryField()[inletPatchID]; forAll(inletT, pointI) { inletT[pointI] = faceCentresInlet[pointI].y()*(1.0 - faceCentresInlet[pointI].y()); } best, Mahdi

 November 10, 2008, 21:27 Hi Mahdi, So Inlet function #12 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 Hi Mahdi, So Inlet function now is Dirichlet B.C given as (y*(1-y)). This is what you have done right. For Gradient (Neumann) B.C how exactly should we go about it, as we have to consider the flux form right for this case. Could you elaborate a little bit regard to this. Regards Vishal

 November 11, 2008, 17:06 Hi, You should add them to y #13 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 Hi, You should add them to your creatField.H file. It does not matter what you set in 0 directory, just set a fixed value of zero for example. Because, then you will modify it by those lines in creatField.H file. About Neumann bounday condition, I never used it. I searched the forum, there is one called wallHeatTransfer, but I have no idea how to set it. Cheers, Mahdi

 November 11, 2008, 17:30 Hi Mahdi, I incorporated th #14 Senior Member   Vishal Nandigana Join Date: Mar 2009 Location: Champaign, Illinois, U.S.A Posts: 208 Rep Power: 17 Hi Mahdi, I incorporated them and tried it again. But this time I got the non uniform values for inlet alone and repeated the same procedure for outlet, lowerwall and upperwall, but didn't see any change. I know I would have made a simple error. Could you point out where exactly I am doing a mistake. I however incorporated the values as a list and tried it out and it is working fine. (I know this is a crude way of doing it). Here is my poissonfoam.c file. #include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" # include "createFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // label inletPatchID = mesh.boundaryMesh().findPatchID("inlet"); fvPatchScalarField& inletT = T.boundaryField()[inletPatchID]; const vectorField& faceCentresInlet = mesh.Cf().boundaryField()[inletPatchID]; forAll(inletT, pointI) { inletT[pointI] = faceCentresInlet[pointI].y()*(1.0 - faceCentresInlet[pointI].y()); } label outletPatchID = mesh.boundaryMesh().findPatchID("outlet"); fvPatchScalarField& outletT = T.boundaryField()[outletPatchID]; const vectorField& faceCentresOutlet = mesh.Cf().boundaryField()[outletPatchID]; forAll(outletT, pointI) { outletT[pointI] = 5.0 - faceCentresOutlet[pointI].y()*(2.0 - faceCentresOutlet[pointI].y()); } label upperwallPatchID = mesh.boundaryMesh().findPatchID("upperwall"); fvPatchScalarField& upperwallT = T.boundaryField()[upperwallPatchID]; const vectorField& faceCentresUpperwall = mesh.Cf().boundaryField()[upperwallPatchID]; forAll(upperwallT, pointIII) { upperwallT[pointIII] = faceCentresUpperwall[pointIII].y()*(4.0 - faceCentresUpperwall[pointIII].y()); } label lowerwallPatchID = mesh.boundaryMesh().findPatchID("lowerwall"); fvPatchScalarField& lowerwallT = T.boundaryField()[lowerwallPatchID]; const vectorField& faceCentresLowerwall = mesh.Cf().boundaryField()[lowerwallPatchID]; forAll(lowerwallT, pointIV) { lowerwallT[pointIV] = (2.0 -faceCentresLowerwall[pointIV].y())*(1.0 - faceCentresLowerwall[pointIV].y()); } T.write(); Info<< "\nCalculating temperature distribution\n" << endl; for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; # include "readSIMPLEControls.H" for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { solve (fvm::laplacian(DT, T)== 1000 * S); } # include "write.H" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return(0); } // ************************************************** *********************** // My T file is giving values for inlet condition but not for the outlet, lowerwall, etc... Its just taking it as zero gradient. Thanks Vishal mm.abdollahzadeh likes this.

 November 11, 2008, 18:04 I cannot find the error #15 Member   M. Mahdi Salehi Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 50 Rep Power: 16 I cannot find the error

 July 14, 2013, 14:21 #16 Member   Join Date: Jun 2011 Posts: 80 Rep Power: 14 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!

 August 28, 2014, 13:02 nandiganavishal and Mahdi #17 New Member   Shakil Masum Join Date: Aug 2014 Posts: 11 Rep Power: 11 Hi nandiganavishal, what did you do to the initial condition for the source "S"? did you add another file in "0" folder named "S"? If yes, is it possible to get that file? Because, i have a constant Source. the problem is I have to create a file in 0 folder named S with internal field and boundary field. since it is a constant value, i do not need to have any information for internal and boundary field but I am not sure about it. Any suggestion would be appreciated. thanks Shakil

January 21, 2017, 04:24
#18
Senior Member

A. Min
Join Date: Mar 2015
Posts: 305
Rep Power: 11
Quote:
 Originally Posted by smehdi609 Hi, do the following changes: volScalarField S ( IOobject ( "S", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); forAll(S,counter) { S[counter] = Foam::sin(M_PI*xPos[counter]) + Foam::sin(M_PI*yPos[counter])); } The first change is to be able to read and write your source term and see if you evaluate it right or wrong. Also to be able to set the dimension. You may also be able to do it in a better way. The second change was to have a source term which is compatible with my boundary conditions. I use zero boudary condition in a 1 by 1 square. Also, I added the "Foam::" because otherwise I cannot compile it. Maybe that's not the case in your version. I use 1.4.1. Also M_PI is pi number. have fun, Mahdi
Hi
could you plz give me procedure to solve this equation:
laplacian(y)=y0
in a rectangular domain that y0 is a velocity distribution { u(x,y) } that is obtained by solving the other problem in this domain and I should import it to this solution.
I saw the probelem of

laplacian(y)=0

in tutorial but I have to change it to laplacian(y)=y0.
How can I do that?
How should I import a u(x,y) as a y0?

 November 14, 2022, 09:12 #19 New Member     Join Date: Feb 2015 Posts: 9 Rep Power: 10 Seems like openfoam can't solve a general Poission Equation with the fvm::laplacian operator. Things are better with the fvc operator, but it is explicit. fvm::laplacian(field) == S gives total garbage.

 Thread Tools Search this Thread Search this Thread: Advanced Search 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 Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post yong Main CFD Forum 2 February 2, 2007 09:31 luckyluke Phoenics 3 December 13, 2004 02:52 Mehdi Main CFD Forum 1 February 17, 2004 23:40 windhair CFX 3 January 27, 2004 21:09 Linfeng BI Main CFD Forum 4 November 20, 2002 23:46

All times are GMT -4. The time now is 01:13.

 Contact Us - CFD Online - Privacy Statement - Top