CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Error with laplacian of two volScalarFields

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

Like Tree1Likes
  • 1 Post By Hisham

Reply
 
LinkBack Thread Tools Display Modes
Old   December 1, 2011, 13:22
Default Error with laplacian of two volScalarFields
  #1
Senior Member
 
Hisham's Avatar
 
Hisham El Safti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 247
Blog Entries: 10
Rep Power: 8
Hisham is on a distinguished road
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 
      );
solidRho and TwoMuPlusLambda are both volScalarFields .... The code compiles but it gives infinity for initial_D .... If I use a scalar instead of TwoMuPlusLambda, it gives rather logical output.
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
Hisham is offline   Reply With Quote

Old   December 2, 2011, 10:57
Default
  #2
Senior Member
 
Hisham's Avatar
 
Hisham El Safti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 247
Blog Entries: 10
Rep Power: 8
Hisham is on a distinguished road
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;
I applied the solver on the geometry of the plate with a hole tutorial with a fixed value for all patches of a zero vector value.

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.
Hisham is offline   Reply With Quote

Old   December 5, 2011, 06:11
Default
  #3
Senior Member
 
Hisham's Avatar
 
Hisham El Safti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 247
Blog Entries: 10
Rep Power: 8
Hisham is on a distinguished road
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.
Hisham is offline   Reply With Quote

Old   December 5, 2011, 06:46
Default
  #4
Senior Member
 
Hisham's Avatar
 
Hisham El Safti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 247
Blog Entries: 10
Rep Power: 8
Hisham is on a distinguished road
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
Hisham is offline   Reply With Quote

Old   February 12, 2013, 08:04
Default
  #5
Member
 
Hrushi
Join Date: Jan 2013
Posts: 57
Rep Power: 4
hrushi.397 is on a distinguished road
Quote:
Originally Posted by Hisham View Post
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
How exactly did you do the update?
hrushi.397 is offline   Reply With Quote

Old   February 12, 2013, 08:46
Default
  #6
Senior Member
 
Hisham's Avatar
 
Hisham El Safti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 247
Blog Entries: 10
Rep Power: 8
Hisham is on a distinguished road
Quote:
Originally Posted by hrushi.397 View Post
How exactly did you do the update?
Hello

First you do:
Code:
forAll(variable.boundaryField(), patchi)
{
      variable.boundaryField()[patchi] == variable.boundaryField()[patchi].patchInternalField();
}
Then:
Code:
D.correctBoundaryConditions();
Hisham
ChrisA likes this.
Hisham is offline   Reply With Quote

Old   February 12, 2013, 10:07
Default
  #7
Member
 
Hrushi
Join Date: Jan 2013
Posts: 57
Rep Power: 4
hrushi.397 is on a distinguished road
Thanks Hisham . Appreciate your help.
hrushi.397 is offline   Reply With Quote

Old   February 20, 2013, 17:48
Default
  #8
Member
 
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 4
ChrisA is on a distinguished road
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)
ChrisA is offline   Reply With Quote

Old   August 12, 2013, 03:07
Default
  #9
Member
 
Hrushi
Join Date: Jan 2013
Posts: 57
Rep Power: 4
hrushi.397 is on a distinguished road
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
hrushi.397 is offline   Reply With Quote

Old   August 12, 2013, 06:27
Default
  #10
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 565
Rep Power: 19
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by hrushi.397 View Post
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
Hi Hrushi,

Just to clarify, you changed the equation from implicit:
Code:
                fvm::ddt(T) - fvm::laplacian(DT, T)
to explicit like this:
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 1-D heat conduction, the time-step is limited:
dt < rho*c*dx*dx/(2*k)
where dt is the time-step, rho is the density, c is the specific heat, dx is the mesh spacing and k is the conductivity.

Philip
bigphil is offline   Reply With Quote

Old   August 12, 2013, 06:50
Default
  #11
Member
 
Hrushi
Join Date: Jan 2013
Posts: 57
Rep Power: 4
hrushi.397 is on a distinguished road
Thanks for the reply, Phil. I was thinking on similar lines today. I shall definitely try this too.
hrushi.397 is offline   Reply With Quote

Reply

Tags
laplacian, volscalarfield

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 17:45.