[edmund32]

Hi,
sorry, I feel like I'm just doing this incorrectly, but I'm trying to add a description of the magnetic field to one of the compressible solvers, and I have modified the energy equation like so:

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      - 2*DMU*fvc::div(B*(B & U))
      + (DDB*fvc::div(fvc::div(B*B)))
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
      - DDB*fvm::laplacian(sqr(B))
     ==
        fvOptions(rho, he)
    );

Now, DDB is a dimensionedScalar, he a volScalarField, U a volVectorField, and B a volVectorField. The line in particular that causes issue (code compiles fine without it) is:

      - DDB*fvm::laplacian(sqr(B))

and it throws the issue:

In file included from mhdFoam.C:149:0:
EEqn.H: In function 'int main(int, char**)':
EEqn.H:21:34: error: no matching function for call to 'laplacian(Foam::tmp<Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh> >)'
       - DDB*fvm::laplacian(sqr(B))
                                  ^
In file included from /home/edders/OpenFOAM-dev/src/finiteVolume/lnInclude/fvmLaplacian.H:188:0,
                 from /home/edders/OpenFOAM-dev/src/finiteVolume/lnInclude/fvm.H:46,
                 from /home/edders/OpenFOAM-dev/src/finiteVolume/lnInclude/fvCFD.H:10,
                 from mhdFoam.C:36:
/home/edders/OpenFOAM-dev/src/finiteVolume/lnInclude/fvmLaplacian.C:45:1: note: candidate: template<class Type> Foam::tmp<Foam::fvMatrix<Type> > Foam::fvm::laplacian(const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&, const Foam::word&)
 laplacian
 ^

and so on... I have tried all feasible changes to this operator ( DDB*fvm::laplacian(sqr(B)), DDB*fvm::laplacian(B*B)), DDB*fvm::laplacian(B&B)) and nothing seems to work. What am I doing wrong here?
Thanks

[Tobi]

Hi,

is your B field really a volVectorField?

EEqn.H:21:34: error: no matching function for call to 'laplacian(Foam::tmp<Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh> >)'

As you can see, the B field is a symmTensor field (if I get it correct).

[edmund32]

Hi, yeah the B field is defined in createFields.H as

