Problem with operator fvm::laplacian

palarcon · April 27, 2022, 15:33

Dear all

I'm trying to perform a proof of concept by solving a Poisson equation.
This part is purely ''mathematical'' but because if the computation works I will add it to another OF solver.
The mathematical part inside the code is as follows:

Code:

scalar vVoid 		= 0.1;
vector cosSin(0.7071068,0.7071068,0);
Cxy = vVoid + (1-vVoid)*rho;	// volScalarField
Coh = Cxy * cosSin;				// volVectorField

fvScalarMatrix TarrivalEqn
(
	fvm::laplacian(Coh,Tarrival)
);
TarrivalEqn.solve();

Unfortunately, when I try to compile it I get the following error

Code:

Compiling enabled on 8 cores
Making dependency list for source file overhangFoam.C
g++ -std=c++14 -m64 -pthread -DOPENFOAM=2012 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/usr/lib/openfoam/openfoam2012/src/finiteVolume/lnInclude -I/usr/lib/openfoam/openfoam2012/src/meshTools/lnInclude -iquote. -IlnInclude -I/usr/lib/openfoam/openfoam2012/src/OpenFOAM/lnInclude -I/usr/lib/openfoam/openfoam2012/src/OSspecific/POSIX/lnInclude   -fPIC -c overhangFoam.C -o Make/linux64GccDPInt32Opt/overhangFoam.o
g++ -std=c++14 -m64 -pthread -DOPENFOAM=2012 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/usr/lib/openfoam/openfoam2012/src/finiteVolume/lnInclude -I/usr/lib/openfoam/openfoam2012/src/meshTools/lnInclude -iquote. -IlnInclude -I/usr/lib/openfoam/openfoam2012/src/OpenFOAM/lnInclude -I/usr/lib/openfoam/openfoam2012/src/OSspecific/POSIX/lnInclude   -fPIC -Xlinker --add-needed -Xlinker --no-as-needed  Make/linux64GccDPInt32Opt/overhangFoam.o -L/usr/lib/openfoam/openfoam2012/platforms/linux64GccDPInt32Opt/lib \
    -lfiniteVolume -lfvOptions -lmeshTools -lOpenFOAM -ldl  \
     -lm -o /home/pablo/OpenFOAM/pablo-v2012/platforms/linux64GccDPInt32Opt/bin/overhangFoam
