
[Sponsors] 
Solve poisson equation just add a source term 

LinkBack  Thread Tools  Search this Thread  Display Modes 
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.cfdonline.com/Ope

#6 
Senior Member
Vishal Nandigana
Join Date: Mar 2009
Location: Champaign, Illinois, U.S.A
Posts: 208
Rep Power: 17 

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) = 32*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*(1y)). 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 

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:
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 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
How to solve for transport with large source term  yong  Main CFD Forum  2  February 2, 2007 09:31 
How to add a source term in u or v equation?  luckyluke  Phoenics  3  December 13, 2004 02:52 
Can we solve only convective and source term  Mehdi  Main CFD Forum  1  February 17, 2004 23:40 
How to solve Poisson equation in CFX 4?  windhair  CFX  3  January 27, 2004 21:09 
Solve Poisson equation  Linfeng BI  Main CFD Forum  4  November 20, 2002 23:46 