CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Solving Coupled Set of PDEs

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

Like Tree5Likes
  • 4 Post By Tron
  • 1 Post By yhaomin2007

Reply
 
LinkBack Thread Tools Display Modes
Old   May 24, 2011, 15:16
Default Solving Coupled Set of PDEs
  #1
New Member
 
Nick Mai
Join Date: Feb 2011
Posts: 5
Rep Power: 6
nickmai123 is on a distinguished road
Hi Everyone,

I'm trying to solve the following coupled system with OpenFOAM:

A) \partial_t u=\nabla^2 v
B) \partial_t v=-\nabla^2 u

I've realized that you can use matrices to represent the switch and declare a tensor as follows:

Code:
tensor j(0,1,0,-1,0,0,0,0,0);
solve
{
   fvm::ddt(U)
 ==
   j & fvm::laplacian(U)
}
However, the inner product is not defined between a primitive tensor and a fvMatrix. Does anyone here have any recommendations? I need to use the fvm:: namespace because it does implicit discretization of the Laplacian.

Thanks for reading!

EDIT: I've found coupledFvScalarMatrix in the OpenFOAM 1.6-ext version. However, from the examples that I see, it seems as though it is intended for multiple domains/meshes. Is it possible to use the coupledFvScalarMatrix for a single mesh? If so, then my problem is solved.

--Nick
nickmai123 is offline   Reply With Quote

Old   May 24, 2011, 16:05
Default try blockCoupledScalarTransportFoam
  #2
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18
chegdan will become famous soon enough
What about the blockCoupledScalarTransportFoam? I know that is on the same mesh, with two equations coupled to eachother through a source term


Code:
fvScalarMatrix TEqn
            (
                fvm::div(phi, T)
              - fvm::laplacian(DT, T)
             ==
                 alpha*Ts
              - fvm::Sp(alpha, T)
            );

            TEqn.relax();

            fvScalarMatrix TsEqn
            (
              - fvm::laplacian(DTs, Ts)
             ==
                 alpha*T
              - fvm::Sp(alpha, Ts)
            );
With some additional methods to form a block matrix. Maybe that will help.

Dan
chegdan is offline   Reply With Quote

Old   May 25, 2011, 02:15
Default
  #3
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Hi,

there is a coupled solver under development in -extend, as pointed out in comment 2 (warning: it is a bit complicated at a first glance).

However, you could define j as a volTensorField, and initialise it with your value of j.

Code:
 volTensorField j
 (
    IOobject
    (
        "j",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedTensor("jTensor", dimless, Foam::Tensor(0,1,0,-1,0,0,0,0,0))
 );
where "dimless" should become

dimensionSet(.....)

if the tensor has dimensional units. Keep in mind however that the solution won't be coupled but sequential: meaning component by component.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; May 25, 2011 at 02:16. Reason: Added comment
alberto is offline   Reply With Quote

Old   November 24, 2011, 09:09
Default blockCoupledScalarTransportFoam - coupling via operator
  #4
New Member
 
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8
Tron is on a distinguished road
Hi!

I am trying to implement a solver for the following two coupled equations:

Equation 1: laplacian(D1, phi_real) - laplacian(D2, phi_img) = 0
Equation 2: laplacian(D1, phi_img) + laplacian(D2, phi_real) = 0

First, I implemented a segregated version:
solve( fvm::laplacian(D1, phi_real) - fvc::laplacian(D2, phi_img));
solve( fvm::laplacian(D1, phi_img) + fvc::laplacian(D2, phi_real));
which seems to work for weakly couplings. Naturally, this approach is rather useless if the coupling grows stronger.

Hence, I had a look at 1.6.-ext's solver "blockCoupledScalarTransportFoam". I guess I understand the main idea of it. I think one of the critical parts is the manipulation of the off-diagonal elements and the block source in order to create the correctly coupled system which looks like that in blockCoupledScalarTransportFoam.C:

forAll(d, i)
{
d[i](0,1) = -alpha.value()*mesh.V()[i];
d[i](1,0) = -alpha.value()*mesh.V()[i];

blockB[i][0] -= alpha.value()*blockX[i][1]*mesh.V()[i];
blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i];
}

