Gradient operator implicit discretization

diegon · November 29, 2006, 06:54

In order to build a solver for combustion I need to discretize the radiative transfer equation.
The RTE contains the term:

s*grad(I)

where I is a scalar and s a vector.
My problem comes from the gradient operator for which is avaible only the explicit representation but I need an implicit form to solve for I.
So I have thought to use some math and rewrite it in a form containg operator having an implicit discretization such as divergence:

s*grad(I)=div(s*I)-I*div(s)

Should it work?

hartinger · November 29, 2006, 08:32

Yep, looks correct, should work.

surfaceScalarField sf = fvc::surfaceInterpolate(s)& mesh.Sf();
& mesh.Sf()
s * grad(I) := fvm::div(sf, I) - fvm::Sp(fvc::div(fs), I);

the 'Sp' - thing adds a coefficient to the diagonal of the implicit matrix.

Taking this opportunity, why is there no implicit grad implementation? Anybody?

pierre and markus

hartinger · November 29, 2006, 08:34

ignore second & mesh.Sf()

P & M

hjasak · November 29, 2006, 08:50

Implicit gradient operator:

- firstly, the diagonal would be zero.
- secondly, the matrix coefficients would be vectors for a gradient and vectors transpose for a divergence
- thirdly, you cannot solve the equation

grad(thingy) = rhs

beucase the diagonal of the gradient matrix equals zero for a uniform mesh

Implicit gradient matrix makes sense only for implicit block coupled (e.g. pressure velocity) algorithms, and I'm pretty sure noone is quite there yet with OpenFOAM.

Hrv

diegon · November 29, 2006, 09:54

The equation I need to discretize is not

s*grad(I)=div(s*I)-I*div(s)

but it contains the "s*grad(I)" that I have thought to sobstitute it with "div(s*I)-I*div(s)".

hartinger · November 29, 2006, 10:24

Hrv, thanks, does make sense.

Diego, we decribed the implementation of "s*grad(I)" as "div(s*I)-I*div(s)", which would be one of the terms for the matrix setup (:= means defined as)

fvScalarMatrix yourEqn
(
...
+ fvm::div(sf, I) - fvm::Sp(fvc::div(fs), I)
...
);

PM

diegon · November 29, 2006, 10:35

Ok thanks

so it seems it could not work.

hartinger · November 29, 2006, 10:48

yes, it can

you can't have the term "s*grad(I)" implicitly, but you can replace that with the term you suggested "div(s*I)-I*div(s)", for which we gave the actual implementation.
The "fvm::"-prefix means in Foam-speak implicit. More precise, it is the "fvm" namespace in which all implicit functions for the Finite Volume Method (fvm) are defined.

"fvc::" denotes "Finite Volume Calculus", all explicit stuff.

So again, your reasoning is right, you can do it as you suggested.

P & M

diegon · November 29, 2006, 10:55

Sorry but I did not get what Jasak was writing so I guessed it would not have worked.

Thank you again.

matteoc · December 27, 2006, 07:23

Hi Diego,

I think that "s" is a const vector, once u have decided the direction of the radiation...
so div(s) must be equal to zero.

So, I think u can write:
s&grad(I) = div(sI)

bye
M

Erik · May 2, 2007, 06:52

Hi

This might be a dumb question but is this why the pressure is solved in a semi-discretised form of the momentum equation (A and H decompositions and solving through Jacobi metod) i.e. to find another way of implementing an implicit form of grad(p)?

/Erik

hjasak · May 2, 2007, 06:58

Do you know CFX? They implement a pressure-based block solver and they indeed have an implicit grad (and div!) to form a 2x2 block matrix system.

No such thing in OpenFOAM at the moment.

Hrv

Erik · May 2, 2007, 07:22

Thank you Hrv!

I guess this is the reason for the special treatment of grad(p) then.

I dont know about CFX. I am fairly new to the field of CFD and keen on using and learning more about OpenFOAM.

My problem is that I am trying to implement a different momentum equation involving gradients of density as well as the gradient of pressure. I would like to know how to formulate this in a similar manner to the one done in the PISO-loop. Ive looked through your Ph.D and found some information on the subject but I would like to see some DOC (if available) on how to get the momentum equation in the semi-discretised form. Is such DOC available to your knowledge?

regards
/Erik

MdoNascimento · November 5, 2015, 04:51

Quote:

Originally Posted by hjasak

Implicit gradient operator:

- firstly, the diagonal would be zero.
- secondly, the matrix coefficients would be vectors for a gradient and vectors transpose for a divergence
- thirdly, you cannot solve the equation

grad(thingy) = rhs

