[wernero]

Hi all!
I'm upgrading a solver I wrote (which is numerically good, but not enough practical), by changing some condition in the case, and adding another equation for conservation in the solver.
the solver is for Maxwell inductance equation. The equation for conservation is given by div(Ject)=0. Here Jext is current density.

the solver did run good before adding the divergence part. the whole code until the div part is below:

#include "fvCFD.H"


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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"


    
    #include "createFields.H"


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

    Info<< nl << "Starting time loop" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        
            fvVectorMatrix AEqn
            (
                fvm::ddt(A)
              - lambda*fvm::laplacian(A)
              - (U^(fvc::curl(A)))
              ==
                Jext/sigma
            );
        // solve the induction equation
          AEqn.solve();
          
          surfaceScalarField phi
                     (
                         "phi",
                         (fvc::interpolate(Jext) & mesh.Sf())
                     );

          fvScalarMatrix JextEqn
          (
              fvc::div(phi)==0
          );
          JextEqn.solve();

when trying to use wmake, the error below show up:

                 from inductionAC2_Foam.C:51:
/home/orw/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:647:13: note: template<class Form, class Cmpt, int nCmpt> bool Foam::operator==(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, const Foam::VectorSpace<Form, Cmpt, nCmpt>&)
 inline bool operator==
             ^
/home/orw/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:647:13: note:   template argument deduction/substitution failed:
inductionAC2_Foam.C:94:30: note:   ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ is not derived from ‘const Foam::VectorSpace<Form, Cmpt, nCmpt>’
               fvc::div(phi)==0
                              ^
make: *** [Make/linux64GccDPInt64Opt/inductionAC2_Foam.o] Error 1

the error is actually endless, so I hope it could be recognized from it's end, above.

I whould really appreciate any help. I already search the internet, saw the conservation term fvc::div(phi,u), but that is not the case for me. Couldn't find any explanation to fix my code.
thanks,

Or Werner.

[blais.bruno]

I am not 100% sure, but I think that fvc:: operators are explicit operator, that is they do not generate any matrix system.

Only fvm:: operators generate a matrix system. Therefore you are getting an error because you are trying to solve something which is not a system.

fvc::div(J)==0 does not mean anything because this is not an equation.
If I were you, i would try to rewrite mathematically your problem to get a laplacian equation for J, juste like we get in traditional CFD

BB

Originally Posted by wernero

Hi all!
I'm upgrading a solver I wrote (which is numerically good, but not enough practical), by changing some condition in the case, and adding another equation for conservation in the solver.
the solver is for Maxwell inductance equation. The equation for conservation is given by div(Ject)=0. Here Jext is current density.

the solver did run good before adding the divergence part. the whole code until the div part is below:

#include "fvCFD.H"


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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"


    
    #include "createFields.H"


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

    Info<< nl << "Starting time loop" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        
            fvVectorMatrix AEqn
            (
                fvm::ddt(A)
              - lambda*fvm::laplacian(A)
              - (U^(fvc::curl(A)))
              ==
                Jext/sigma
            );
        // solve the induction equation
          AEqn.solve();
          
          surfaceScalarField phi
                     (
                         "phi",
                         (fvc::interpolate(Jext) & mesh.Sf())
                     );

          fvScalarMatrix JextEqn
          (
              fvc::div(phi)==0
          );
          JextEqn.solve();

when trying to use wmake, the error below show up:

                 from inductionAC2_Foam.C:51:
/home/orw/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:647:13: note: template<class Form, class Cmpt, int nCmpt> bool Foam::operator==(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, const Foam::VectorSpace<Form, Cmpt, nCmpt>&)
 inline bool operator==
             ^
/home/orw/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/VectorSpaceI.H:647:13: note:   template argument deduction/substitution failed:
inductionAC2_Foam.C:94:30: note:   ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ is not derived from ‘const Foam::VectorSpace<Form, Cmpt, nCmpt>’
               fvc::div(phi)==0
                              ^
make: *** [Make/linux64GccDPInt64Opt/inductionAC2_Foam.o] Error 1

the error is actually endless, so I hope it could be recognized from it's end, above.

I whould really appreciate any help. I already search the internet, saw the conservation term fvc::div(phi,u), but that is not the case for me. Couldn't find any explanation to fix my code.
thanks,

Or Werner.

[wernero]

Dear Bruno, first thank you for your comment.

In OpenFOAM programmer's guide (page 35) its written that divergence is an explicit operation. Is there a problem with the way I define my equation (or not an equation, as you say)? In that case, maybe correct the syntex could solve the problem?

I'm not familliar with conservation law using laplacian instead of divergence, but I will search for that. The conservation for current density is not my invention but very common, though. I will give it a thought.

thanks again,
Or.

[blais.bruno]

If it is an explicit operator, then you cannot solve an equation with it. You really need an FVM:: method to create a matrix system.

More often than not, when a quantity is conserved, there is a way to generate a laplacian equation to enforce the conservation of this quantity. The member of this laplacian equation becomes a lagrangian multiplier to enforce conservation.

For instance, for velocity, we have div(u) = 0
but we end up solving laplacian (p) = div (u) (minus or +, I never recall). and p is nothing more than a lagrangian multiplier that corrects the velocity so that it is divergence free.
Since the Maxwell equation share a lot of similarity with the Navier-Stokes equation, it is most likely that a similar process exists.

Cheers!
BB

Originally Posted by wernero

Dear Bruno, first thank you for your comment.

In OpenFOAM programmer's guide (page 35) its written that divergence is an explicit operation. Is there a problem with the way I define my equation (or not an equation, as you say)? In that case, maybe correct the syntex could solve the problem?

I'm not familliar with conservation law using laplacian instead of divergence, but I will search for that. The conservation for current density is not my invention but very common, though. I will give it a thought.

thanks again,
Or.

[wernero]

thank you again, Bruno. I still haven't got to the bottom of what you wrote, but it gave me some clue and I will read about that.

Meanwhile, about the first thing you wrote..
It is now understood that I can't solve with fvc only, it wasn't clear to me before but it makes sense.
Actually, the equation for Jext is only a secondary, and it should balance my source term (Jext is the source for the first equation). The physical source is in fact the current, which is of course the integral of J over the surface. It is quite difficult to divide homogeneous current into the real Jext, but homogeneous current density as a source is not exact. I'm telling that because I thought I could solve the first equation, and than let it run more iterations until both equation fulfill. As I understand what you said earlier, this cannot be done this way at all.

You gave me food for thought, I will have to go study some more.
Or