The above approach is rather simple for a coupling of the equations via the scalar "alpha". However, can something similar also be done if the coupling is done via an operator such as laplacian()? I tried to figure that out, but I am stuck. Any helpful comments would be very appreciated.

Cheers,
Christian
Tron is offline   Reply With Quote

Old   July 26, 2012, 10:41
Default
  #5
New Member
 
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5
tiat is on a distinguished road
hi,

is there anybody still following this topic?

i'm having exactly the same interests as christian has pointed out. The blockCoupledScalarTransportFoam is very understandable for the idea of block matrix construction. However, when it comes to the coupling of transport terms (like grad(p) in the momentum equation, not as simple as the source term coupling in the blockCoupledScalarTransportFoam case), it is hard to construct the block because no implict gradient scheme is available

does anyone have idea about a fvm::grad() scheme?
tiat is offline   Reply With Quote

Old   July 27, 2012, 01:20
Default
  #6
New Member
 
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8
Tron is on a distinguished road
Hi Tian,

Quote:
Originally Posted by tiat View Post

is there anybody still following this topic?
Yes, I do. At least I did. I managed to solve my problem. If my memory servers me right then I wrote a solver that could solve the following coupled equations
grad * (sigma * grad(phi_real)) - grad * ((E0*Er*w) * grad(phi_imag)) = 0
grad * (sigma * grad(phi_imag)) + grad * ((E0*Er*w) * grad(phi_real)) = 0

using the following code:

#include "fvCFD.H"
#include "fieldTypes.H"
#include "Time.H"
#include "fvMesh.H"

#include "blockLduSolvers.H"

#include "VectorNFieldTypes.H"
#include "volVectorNFields.H"
#include "blockVectorNMatrices.H"

#include "blockMatrixTools.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


int main(int argc, char *argv[])

{

# include "setRootCase.H"

# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


// # include "CourantNo.H"
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readSIMPLEControls.H"


for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

{
fvScalarMatrix phi_realEqn
(
- fvm::laplacian(EpsilonTimesAngVel, phi_real)
);

phi_realEqn.relax();


fvScalarMatrix phi_realEqn2

(
- fvm::laplacian(Sigma, phi_real)
);

phi_realEqn2.relax();


fvScalarMatrix phi_imgEqn

(
- fvm::laplacian(EpsilonTimesAngVel, phi_img)
);

phi_imgEqn.relax();


fvScalarMatrix phi_imgEqn2

(
- fvm::laplacian(Sigma, phi_img)
);

phi_imgEqn2.relax();


// Prepare block system

BlockLduMatrix<vector2> blockM(mesh);
BlockLduMatrix<vector2> blockM2(mesh);

// Grab block diagonal and set it to zero

Field<tensor2>& d = blockM.diag().asSquare();
d = tensor2::zero;
Field<tensor2>& d2 = blockM2.diag().asSquare();
d2 = tensor2::zero;

// Grab linear off-diagonal and set it to zero

Field<tensor2>& u = blockM.upper().asSquare();
Field<tensor2>& l = blockM.lower().asSquare();
u = tensor2::zero;
l = tensor2::zero;
Field<tensor2>& u2 = blockM2.upper().asSquare();
Field<tensor2>& l2 = blockM2.lower().asSquare();
u2 = tensor2::zero;
l2 = tensor2::zero;

// forAll(u, j)

// {
// u[j](0,0) = 0.0;
// u[j](1,0) = 0.0;
// u[j](0,1) = 0.0;
// u[j](1,1) = 0.0;
// u2[j](0,0) = 0.0;
// u2[j](1,0) = 0.0;
// u2[j](0,1) = 0.0;
// u2[j](1,1) = 0.0;
// }

// forAll(l, j)
// {
// l[j](0,0) = 0.0;
// l[j](1,0) = 0.0;
// l[j](0,1) = 0.0;
// l[j](1,1) = 0.0;
// l2[j](0,0) = 0.0;
// l2[j](1,0) = 0.0;
// l2[j](0,1) = 0.0;
// l2[j](1,1) = 0.0;
// }

vector2Field& blockX = blockT.internalField();

vector2Field& blockX2 = blockT2.internalField();
vector2Field blockB(mesh.nCells(), vector2::zero);
vector2Field blockB2(mesh.nCells(), vector2::zero);

//- Inset equations into block Matrix

blockMatrixTools::insertEquation(0, phi_realEqn, blockM, blockX, blockB);
blockMatrixTools::insertEquation(1, phi_imgEqn, blockM, blockX, blockB);
blockMatrixTools::insertEquation(0, phi_realEqn2, blockM2, blockX2, blockB2);
blockMatrixTools::insertEquation(1, phi_imgEqn2, blockM2, blockX2, blockB2);

forAll(d, i)

{
d[i](0,1) = d2[i](1,1);
d[i](1,0) = -d2[i](0,0);
blockB[i][0] += blockB2[i][1];
blockB[i][1] -= blockB2[i][0];
}

forAll(u, j )

{
u[j](0,1) = u2[j](1,1);
u[j](1,0) = -u2[j](0,0);
}

forAll(l, j )

{
l[j](0,1) = l2[j](1,1);
l[j](1,0) = -l2[j](0,0);
}

//- Block coupled solver call

BlockSolverPerformance<vector2> solverPerf =
BlockLduSolver<vector2>::New
(
word("blockVar"),
blockM,
mesh.solver("blockVar")
)->solve(blockX, blockB);

solverPerf.print();


// Retrieve solution

blockMatrixTools::blockRetrieve(0, phi_real.internalField(), blockX);
blockMatrixTools::blockRetrieve(1, phi_img.internalField(), blockX);

phi_real.correctBoundaryConditions();

phi_img.correctBoundaryConditions();
}

volVectorField gradPhi_real = fvc::grad(phi_real);

volVectorField gradPhi_img = fvc::grad(phi_img);

# include "write.H"
}