beucase the diagonal of the gradient matrix equals zero for a uniform mesh

Hello, could anyone explain these arguments more detailed?

Thanks

allett02015 · October 9, 2017, 13:08

Even if it has passed quite a few time since this post I'll answer anyway.

If you discretize your gradient with a finite Volume Method and sum the values over the faces of your control volume, the value at your control volume center cancels out for a regular grid.

Novel · November 2, 2017, 05:50

I came across the same problem after discretizing of the gradient of an scalar.

The gradient of a scalar is discretized as

$\nabla \phi = \frac{1}{V_P} \sum\limits^{n_{faces}}_{f} \phi_f \mathbf{S}_f = RHS$

Where the lower f indicates the value of the scalar at the faces

$\phi_f = \omega_f \phi_P + (1-\omega_P)\phi_F$

F represents the neighboring nodes

For an weight of 0.5

$\phi_f = \frac{\phi_P+\phi_F}{2}$

So the gradient can be written as

$\nabla \phi = \frac{S}{V_P} ( \frac{\phi_P + \phi_N}{2}\mathbf{n}_n + \frac{\phi_P + \phi_E}{2}\mathbf{n}_e + \frac{\phi_P + \phi_S}{2}\mathbf{n}_s + \frac{\phi_P + \phi_W}{2}\mathbf{n}_w + ... ) = RHS$

For a structured grid we can assume that the normal vector of one face must be the normal vector of another face just in the opposite direction

$\mathbf{n}_n = -\mathbf{n}_s$

We can see that in the discretized form the center point of the owner vanishes

$\nabla \phi = \frac{S}{V_P} ( \frac{\phi_N - \phi_S}{2}\mathbf{n}_n + \frac{\phi_E - \phi_W}{2}\mathbf{n}_e + ... ) = RHS$

The obtained matrix wont have any diagonal elements.

mboesi · July 22, 2020, 06:59

I stumbled across this thread while searching for a way to discretize a somehow similar problem. I need to discretize the following term: $c\nabla\cdot\vec{U}$

Rearranging the divergence term: $\nabla\cdot\left(c\vec{U}\right)$ gives results in the following discretization: $c\nabla\cdot\vec{U} = \nabla\cdot\left(c\vec{U}\right) - \nabla c\cdot\vec{U}$

The implementation of the first term of the RHS should be:

Code:

fvm::div(phi,c)

Based on the example above, I've tried to implement the second term of the RHS as follows:

Code:

 fvm::Sp(fvc::grad(c),phi)

However, the second term doesn't compile, since there is no suitable Sp method available for the supplied arguments. I'm missing something on this term.

Has anybody of you an idea on how to correctly discretize RHS?

CFD_Freemountain · February 25, 2022, 13:26

Hello mboesi,

did you already solve your problem?

I'm facing the same issue.

With kind regards
Christian

mAlletto · February 27, 2022, 10:02

Quote:

Originally Posted by mboesi

I stumbled across this thread while searching for a way to discretize a somehow similar problem. I need to discretize the following term: $c\nabla\cdot\vec{U}$

Rearranging the divergence term: $\nabla\cdot\left(c\vec{U}\right)$ gives results in the following discretization: $c\nabla\cdot\vec{U} = \nabla\cdot\left(c\vec{U}\right) - \nabla c\cdot\vec{U}$

The implementation of the first term of the RHS should be:

Code:

fvm::div(phi,c)

Based on the example above, I've tried to implement the second term of the RHS as follows:

Code:

 fvm::Sp(fvc::grad(c),phi)

However, the second term doesn't compile, since there is no suitable Sp method available for the supplied arguments. I'm missing something on this term.

Has anybody of you an idea on how to correctly discretize RHS?

I think here the problem is that phi is a surfaceScajarfield and not a volScalarField. The second argument of the fvm:: Sp operator should always be the variable you are solving for and the variables are defined at the cell centers. The OpenFOAM propramer guide gives some advice how to use it

FengShawn · August 27, 2022, 11:44

has anybody solved this problem？

November 29, 2006, 06:54	In order to build a solver for	#1
diegon New Member diego n. Join Date: Mar 2009 Posts: 17 Rep Power: 17	In order to build a solver for combustion I need to discretize the radiative transfer equation. The RTE contains the term: sgrad(I) where I is a scalar and s a vector. My problem comes from the gradient operator for which is avaible only the explicit representation but I need an implicit form to solve for I. So I have thought to use some math and rewrite it in a form containg operator having an implicit discretization such as divergence: sgrad(I)=div(sI)-Idiv(s) Should it work? Kummi likes this.

