CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Problem with operator fvm::laplacian (https://www.cfd-online.com/Forums/openfoam-programming-development/242560-problem-operator-fvm-laplacian.html)

palarcon April 27, 2022 14:33
Problem with operator fvm::laplacian
 
1 Attachment(s)
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 https://www.cfd-online.com/Forums/op...tml#post467411, answer
https://www.cfd-online.com/Forums/op...tml#post468691)

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 07: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 07:15
Run scalarTransportFoam with

Code:
ddtSchemes
{
    default        steadyState;
}

palarcon April 30, 2022 14:38
Quote:
Originally Posted by dlahaye (Post 827227)
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 14:46
Quote:
Originally Posted by snak (Post 827111)
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 14: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 03:27
Quote:
Originally Posted by Voulet (Post 827271)
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 07: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.


All times are GMT -4. The time now is 09:33.