Info<< "End\n" << endl;


return(0);

}


I do not remember any details anymore. However, I remember that the solver worked nicely.

I hope this helps,
Christian
Tron is offline   Reply With Quote

Old   July 27, 2012, 02:39
Default
  #7
New Member
 
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5
tiat is on a distinguished road
hi, christian,

thanks soooo much for your kind reply!
i understand your code generally, but i just have a question, in my equation, the coupling term is the gradient of p (in your equation, the coupling is the laplacian, and luckily OpenFOAM has the fvm::laplacian(), but they don't have fvm::grad()...that's my big problem!).

do you also have any idea about this?
tiat is offline   Reply With Quote

Old   July 27, 2012, 10:10
Default
  #8
New Member
 
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8
Tron is on a distinguished road
Quote:
Originally Posted by tiat View Post
do you also have any idea about this?
Unfortunately not. It is quite a while ago and I didn't work on this project since then. Sorry.
Tron is offline   Reply With Quote

Old   July 29, 2012, 04:53
Default
  #9
New Member
 
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5
tiat is on a distinguished road
hi, chritian,

thanks all the same

Tian
tiat is offline   Reply With Quote

Old   September 5, 2012, 16:05
Default
  #10
Member
 
Haomin Yuan
Join Date: Jan 2012
Location: Madison, Wisconsin, USA
Posts: 54
Rep Power: 5
yhaomin2007 is on a distinguished road
HI, All

As I read your guys' reply, you are couple two equations maximum. Is there any example that we can couple more equations? like 3 or 4?

thank you in advance
Annier likes this.
yhaomin2007 is offline   Reply With Quote

Old   March 27, 2014, 10:05
Default
  #11
New Member
 
Join Date: Mar 2014
Posts: 20
Rep Power: 3
alisina-s is on a distinguished road
Quote:
Originally Posted by Tron View Post
Hi Tian,


Yes, I do. At least I did. I managed to solve my problem. If my memory servers me right then I wrote a solver that could solve the following coupled equations
grad * (sigma * grad(phi_real)) - grad * ((E0*Er*w) * grad(phi_imag)) = 0
grad * (sigma * grad(phi_imag)) + grad * ((E0*Er*w) * grad(phi_real)) = 0

using the following code:

#include "fvCFD.H"
#include "fieldTypes.H"
#include "Time.H"
#include "fvMesh.H"