November 29, 2006, 08:32	Yep, looks correct, should wor	#2
hartinger Senior Member Markus Hartinger Join Date: Mar 2009 Posts: 102 Rep Power: 17	Yep, looks correct, should work. surfaceScalarField sf = fvc::surfaceInterpolate(s)& mesh.Sf(); & mesh.Sf() s * grad(I) := fvm::div(sf, I) - fvm::Sp(fvc::div(fs), I); the 'Sp' - thing adds a coefficient to the diagonal of the implicit matrix. Taking this opportunity, why is there no implicit grad implementation? Anybody? pierre and markus babala and Kummi like this.

November 29, 2006, 08:34	ignore second & mesh.Sf()	#3
hartinger Senior Member Markus Hartinger Join Date: Mar 2009 Posts: 102 Rep Power: 17	ignore second & mesh.Sf() P & M

November 29, 2006, 08:50	Implicit gradient operator:	#4
hjasak Senior Member Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,907 Rep Power: 33	Implicit gradient operator: - firstly, the diagonal would be zero. - secondly, the matrix coefficients would be vectors for a gradient and vectors transpose for a divergence - thirdly, you cannot solve the equation grad(thingy) = rhs beucase the diagonal of the gradient matrix equals zero for a uniform mesh Implicit gradient matrix makes sense only for implicit block coupled (e.g. pressure velocity) algorithms, and I'm pretty sure noone is quite there yet with OpenFOAM. Hrv chegdan and Teresa.Z like this. __________________ Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk

November 29, 2006, 09:54	The equation I need to discret	#5
diegon New Member diego n. Join Date: Mar 2009 Posts: 17 Rep Power: 17	The equation I need to discretize is not sgrad(I)=div(sI)-Idiv(s) but it contains the "sgrad(I)" that I have thought to sobstitute it with "div(sI)-Idiv(s)".

November 29, 2006, 10:24	Hrv, thanks, does make sense.	#6
hartinger Senior Member Markus Hartinger Join Date: Mar 2009 Posts: 102 Rep Power: 17	Hrv, thanks, does make sense. Diego, we decribed the implementation of "sgrad(I)" as "div(sI)-I*div(s)", which would be one of the terms for the matrix setup (:= means defined as) fvScalarMatrix yourEqn ( ... + fvm::div(sf, I) - fvm::Sp(fvc::div(fs), I) ... ); PM

November 29, 2006, 10:35	Ok thanks so it seems it c	#7
diegon New Member diego n. Join Date: Mar 2009 Posts: 17 Rep Power: 17	Ok thanks so it seems it could not work.

November 29, 2006, 10:48	yes, it can you can't have	#8
hartinger Senior Member Markus Hartinger Join Date: Mar 2009 Posts: 102 Rep Power: 17	yes, it can you can't have the term "sgrad(I)" implicitly, but you can replace that with the term you suggested "div(sI)-I*div(s)", for which we gave the actual implementation. The "fvm::"-prefix means in Foam-speak implicit. More precise, it is the "fvm" namespace in which all implicit functions for the Finite Volume Method (fvm) are defined. "fvc::" denotes "Finite Volume Calculus", all explicit stuff. So again, your reasoning is right, you can do it as you suggested. P & M Kummi likes this.

November 29, 2006, 10:55	Sorry but I did not get what J	#9
diegon New Member diego n. Join Date: Mar 2009 Posts: 17 Rep Power: 17	Sorry but I did not get what Jasak was writing so I guessed it would not have worked. Thank you again.

December 27, 2006, 07:23	Hi Diego, I think that "s"	#10
matteoc New Member matteo cerutti Join Date: Mar 2009 Location: Florence, Tuscany, Italy Posts: 10 Rep Power: 17	Hi Diego, I think that "s" is a const vector, once u have decided the direction of the radiation... so div(s) must be equal to zero. So, I think u can write: s&grad(I) = div(sI) bye M

May 2, 2007, 06:52	Hi This might be a dumb que	#11
Erik Member Erik Arlemark Join Date: Mar 2009 Location: Eindhoven, Netherlands Posts: 47 Rep Power: 17	Hi This might be a dumb question but is this why the pressure is solved in a semi-discretised form of the momentum equation (A and H decompositions and solving through Jacobi metod) i.e. to find another way of implementing an implicit form of grad(p)? /Erik

May 2, 2007, 06:58	Do you know CFX? They impleme	#12
hjasak Senior Member Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,907 Rep Power: 33	Do you know CFX? They implement a pressure-based block solver and they indeed have an implicit grad (and div!) to form a 2x2 block matrix system. No such thing in OpenFOAM at the moment. Hrv __________________ Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk

May 2, 2007, 07:22	Thank you Hrv! I guess this	#13
Erik Member Erik Arlemark Join Date: Mar 2009 Location: Eindhoven, Netherlands Posts: 47 Rep Power: 17	Thank you Hrv! I guess this is the reason for the special treatment of grad(p) then. I dont know about CFX. I am fairly new to the field of CFD and keen on using and learning more about OpenFOAM. My problem is that I am trying to implement a different momentum equation involving gradients of density as well as the gradient of pressure. I would like to know how to formulate this in a similar manner to the one done in the PISO-loop. Ive looked through your Ph.D and found some information on the subject but I would like to see some DOC (if available) on how to get the momentum equation in the semi-discretised form. Is such DOC available to your knowledge? regards /Erik

October 9, 2017, 13:08		#15
allett02015 New Member Michael Join Date: Feb 2015 Posts: 18 Rep Power: 11	Even if it has passed quite a few time since this post I'll answer anyway. If you discretize your gradient with a finite Volume Method and sum the values over the faces of your control volume, the value at your control volume center cancels out for a regular grid.

November 2, 2017, 05:50		#16
Novel New Member Roman G. Join Date: Apr 2017 Posts: 16 Rep Power: 9	I came across the same problem after discretizing of the gradient of an scalar. The gradient of a scalar is discretized as $\nabla \phi = \frac{1}{V_P} \sum\limits^{n_{faces}}_{f} \phi_f \mathbf{S}_f = RHS$ Where the lower f indicates the value of the scalar at the faces $\phi_f = \omega_f \phi_P + (1-\omega_P)\phi_F$ F represents the neighboring nodes For an weight of 0.5 $\phi_f = \frac{\phi_P+\phi_F}{2}$ So the gradient can be written as $\nabla \phi = \frac{S}{V_P} ( \frac{\phi_P + \phi_N}{2}\mathbf{n}_n + \frac{\phi_P + \phi_E}{2}\mathbf{n}_e + \frac{\phi_P + \phi_S}{2}\mathbf{n}_s + \frac{\phi_P + \phi_W}{2}\mathbf{n}_w + ... ) = RHS$ For a structured grid we can assume that the normal vector of one face must be the normal vector of another face just in the opposite direction $\mathbf{n}_n = -\mathbf{n}_s$ We can see that in the discretized form the center point of the owner vanishes $\nabla \phi = \frac{S}{V_P} ( \frac{\phi_N - \phi_S}{2}\mathbf{n}_n + \frac{\phi_E - \phi_W}{2}\mathbf{n}_e + ... ) = RHS$ The obtained matrix wont have any diagonal elements.

August 27, 2022, 11:44		#20
FengShawn New Member FengShawn Join Date: Jul 2022 Location: USA Posts: 1 Rep Power: 0	has anybody solved this problem？ Last edited by FengShawn; August 31, 2022 at 12:49.

July 22, 2020, 06:59		#17
mboesi New Member Markus Bösenhofer Join Date: Mar 2016 Posts: 2 Rep Power: 0	I stumbled across this thread while searching for a way to discretize a somehow similar problem. I need to discretize the following term: $c\nabla\cdot\vec{U}$ Rearranging the divergence term: $\nabla\cdot\left(c\vec{U}\right)$ gives results in the following discretization: $c\nabla\cdot\vec{U} = \nabla\cdot\left(c\vec{U}\right) - \nabla c\cdot\vec{U}$ The implementation of the first term of the RHS should be: Code: fvm::div(phi,c) Based on the example above, I've tried to implement the second term of the RHS as follows: Code: fvm::Sp(fvc::grad(c),phi) However, the second term doesn't compile, since there is no suitable Sp method available for the supplied arguments. I'm missing something on this term. Has anybody of you an idea on how to correctly discretize RHS?

February 25, 2022, 13:26		#18
CFD_Freemountain New Member Christian Ernst Join Date: Oct 2021 Posts: 1 Rep Power: 0	Hello mboesi, did you already solve your problem? I'm facing the same issue. With kind regards Christian

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
FvPatchScalarField access operator	maka	OpenFOAM	3	July 1, 2008 09:32
Implicit Formulation of Velocity Gradient	stefan82	OpenFOAM Running, Solving & CFD	1	August 9, 2007 12:09
Operator Splitting.	Maria.	Main CFD Forum	5	September 17, 2005 23:10
Operator precedence	hemph	OpenFOAM	1	September 13, 2005 13:40
a math operator in UDF	lichun Dong	FLUENT	7	June 18, 2005 23:04