[praphul]

Hi,


I was trying to modify the V momentum equation with a source term. I tried it with the vectorCodedSource term specification as given in the below code snippet. The issue is that the case runs but the source terms are not being included while solving the problem. i.e the solution I get is similar to the solution I get without employing the source term.

I am using simpleFoam (openFOAM v1906) as my solver and the test case I am solving corresponds to a lid driven cavity with manufactured solutions. The reference paper is :
Shih, T. M., Tan, C. H., & Hwang, B. C. (1989). Effects of grid staggering on numerical schemes. International Journal for Numerical Methods in Fluids, 9(2), 193–212. https://doi.org/10.1002/fld.1650090206

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system" ;
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


MomentumSource 

{ type vectorCodedSource ;
  name MomentumSource ;
  active on ;  
  fields   (U) ;
  selectionMode all ;
  redirectType velocitySource; 
   
  codeInclude
  #{

  #};

  codeCorrect
  #{

    Info<< "CodeCorrect" << endl << "\n" ; 
  
   #};

  codeConstrain
  #{
   
    Info<< "CodeConstrain" << endl << "\n" ; 

  #};


  codeAddSup
  #{ 
     
     Info<< "codeAddSup" << endl << "\n" ; 

     const scalarField&    V  = mesh_.V() ;  // Cell volume  

     vectorField& VSource = eqn.source() ;   // Equation source definition

     const vectorField& cell  = mesh_.C()  ;
     const scalarField& cellx = cell.component(0) ; 
     const scalarField& celly = cell.component(1) ; 
     

     int n = cell.size() ;  
     double arg[n] = {} ; 
     double v1[n] = {} ;   
     double v2[n] = {} ;  
     double v3[n] = {} ;  
     double v4[n] = {} ;  
     double v5[n] = {} ;  
     double v6[n] = {} ;  
     double v7[n] = {} ;  
     double v8[n] = {} ;  
     double v9[n] = {} ;  
     double v10[n] = {} ;  

Info<<"VSource"<<VSource<<"\n" ; 

     forAll(cell, i) 
     {

     v1[i] = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4) - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;


     v2[i] = -768*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 8*pow(cellx[i],6) - 10*pow(cellx[i],5) + 8*pow(cellx[i],4) - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],5) ; 


     v3[i] = 38.4*pow(cellx[i],5) ;  

     v4[i] = -96*pow(cellx[i],4) ; 

     v5[i] = 256*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 8*pow(cellx[i],6) - 10*pow(cellx[i],5) + 8*pow(cellx[i],4) - 4*pow(cellx[i],3) +  pow(cellx[i],2) ) * pow(celly[i],3) ;

     v6[i] = -64*pow(cellx[i],3) ; 


     v7[i] =  96*(8*pow(cellx[i],3) - 12*pow(cellx[i],2) + 2*pow(cellx[i],1) + 1)*pow(celly[i],2) ;

     v8[i] = 192*pow(cellx[i],2) ;

     v9[i] = -128*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 6*pow(cellx[i],6) - 4*pow(cellx[i],5) + pow(cellx[i],4))*pow(celly[i],1);

     v10[i] = -64*cellx[i] ; 
 
      arg[i] = (v1[i] + v2[i] + v3[i] + v4[i] + v5[i] + v6[i] + v7[i] + v8[i] + v9[i] + v10[i])*V[i]  ; 

      VSource[i] = vector(0,arg[i],0) ;

     } 



  #} ; 


code
        #{
            $codeInclude
            $codeCorrect
            $codeConstrain
            $codeAddSup
           $codeSetValue
        #};
  }

Any ideas on tackling this issue ?

Thanks in advance,
Praphul

[praphul]

I tried the same with the following command. Just to see whether it invokes the source term but it doesnt do so. Am I doing anything wrong?

 vectorField& Usource = eqn.source();

       Usource += vector(0, 10000, 0);

[pbachant]

You may want to try creating a field an adding it to eqn directly, i.e., eqn += myField, or you could try adding to the equation's source rather than setting the values explicitly. I don't know enough about the inner workings to provide more insight than that!

[praphul]

Thank you pbachant. I tried setting the field directly to the equation as you had mentioned. Unfortunately, It does not work. But, I Thank you for your time and consideration.

[ELwardi]

Well, I see you tagged the thread with openfoamv2019; I don't have OF installed (Just fragments of foam-extend 4 actually ). So I'll try to help with no promise that the code will even compile:


1. Avoid using C ARRAYS. At least, use Foam::List if absolutely required .


2. You create way too much temp vars; Try getting to the point faster: You want to modify the second component in the source of UEqn? get to it immediately in the forAll loop.


3. Always ADD values to source, don't assign to it !!!


It has been a while since I last medeled with fvOptions, but I will share what I would write to accomplish this (Of course, I welcome any improvements/suggestions):

