
[Sponsors] 
December 1, 2011, 13:22 
Error with laplacian of two volScalarFields

#1 
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 249
Blog Entries: 10
Rep Power: 9 
Dear Foamers,
I have this equation: Code:
const volScalarField TwoMuPlusLambda = 2 * solidShearMod + solidLambdaMod; //solve fvVectorMatrix gEqn ( fvm::d2dt2(solidRho, initial_D) == fvm::laplacian( TwoMuPlusLambda, initial_D, "laplacian(DD,D)") + divSigmaExp ); It seems there is something wrong with laplacian for a volScalarField Gamma ... Can anyone give me feedback on that. I use OF 2.0.x on Ubuntu 11.04 Best regards, Hisham 

December 2, 2011, 10:57 

#2 
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 249
Blog Entries: 10
Rep Power: 9 
I tried to compare two equations one with a scalar and another with a volScalarField for comparisons and they yield different results for the same numerical input ... Results for the scalar are the logical ones:
Code:
volVectorField D ( IOobject ( "D", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "\nReading field initial_D\n" << endl; // Initial displacement volVectorField initial_D ( IOobject ( "initial_D", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); volSymmTensorField sigmaD ( IOobject ( "sigmaD", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), solidShearMod * twoSymm(fvc::grad(D)) + solidLambdaMod * (I*tr(fvc::grad(D))) ); volSymmTensorField initialSigmaD ( IOobject ( "initialSigmaD", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), sigmaD ); volVectorField divSigmaExp ( IOobject ( "divSigmaExp", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::div(sigmaD) ); Info<< "Initializing field of explicit plastic div(sigma) divSigmaExpPlastic\n" << endl; volVectorField divSigmaExpPlastic ( IOobject ( "divSigmaExpPlastic", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::div(sigmaD) ); volTensorField gradInitialD = fvc::grad(initial_D); volTensorField gradD = fvc::grad(D); Info << "gradD " << gradD << endl; Info << "gradInitialD " << gradInitialD << endl; initialSigmaD = solidShearMod * twoSymm(gradInitialD) + solidLambdaMod * (I * tr(gradInitialD)); sigmaD = mu * twoSymm(gradD) + lambda * (I * tr(gradD)) ; Info << "sigmaD " << sigmaD << endl; Info << "initialSigmaD " << initialSigmaD << endl; divSigmaExp = fvc::div ( initialSigmaD  (2 * solidShearMod + solidLambdaMod ) * gradInitialD, "div(sigmaD)" ) + solidRho * Gravity; Info << "divSigmaExp " << divSigmaExp << endl; divSigmaExpPlastic = fvc::div ( sigmaD  (2 * mu + lambda ) * gradD, "div(sigmaD)" ) + rho * Gravity; Info << "divSigmaExpPlastic " << divSigmaExpPlastic << endl; // + solidRho * Gravity ; const volScalarField TwoMuPlusLambda = 2 * solidShearMod + solidLambdaMod; //solve fvVectorMatrix gEqn ( fvm::d2dt2(solidRho, initial_D) == fvm::laplacian( 2* solidShearMod + solidLambdaMod, initial_D, "laplacian(DD,D)") + divSigmaExp ); fvVectorMatrix eEqn ( fvm::d2dt2(rho, D) == fvm::laplacian( 2* mu + lambda, D, "laplacian(DD,D)") + divSigmaExpPlastic ); initialResidual = eEqn.solve().initialResidual(); initialResidual = gEqn.solve().initialResidual(); Info << "D " << D << endl; Info << "initial_D " << initial_D << endl; The results of "D" are logical but for "initial_D" they are not .... the numbers are the same until the gEq and eEq are first solved. 

December 5, 2011, 06:11 

#3 
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 249
Blog Entries: 10
Rep Power: 9 
To summarize:
When I do a laplacian of a volVectorField multiplied by a constant scalar it runs OK and values are logical But when I replace the constant scalar with a constant volScalarField with a uniform value of the same value as the scalar in the former case, it compiles but gives very large values for the volVectorField solved for that paraFoam views as "inf" or very high values when I change the values. I need a volScalarField to introduce nonuniform values in practice. Can someone please help me? It seems like a bug, but I want to make sure I'm not doing it wrong before sending a bug report. 

December 5, 2011, 06:46 

#4 
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 249
Blog Entries: 10
Rep Power: 9 
I found out what the problem was:
The values of the volScalarField at patches needed to be updated to be equal to the values of cells with faces at patch. I tried that and it gave same results as a the scalar. I do not know if this fully solved the problem or not .... Until other bugs appear, I would assume it debugged 

February 12, 2013, 08:04 

#5  
Member
Hrushi
Join Date: Jan 2013
Posts: 58
Rep Power: 6 
Quote:


February 12, 2013, 10:07 

#7 
Member
Hrushi
Join Date: Jan 2013
Posts: 58
Rep Power: 6 
Thanks Hisham . Appreciate your help.


February 20, 2013, 17:48 

#8 
Member
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 6 
Thank you very much for following up with the solution, I never would have figured out my issue otherwise. (I didn't even think of this as an issue)


August 12, 2013, 03:07 

#9 
Member
Hrushi
Join Date: Jan 2013
Posts: 58
Rep Power: 6 
Hi Hisham,
I am facing a somewhat similar issue now. It has to do with laplacian. I am trying to play around with laplacianFoam and I tried to change the equation from implicit to explicit(fvm::laplacian(DT,T) to fvc::laplacian(DT,T)) just to see if I get a bit different answer. To my surprise the laplacianFoam crashed. I checked the laplacian values and it seems laplacian explodes and runs into very high values. What could have gone wrong? Is this a bug in Openfoam? Thanks, Hrushi 

August 12, 2013, 06:27 

#10  
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 629
Rep Power: 23 
Quote:
Just to clarify, you changed the equation from implicit: Code:
fvm::ddt(T)  fvm::laplacian(DT, T) Code:
fvm::ddt(T)  fvc::laplacian(DT, T) If you did this then, this method will be conditionally stable (as with all explicit methods). Therefore the results will 'blow up' unless the stability criterion is satisfied. For 1D heat conduction, the timestep is limited: dt < rho*c*dx*dx/(2*k) where dt is the timestep, rho is the density, c is the specific heat, dx is the mesh spacing and k is the conductivity. Philip 

August 12, 2013, 06:50 

#11 
Member
Hrushi
Join Date: Jan 2013
Posts: 58
Rep Power: 6 
Thanks for the reply, Phil. I was thinking on similar lines today. I shall definitely try this too.


Tags 
laplacian, volscalarfield 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
How to calculate laplacian of a scalar in cfx?  Jun  CFX  16  October 24, 2013 17:21 
Modifying the laplacian operator  mlawson  OpenFOAM Running, Solving & CFD  11  September 7, 2011 12:30 
laplacian problem  ziemowitzima  OpenFOAM Running, Solving & CFD  0  August 20, 2011 16:56 
laplacian  nimasam  OpenFOAM  7  May 9, 2011 16:06 
why laplacian() failed  OFCrash  OpenFOAM Running, Solving & CFD  1  February 1, 2010 08:32 