Info<< "Reading field B\n" << endl;
volVectorField B
(
    IOobject
    (
        "B",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

But I fixed my problem anyway, I guess because that equation is defined as fvScalarMatrix, it required the terms in the Laplacian to be volScalarField, and the dot product of the B field didn't produce that type for some reason. So i defined a new scalar field above it:

volScalarField BB
(
    IOobject
    (
        "BB",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    B&B
);

and then used it in the equation:

fvScalarMatrix EEqn
(
    fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
- 2*DMU*fvc::div(B*(B & U))
+ (DDB*fvc::div(fvc::div(B*B)))
+ (
        he.name() == "e"
    ? fvc::div
        (
            fvc::absolute(phi/fvc::interpolate(rho), U),
            (p+DMU*(BB)),
            "div(phiv,p)"
        )
    : -dpdt
    )
- fvm::laplacian(turbulence->alphaEff(), he)
- DDB*fvc::laplacian((BB))
==
    fvOptions(rho, he)
);

and now it seems to work.

[edmund32]

also yes, I think the output from the Laplacian, when the input was a volVectorField, is a symmTensor field

[Gerry Kan]

Howdy Folks:

At the risk of hijacking the thread ......

I am trying to perform a Laplacian of a volVectorField U that contains only the Uz component (i.e., Ux and Uy will be set to zero). It should be implemented in UEqn.H like this:

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + turbulence->divDevRhoReff(U)
 ==
    fvOptions(rho, U)
  + fvm::laplacian(rho*50.0,Uz)
);

Never mind the 50. I will spend the effort to code this properly if this works.

To produce the Uz volVectorField, I tried the following method:

volVectorField& Uz = U.component(vector::Z)

where I got a compiler error of "invalid initialization of non-const reference ... from from an rvalue ...". Since U is not a constant, I am a bit surprised I got this error.

Then I tried the following:

volVectorField Uz(U.component(vector::Z), mesh)

This time it compiles, but at runtime the solver is looking for a file called "U.component(2)" in my case directory, which does not exist.

I am not quite sure how I should proceed. Any pointer would be greatly appreciated.

Sincerely, Gerry.

[edmund32]

Hi, yeah it's looking for that file because it want you to set boundary conditions for the new vector field you've defined, so you can either add a Uz file to the ./0/ folder, with the boundary conditions taken from U, or you can define this vector field in a similar way to how I defined the BB field above, i.e. with NO_READ defined so that it doesn't look for BCs

[Tobi]

Originally Posted by edmund32

Hi, yeah it's looking for that file because it want you to set boundary conditions for the new vector field you've defined, so you can either add a Uz file to the ./0/ folder, with the boundary conditions taken from U, or you can define this vector field in a similar way to how I defined the BB field above, i.e. with NO_READ defined so that it doesn't look for BCs

Hi Edmund, I guess your statement is incorrect in terms of the idea of Gerry.

The first way of his implementation is not working because a volvectorField is not a scalar Field. The return of U.x() does return the x component as a scalarField (out of the box, check doxygen).

The second one calls the constructor of Geometric Field and uses the provided information. However, if he uses the second approach and a new file in the 0 folder needs to be created as you already pointed out (never mind the name of the file). However, he does have again the same problem, if I get his question correct and also the code snippet.

I guess he wants the velocity of z only to be diffused while the other components should only feel the convection and time term't. No idea for which purpose one would use this but this would make at least a bit sense to me, otherwise the equation is somehow wired (and I am sure it does not work with fvm)

This assumption comes from the fact that he tries to reference the U.x() component rather than copying it or just use the second way you provided.

The first question we should know is, is my assumption correct? If not, use the approach given by Edmund and the explicit Laplace operator.

[Tobi]

Originally Posted by edmund32

the BB field above, i.e. with NO_READ defined so that it doesn't look for BCs

This is wrong. Setting No_Read probably needs the boundaries. In your case it is taken from the fields you set here. If you don't use the B field here you are forced to give information about the boundaries. You just had luck.

[Gerry Kan]

Originally Posted by Tobi

I guess he wants the velocity of z only to be diffused while the other components should only feel the convection and time term't. No idea for which purpose one would use this but this would make at least a bit sense to me, otherwise the equation is somehow wired (and I am sure it does not work with fvm).

Dear Tobias:

You guessed correct. I would like to introduce a numerical diffusion term to the z-Momentum equation only. I use a compressible solver and I am seeing spurious oscillation in Uz near the corresponding boundaries. I ruled out pressure-momentum coupling, as well as advection scheme. My next guess is the sysyem's sensitivity to boundary conditions. This is ehy I want to see if adding the diffusion term in z alone, at least at the start of the simulation, would dampen the oscillation sufficiently.

I have posted my original question here for reference.

Gerry.

[Gerry Kan]

Now that you mentioned .... you're both right. fvm::lsplacian won't work on U in UEwn because it must be operated on a volScalarField, not a volVectorField. The closest. thing I can think of is something similar to the turbulent stress tensor.

Having said that, there is a part of me that tells me the offending velocity oscillation could be remedied without artificial diffusion; I just don't know how exactly.

Gerry.

[Tobi]

Originally Posted by Gerry Kan

Now that you mentioned .... you're both right. fvm::lsplacian won't work on U in UEwn because it must be operated on a volScalarField, not a volVectorField.

No, the problem here is that you solve for the field U. You build a matrix for that quantity and while you add a new variable as fvm:: (implicit), the field goes to the matrix diagonal elements which cannot be because here, you have the values of the cells of the field U and not a other field. Thus, you need to use the explicit form (fvc::).

However, I guess the easiest is to use a vector as diffusivity rather than modifying the field here:

// Correct the dimensions, should be m²/s I guess
const dimensionedVector diffusivity("diffusivity", dimensionSet(0,0,0,0,0,0,0), vector(50,0,0));

// In your equation
fvm::laplacian(diffusivity, U);

https://cpp.openfoam.org/v7/namespac...b464d008b3fb12

[Gerry Kan]

Hallo Tobias:

I suppose I can try first with the entire U first. If I get the desired damping effects in Uz, I can think about isolating the Laplacian for the z-momentum only.

Thanks, Gerry.

[Gerry Kan]

Dear Tobias:

Now I have a bigger screen, I managed to read your response more carefully.

Actually, declaring a diffusivity vector and setting the other two components to zero will serve the purpose of introducing one-dimension diffusion without having to reconstruct a new volVectorField. Thank you very much for this.

Further, based on this thread, I realize that it is still possible to call a Laplacian of U inside UEqn, which is still good.

I implemented it quickly to test, and the code compiled. However, there was a linker error

undefined reference to
`Foam::fv::laplacianScheme<
   Foam::Vector<double>, 
   Foam::Vector<double>>::IstreamConstructorTablePtr_'

What I don't quite get, is why all in a sudden U is recognized as a vector<double>, when I was expecting (at least from createFields.H) that it should be a volVectorField. Any comments would be appreciated.

Thanks, Gerry.

P.S. - Here is my implementation, in case there is a typo that I did not catch ...

// in the main foam solver C file before the time loop

const dimensionedVector GAMMA  (
   "GAMMA",
   dimensionSet ( 1, -1, -1, 0, 0, 0, 0 ),
   vector ( 0.0, 0.0, 1.0 )
);

// in UEqn.H

tmp<fvVectorMatrix> tUEqn  (
   fvm::ddt(rho,U) + fvm::div(phi,U) +
   turbulence->divDevRhoReff(U)
   ==
   fvOptions(rho,U) +
   fvm::laplacian(GAMMA,U)
);

[Gerry Kan]

It looks like the first argument of the laplacian operator is to be a volScalarField, at least based on what I saw from scouring the OpenFOAM source files. Changing the definition of the diffusivity term accordingly solved linker issue for now.

At the moment I am looking at how the artificial damping will affect the velocity field ... so far the solver runs stably. I also noticed that this also affects the p and T transport processes, possibly because this artificial diffusion coefficient being set a bit too high. The job is still running as we speak ... let's see how it looks like.

Gerry.