codedSource
{
    type            coded;
    selectionMode   all;

    fields          (U);
    name            MomentumSource;

    codeAddSup
    #{
        const scalarField& V = mesh_.V();
        scalarField& USource = eqn.source();
        const scalarField& cellx = mesh_.C().component(0) ; 
        const scalarField& celly = mesh_.C().component(0) ;
        
        forAll(USource, i)
        {
            scalar v1 = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4)
                             - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;
            // .... and the other 9 scalars
            
            // Please ADD to the source, don't overwrite it.
            // Overwriting it means you don't care about source terms in 
            // DESCRETISATION OPERATORS (laplacian, div ...)
            USource[i][1] += (v1+v2+...)*V[i];
        }
    #};
}

I hope it fixes your problem. Maybe you didn't notice the change by overwriting the source.

[praphul]

Elwardi, I tried ur suggestion and unfortunately that doesn't work too. It seems I am stuck with this problem for eternity.

1) I am curious to know on ur suggestion as to why I should not use C arrays and instead I should use list.



2) The Velocity field is supposed to be a vector quantity. Therefore, I am supposed to use vectorCodedSource isnt it ? I guess the way you have specified it might be because you use a different version of OpenFoam.


3) Why is it said that we should always add to the source rather than assigning it ?

I am attaching the zip file of my problem setup with this mail. May be if you or someone else can look into that, it can be helpful .

Regards,
Praphul

[ELwardi]

Ok, let me clarify things a bit:


1. It's not like C arrays are particularly bad for this situation, but it's good to jump out of the habit of using them (No bounds checking, possibility to overflow the stack ...) http://www.cplusplus.com/forum/articles/17108/


2. I meant "vectorField& USource = eqn.source();", that was a copy-paste typo.


3. I was considering that equation source may be already be populated by the time fvOptions kicks in!
Suppose your fvm::laplacian generates an explicit term (fvm::ddt certainly does in transient simulations). That term will be stored in the source vector. By the time your additional source is applied, if you overwrite what was there with your source's value, you're no longer solving the problem you want to solve
But it appears recent fvOptions generate their own hole "eqn" then adds it the main equation. I almost always access sources directly so i got confused


4. I had problems compiling your U boundary conditions (writeEntry seems to be dropped in OF7) but I managed to get it working (At least there is a change) with basic U BCs. Let me know if the problem persists for you

[praphul]

Thank you for clarifying my doubts Elwardi.



1) However, in the file 0/U you had attached I see that you had specified zeroGradient boundary condition for the walls which I think should be changed to noSlip condition.

2) I do not understand why the celly component was specified like this is in the fvOptions.

        const scalarField& celly = mesh_.C().component(0) ;

y component should have an index 1 isnt it ? i.e

        const scalarField& celly = mesh_.C().component(1) ;

I changed both these criterion and ran the case. I found that the solution I get is exactly the same as with fvOptions and without fvOptions. I dont understand why that is so. I have attached graphs for the same.

[ELwardi]

After fixing the issues you mentioned, I got slower convergence with the source enabled (340 iterations vs 330 iters in a x4 finer mesh). Maybe the source is too week to notice the change, try to multiply USource with 10 in all cells, there will be a clear change to notice. Revise your vi quantities, there might be a bug in there!

[praphul]

In the particular case of manufactured solution which I solve, I assume a solution for u, v and plug it in to the governing equations to obtain the source terms for the respective equation. In this case, I get the x momentum source term to be zero and the y momentum source term to be what I have specified in the fvOptions. Then I solve the equations again to get back the solution I have assumed. Its like a merry-go-round but it is used to verify the code.

To get the same solution as the assumed one, the boundary conditions needs to be in a format I had specified in the zip file earlier.

I dont think the source term comes out to be small as I had tried printing out the value of sources using Info command and it happen to be finite value. The research paper suggests that the source term at x = 0.5 and y = 0.5 and Reynolds number of 1 equals 3.356250 and I get the same value. In this example, this value equals source term divded by the volume of the cell. I am supposed to get back the original solution I had assumed unless there is some problem with the code.

I had cross checked the vi terms too. Infact I had missed a term i.e 96*(2x-1)*y^4 in the fvOption file I had specified, I corrected it but the solution remains same. Still wondering what the issue could be.

[ELwardi]

1. Double check your units! I think UEqn has [0 4 -2 0 0 0 0]



2. You can use:
debugSwitches{
vectorCodedSource 1;

}
in controlDict to verify when ::addSup method is called and whether it does what you want when you want.

[CFD_10]

Originally Posted by ELwardi

1. Double check your units! I think UEqn has [0 4 -2 0 0 0 0]



2. You can use:
debugSwitches{
vectorCodedSource 1;

}
in controlDict to verify when ::addSup method is called and whether it does what you want when you want.

I cannot see vectorCodedSource in the global controlDict of (OF5, OF6, and OF 1906).

grep -i 'codedVectorSource' $FOAM_SRC/controlDict

This will return nothing.

[ELwardi]

`simpleFoam -listSwitches | grep Source`
Tells me there is one in OF7; there should be something similar eg. codedSource, run the command and see if you can find any useful switch

[CFD_10]

Originally Posted by ELwardi

`simpleFoam -listSwitches | grep Source`
Tells me there is one in OF7; there should be something similar eg. codedSource, run the command and see if you can find any useful switch

According to the solver help, the option -listSwitches : List switches declared in libraries but not set in etc/controlDict