CFD Online Discussion Forums - How do I execute solve w/in a loop over which the solved variable is changing?

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - How do I execute solve w/in a loop over which the solved variable is changing? (https://www.cfd-online.com/Forums/openfoam-programming-development/127357-how-do-i-execute-solve-w-loop-over-solved-variable-changing.html)

How do I execute solve w/in a loop over which the solved variable is changing?

I thought this would be straight forward but it is tripping me up. Any help would be appreciated.

I need to iterate a 2D (geographical space, x and y) equation for each point in wave space (say kx and ky). I thought the following would work:

Code:

for (int m=0; m<nkx; ++m)

{

   for (int n=0; n<nky; ++n)

   {

      forAll(Umn,i)

      {

         Amn[i]      = wasd[m][n][i];

         Umn[i].x() = U[i].x() + some f(m,n);

         Umn[i].y() = U[i].y() + some f(m,n);

         Umn[i].z() = 0.0;

         Smn[i]      = s[m][n][i]

      }



      phi = (fvc::interpolate(Umn) & mesh.Sf());



      solve

      (

         fvm::ddt(Amn)+fvm::div(phi,Amn) == Smn

      )



      forall(Amn,i)

      {

         wasd[m][n][i] = Amn[i];

      }

   }

}

But it seems as though the solve function is only seeing Amn as it was for the first call to solve (although I have confirmed that Amn is being changed for each m and n.)

A simplified version of this problem is

Code:

#include "fvCFD.H"

#include "fvIOoptionList.H"

#include "simpleControl.H"



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



int main(int argc, char *argv[])

{

    #include "setRootCase.H"

    #include "createTime.H"

    #include "createMesh.H"

    #include "createFields.H"

    #include "createFvOptions.H"



    simpleControl simple(mesh);



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



    Info<< "\nCalculating scalar transport\n" << endl;



    #include "CourantNo.H"



    int l;



    scalar DD[64];



    for (l=0; l<64; ++l)

      {

    DD[l] = scalar(l+7);

      }



    for (l=0; l<64; ++l)

      {

    forAll(U,i)

      {

        T[i] = DD[l];

      }



    Info<<"T Before "<<T[56]<<endl;



    solve

      (

       fvm::ddt(T)

       );

    

    Info<<"T After "<<T[56]<<endl;

      }





    return 0;

}

which I adapted from scalarTransportFoam. No matter what T is before the solve (a uniform field with values varying from l=0 to 63), T after the solve is uniformly 7 (which is associated with l=0). Anyone have any ideas on how I am misinterpreting this and how to get around it? Thanks.

(Note that this post is an update to a post I made on 11/19)

Tom

Hi Tom,

What you are running into is a been of underlying OpenFoam thing.

When you are calling

Code:

fvm::ddt(T)

the building of the matrix obviously requires the values of T at the previous time step for building the matrix system. For a given field, you can only update the old time values once per time step (in an automated fashion, there might be manual ways), i.e. as soon as you have constructed the matrix one time, you have to step forward in time in order to be able to reassign the values in T from the previous time step.

Another detail, which affect you, is that if the old time values do not exist, then they are set equal to T itself, so as you are initialising T with 0 in the solving, you are also initialising the previous time step to 0. And it stays zero for the remaining solutions.

You can find the details on old time fields in:

GeometricField.C

Good luck,

Niels

Niels,

Thanks for your post. It was very helpful. Looks like I have some issues to resolve.

Tom