/usr/bin/ld: Make/linux64GccDPInt32Opt/overhangFoam.o: in function `Foam::fv::laplacianScheme<double, Foam::Vector<double> >::New(Foam::fvMesh const&, Foam::Istream&)':
overhangFoam.C:(.text._ZN4Foam2fv15laplacianSchemeIdNS_6VectorIdEEE3NewERKNS_6fvMeshERNS_7IstreamE[_ZN4Foam2fv15laplacianSchemeIdNS_6VectorIdEEE3NewERKNS_6fvMeshERNS_7IstreamE]+0x4b): undefined reference to `Foam::fv::laplacianScheme<double, Foam::Vector<double> >::IstreamConstructorTablePtr_'
/usr/bin/ld: overhangFoam.C:(.text._ZN4Foam2fv15laplacianSchemeIdNS_6VectorIdEEE3NewERKNS_6fvMeshERNS_7IstreamE[_ZN4Foam2fv15laplacianSchemeIdNS_6VectorIdEEE3NewERKNS_6fvMeshERNS_7IstreamE]+0x209): undefined reference to `Foam::fv::laplacianScheme<double, Foam::Vector<double> >::IstreamConstructorTablePtr_'
collect2: error: ld returned 1 exit status
make: *** [/usr/lib/openfoam/openfoam2012/wmake/makefiles/general:150: /home/pablo/OpenFOAM/pablo-v2012/platforms/linux64GccDPInt32Opt/bin/overhangFoam] Error 1

Fields ''rho'', ''Cxy'' qnd ''Tarrival'' are defined as volScalarField with dimensions (Nelem x 1) and field Coh is a volVectorField with dimensions (Nelem x 3).

I know for certain that the problem arises when computing fvm::laplacian, because the types of the fields is not compatible with the operator.

I found a really old discussion that deals with the same problem, but unfortunately I don't understand how can I solve it in my case (original question convert dimensioned<double>’ to ‘double’, answer
convert dimensioned<double>’ to ‘double’)

The mathematical problem I'm trying to solve is in the attached file:

Does anyone knows how to solve it? My gut tells me I'm missing some fundamental stuff, but I can not see it.

Thanks you in advance

snak · April 28, 2022, 08:54

Hi Pablo,

Have you take a look at solidFoam solver?
https://www.openfoam.com/documentati...1d4e7b6b2.html

Laplacian term in anisotropic case will be similar to what you want.
https://www.openfoam.com/documentati...8H_source.html

I guess you need a volSymmTensorField.
https://www.openfoam.com/documentati...ce.html#l00011

Good luck!

dlahaye · April 30, 2022, 08:15

Run scalarTransportFoam with

Code:

ddtSchemes
{
    default         steadyState;
}

palarcon · April 30, 2022, 15:38

Quote:

Originally Posted by dlahaye

Run scalarTransportFoam with

Code:

ddtSchemes
{
    default         steadyState;
}

And how should I define the corresponding fields?
Because scalarTransportFoam uses a velocity field that inside the solver is considered by means of the phi field, but my fields have nothing to do with the velocity or the fluxes... just plain ''mathematical'' equations.
Nonetheless, thanks for the information!

palarcon · April 30, 2022, 15:46

Quote:

Originally Posted by snak

Hi Pablo,

Have you take a look at solidFoam solver?
https://www.openfoam.com/documentati...1d4e7b6b2.html

Laplacian term in anisotropic case will be similar to what you want.
https://www.openfoam.com/documentati...8H_source.html

I guess you need a volSymmTensorField.
https://www.openfoam.com/documentati...ce.html#l00011

Good luck!

I had a look, but couldn't find anything useful unfortunately.

Voulet · May 1, 2022, 15:17

I'm not a profesionnal mathematician but it's only mathematicaly correct if you use a tensor instead of a vector isn't it ?

$\nabla\cdot\left(\lambda\nabla T\right) = \nabla\cdot\begin{pmatrix} \lambda\partial_x T \\ \lambda\partial_y T \\ \lambda\partial_z T \end{pmatrix} = \lambda\left(\partial_x^2 T + \partial_y^2 T + \partial_z^2 T \right)$
= divergence of a vector.

$\nabla\cdot\left(\overline{\overline{\epsilon}}\nabla T\right) = \nabla\cdot\left( \begin{pmatrix} \epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz}\\ \epsilon_{yx} & \epsilon_{yy} & \epsilon_{yz}\\ \epsilon_{zx} & \epsilon_{zy} & \epsilon_{zz} \end{pmatrix} \begin{pmatrix} \partial_x T \\ \partial_y T \\ \partial_z T \end{pmatrix} \right) = \nabla\cdot \begin{pmatrix} \epsilon_{xx}\partial_x T + \epsilon_{xy}\partial_y T + \epsilon_{xz}\partial_z T\\ \epsilon_{yx}\partial_x T + \epsilon_{yy}\partial_y T + \epsilon_{yz}\partial_z T\\ \epsilon_{zx}\partial_x T + \epsilon_{zy}\partial_y T + \epsilon_{zz}\partial_z T \end{pmatrix}$
= divergence of a vector.

Ie in fine you make the divergence of a vector.

Thus you can't do :
$\nabla\cdot\left(\begin{pmatrix} v_x \\ v_y \\ v_z\end{pmatrix}^T\nabla T\right)$
because you then get a scalar inside the divergence.

palarcon · May 2, 2022, 04:27

Quote:

Originally Posted by Voulet

I'm not a profesionnal mathematician but it's only mathematicaly correct if you use a tensor instead of a vector isn't it ?

$\nabla\cdot\left(\lambda\nabla T\right) = \nabla\cdot\begin{pmatrix} \partial_x T \\ \partial_y T \\ \partial_z T \end{pmatrix} = \lambda\left(\partial_x^2 T + \partial_y^2 T + \partial_z^2 T \right)$
= divergence of a vector.

$\nabla\cdot\left(\overline{\overline{\epsilon}}\nabla T\right) = \nabla\cdot\left( \begin{pmatrix} \epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz}\\ \epsilon_{yx} & \epsilon_{yy} & \epsilon_{yz}\\ \epsilon_{zx} & \epsilon_{zy} & \epsilon_{zz} \end{pmatrix} \begin{pmatrix} \partial_x T \\ \partial_y T \\ \partial_z T \end{pmatrix} \right) = \nabla\cdot \begin{pmatrix} \epsilon_{xx}\partial_x T + \epsilon_{xy}\partial_y T + \epsilon_{xz}\partial_z T\\ \epsilon_{yx}\partial_x T + \epsilon_{yy}\partial_y T + \epsilon_{yz}\partial_z T\\ \epsilon_{zx}\partial_x T + \epsilon_{zy}\partial_y T + \epsilon_{zz}\partial_z T \end{pmatrix}$
= divergence of a vector.

Ie in fine you make the divergence of a vector.

Thus you can't do :
$\nabla\cdot\left(\begin{pmatrix} v_x \\ v_y \\ v_z\end{pmatrix}^T\nabla T\right)$
because you then get a scalar inside the divergence.

Thank you very much... your post opened my mind.
The problem was is that I'm trying to implement something that was developed by an student in Matlab, where elementwise operations between vectors or matrices and vectors can be done without any pain.

The solution was to ''simply'' perform the laplacian between a volTensorField (my original matrix, in here Coh but defined properly now) and a volScalarField (Tarrival).

snak · May 2, 2022, 08:49

Hi Pablo,

As you can see the equation by Voulet, Coh will be a volSymmTensorField instead of a volVectorField.

A Laplacian term in solidFoam with anisotropic setting uses volSymmTensorField.

April 30, 2022, 08:15		#3
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 798 Blog Entries: 1 Rep Power: 17	Run scalarTransportFoam with Code: ddtSchemes { default steadyState; }

May 1, 2022, 15:17		#6
Voulet Member Join Date: Jun 2019 Posts: 41 Rep Power: 7	I'm not a profesionnal mathematician but it's only mathematicaly correct if you use a tensor instead of a vector isn't it ? $\nabla\cdot\left(\lambda\nabla T\right) = \nabla\cdot\begin{pmatrix} \lambda\partial_x T \\ \lambda\partial_y T \\ \lambda\partial_z T \end{pmatrix} = \lambda\left(\partial_x^2 T + \partial_y^2 T + \partial_z^2 T \right)$ = divergence of a vector. $\nabla\cdot\left(\overline{\overline{\epsilon}}\nabla T\right) = \nabla\cdot\left( \begin{pmatrix} \epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz}\\ \epsilon_{yx} & \epsilon_{yy} & \epsilon_{yz}\\ \epsilon_{zx} & \epsilon_{zy} & \epsilon_{zz} \end{pmatrix} \begin{pmatrix} \partial_x T \\ \partial_y T \\ \partial_z T \end{pmatrix} \right) = \nabla\cdot \begin{pmatrix} \epsilon_{xx}\partial_x T + \epsilon_{xy}\partial_y T + \epsilon_{xz}\partial_z T\\ \epsilon_{yx}\partial_x T + \epsilon_{yy}\partial_y T + \epsilon_{yz}\partial_z T\\ \epsilon_{zx}\partial_x T + \epsilon_{zy}\partial_y T + \epsilon_{zz}\partial_z T \end{pmatrix}$ = divergence of a vector. Ie in fine you make the divergence of a vector. Thus you can't do : $\nabla\cdot\left(\begin{pmatrix} v_x \\ v_y \\ v_z\end{pmatrix}^T\nabla T\right)$ because you then get a scalar inside the divergence. snak likes this. __________________ « Debugging is what CFD is about. 5 minutes to modify your code, 5 months to find why it does not work anymore. » Last edited by Voulet; May 2, 2022 at 05:42.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
SU2-7.0.1 on ubuntu 18.04	hyunko	SU2 Installation	7	March 16, 2020 05:37
UDF compiling problem	Wouter	Fluent UDF and Scheme Programming	6	June 6, 2012 05:43
Gambit - meshing over airfoil wrapping (?) problem	JFDC	FLUENT	1	July 11, 2011 06:59
natural convection problem for a CHT problem	Se-Hee	CFX	2	June 10, 2007 07:29
Adiabatic and Rotating wall (Convection problem)	ParodDav	CFX	5	April 29, 2007 20:13

April 28, 2022, 08:54		#2
snak Senior Member shinji nakagawa Join Date: Mar 2009 Location: Japan Posts: 113 Blog Entries: 1 Rep Power: 18	Hi Pablo, Have you take a look at solidFoam solver? https://www.openfoam.com/documentati...1d4e7b6b2.html Laplacian term in anisotropic case will be similar to what you want. https://www.openfoam.com/documentati...8H_source.html I guess you need a volSymmTensorField. https://www.openfoam.com/documentati...ce.html#l00011 Good luck!

May 2, 2022, 08:49		#8
snak Senior Member shinji nakagawa Join Date: Mar 2009 Location: Japan Posts: 113 Blog Entries: 1 Rep Power: 18	Hi Pablo, As you can see the equation by Voulet, Coh will be a volSymmTensorField instead of a volVectorField. A Laplacian term in solidFoam with anisotropic setting uses volSymmTensorField.