#include "blockLduSolvers.H"

#include "VectorNFieldTypes.H"
#include "volVectorNFields.H"
#include "blockVectorNMatrices.H"

#include "blockMatrixTools.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


int main(int argc, char *argv[])

{

# include "setRootCase.H"

# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


// # include "CourantNo.H"
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readSIMPLEControls.H"


for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

{
fvScalarMatrix phi_realEqn
(
- fvm::laplacian(EpsilonTimesAngVel, phi_real)
);

phi_realEqn.relax();


fvScalarMatrix phi_realEqn2

(
- fvm::laplacian(Sigma, phi_real)
);

phi_realEqn2.relax();


fvScalarMatrix phi_imgEqn

(
- fvm::laplacian(EpsilonTimesAngVel, phi_img)
);

phi_imgEqn.relax();


fvScalarMatrix phi_imgEqn2

(
- fvm::laplacian(Sigma, phi_img)
);

phi_imgEqn2.relax();


// Prepare block system

BlockLduMatrix<vector2> blockM(mesh);
BlockLduMatrix<vector2> blockM2(mesh);

// Grab block diagonal and set it to zero

Field<tensor2>& d = blockM.diag().asSquare();
d = tensor2::zero;
Field<tensor2>& d2 = blockM2.diag().asSquare();
d2 = tensor2::zero;

// Grab linear off-diagonal and set it to zero

Field<tensor2>& u = blockM.upper().asSquare();
Field<tensor2>& l = blockM.lower().asSquare();
u = tensor2::zero;
l = tensor2::zero;
Field<tensor2>& u2 = blockM2.upper().asSquare();
Field<tensor2>& l2 = blockM2.lower().asSquare();
u2 = tensor2::zero;
l2 = tensor2::zero;

// forAll(u, j)

// {
// u[j](0,0) = 0.0;
// u[j](1,0) = 0.0;
// u[j](0,1) = 0.0;
// u[j](1,1) = 0.0;
// u2[j](0,0) = 0.0;
// u2[j](1,0) = 0.0;
// u2[j](0,1) = 0.0;
// u2[j](1,1) = 0.0;
// }

// forAll(l, j)
// {
// l[j](0,0) = 0.0;
// l[j](1,0) = 0.0;
// l[j](0,1) = 0.0;
// l[j](1,1) = 0.0;
// l2[j](0,0) = 0.0;
// l2[j](1,0) = 0.0;
// l2[j](0,1) = 0.0;
// l2[j](1,1) = 0.0;
// }

vector2Field& blockX = blockT.internalField();

vector2Field& blockX2 = blockT2.internalField();
vector2Field blockB(mesh.nCells(), vector2::zero);
vector2Field blockB2(mesh.nCells(), vector2::zero);

//- Inset equations into block Matrix

blockMatrixTools::insertEquation(0, phi_realEqn, blockM, blockX, blockB);
blockMatrixTools::insertEquation(1, phi_imgEqn, blockM, blockX, blockB);
blockMatrixTools::insertEquation(0, phi_realEqn2, blockM2, blockX2, blockB2);
blockMatrixTools::insertEquation(1, phi_imgEqn2, blockM2, blockX2, blockB2);

forAll(d, i)

{
d[i](0,1) = d2[i](1,1);
d[i](1,0) = -d2[i](0,0);
blockB[i][0] += blockB2[i][1];
blockB[i][1] -= blockB2[i][0];
}

forAll(u, j )

{
u[j](0,1) = u2[j](1,1);
u[j](1,0) = -u2[j](0,0);
}

forAll(l, j )

{
l[j](0,1) = l2[j](1,1);
l[j](1,0) = -l2[j](0,0);
}

//- Block coupled solver call

BlockSolverPerformance<vector2> solverPerf =
BlockLduSolver<vector2>::New
(
word("blockVar"),
blockM,
mesh.solver("blockVar")
)->solve(blockX, blockB);

solverPerf.print();


// Retrieve solution

blockMatrixTools::blockRetrieve(0, phi_real.internalField(), blockX);
blockMatrixTools::blockRetrieve(1, phi_img.internalField(), blockX);

phi_real.correctBoundaryConditions();

phi_img.correctBoundaryConditions();
}

volVectorField gradPhi_real = fvc::grad(phi_real);

volVectorField gradPhi_img = fvc::grad(phi_img);

# include "write.H"
}

