[chyczewski]

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:

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

#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

[ngj]

Hi Tom,

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

When you are calling

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

[chyczewski]

Niels,

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

Tom