Info<< "End\n" << endl;


return(0);

}


I do not remember any details anymore. However, I remember that the solver worked nicely.

I hope this helps,
Christian
Hi christan
In which version of openfoam you compiled your solver?
I want to solve some coupled transport equation and some one suggested me the blockCoupledScalarTransportFoam.C that is similar to your solver and I know this just work in extended version 3 and higher and not in version 2.2
alisina-s is offline   Reply With Quote

Old   March 28, 2014, 03:18
Default used version
  #12
New Member
 
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8
Tron is on a distinguished road
Hello!

Well, that was a long long time ago when I worked on this solver I just checked my source code archive: I used OpenFoam-1.6-ext for this project.

Christian
Tron is offline   Reply With Quote

Old   March 30, 2014, 14:19
Default
  #13
New Member
 
Join Date: Mar 2014
Posts: 20
Rep Power: 3
alisina-s is on a distinguished road
Thanks christian
alisina-s is offline   Reply With Quote

Old   April 6, 2014, 08:35
Default
  #14
New Member
 
Join Date: Mar 2014
Posts: 20
Rep Power: 3
alisina-s is on a distinguished road
Quote:
Originally Posted by alberto View Post
Hi,

there is a coupled solver under development in -extend, as pointed out in comment 2 (warning: it is a bit complicated at a first glance).

However, you could define j as a volTensorField, and initialise it with your value of j.

Code:
 volTensorField j
 (
    IOobject
    (
        "j",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedTensor("jTensor", dimless, Foam::Tensor(0,1,0,-1,0,0,0,0,0))
 );
where "dimless" should become

dimensionSet(.....)

if the tensor has dimensional units. Keep in mind however that the solution won't be coupled but sequential: meaning component by component.

Best,
Hi
I have the same problem somehow
I want to solve
dT/dt +d(Uk)/dx =0
dk/dt +d(UT)/dx =0
that are coupled
and I have to use some tensors also
thanks for the reply but in
dimensionedTensor("jTensor", dimless, Foam::Tensor(0,1,0,-1,0,0,0,0,0)) what is jTensor?
because i have seen "0" somewhere
I sthe "0" means 0 directory?
alisina-s is offline   Reply With Quote

Old   April 24, 2014, 08:01
Default
  #15
New Member
 
Join Date: Mar 2014
Posts: 20
Rep Power: 3
alisina-s is on a distinguished road
Does anyone still following this topic?
I want to solve a coupled equation
dT/dt + d/dx(Uk) = 0
dk/dt + d/dx(UT) = 0
for solving that
volTensorField j
(
IOobject
(
"j",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedTensor( "jTensor", dimless,tensor(0,-1,0,1,0,0,0,0,0) )
);
and solving
fvVectorMatrix TEqn
(
fvm::ddt(f)
==
( j & fvc::div(phi,f) )
);

TEqn.solve();
runTime.write();
and defined f as


volVectorField f
(
IOobject
(
"f",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh

);

forAll(f,I)
{
f[I].component(0)=T[I];
f[I].component(1)=kappa[I];
f[I].component(2)=0;
}

my code compiled correctly but my results are unexpected
I am in doubt with f,I mean I am not sure what kind of boundary conditions I have to implement for a vectorfield that including two scalarfield(in this case T and kappa)
alisina-s is offline   Reply With Quote

Reply

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
icoLagrangianFoam OF1.6 myNewParticleSolver heavy_user OpenFOAM 16 February 11, 2012 06:15
Error while running rhoPisoFoam.. nileshjrane OpenFOAM Running, Solving & CFD 8 August 26, 2010 12:50
ForcesCoeffs ronaldo OpenFOAM 4 September 14, 2009 07:11
Problems with simulating TurbFOAM barath.ezhilan OpenFOAM 13 July 16, 2009 05:55
MRFSimpleFoam amp cyclic patches david OpenFOAM Running, Solving & CFD 36 October 21, 2008 21:55


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