[frobaux]

Hello everyone,
I've a noticed a strange behavior in my code, I'm sure it is my mistake or my lack of understanding but I could understand it:


I have a field that I want to interpolate. I've coded something like

autoPtr<interpolation<vector>> interpV
(
       interpolation<vector>::New(mapMethod2, srcV)
 );
forAll(tgtMesh.C(), tgtCelli)
 {
      const point & p=tgtMesh.C()[tgtCelli];
      const label & srcCelli= srcMesh.findCell(p);

       VfromExternal[tgtCelli]= interpV().interpolate(p, srcCelli, -1);
}

where srcV is my source field and mapMethod2 is "cellPoint". By diving into the source, and written a couple of "Info <<" inside interpolation.H/C and interpolationCellPoint.H/C it seems that the "interpolation::new" is called. Inside this function:

    return autoPtr<interpolation<Type>>(cstrIter()(psi));

Which in turns call interpolationCellPoint constructor:

template<class Type>
Foam::interpolationCellPoint<Type>::interpolationCellPoint
(
    const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
    interpolation<Type>(psi),
    psip_
    (
        volPointInterpolation::New(psi.mesh()).interpolate
        (
            psi,
            "volPointInterpolate(" + psi.name() + ')',
            true        // use cache
        )
    )
    
{
    psip_.write(); // local modif
    Info << "should write pointField" << psip_.name() << endl; //local modif
    // Uses cellPointWeight to do interpolation which needs tet decomposition
    (void)psi.mesh().tetBasePtIs();
}

But psip_ is not computed, except at the first time step! (but the Info in blue prints out, so the contructor is called, but volPointInterpolation is definetly not called => the psip_ written in blue above stays the same during to consecutive time step, except in the openFoam header, even with psi different(checked)!!). The really troubling thing is that the same procedure applied to the pressure works well and its "psip_" seems to be recomputed at every time step.


My guess is the new srcV is so close (~10^-3) from the previous one that autoPtr does not actually recompute the volPointInterpolation needed : psip_ is not recomputed inside interpolationCellPoint.H. The problem is I need to and the final results are really sensitive to that! should I manually destroy the interpV at the end? (and how do I do that?) Am I doing something wrong?


Edit: The same behavior is obtained when I tried to define interpV as a const reference of directly interpolationCellPoint

const interpolationCellPoint<vector>  & interpV (srcV);

=> even if srcV changes (slightly) from 1 instant to the next, interpV.psip_ doesn't!!!



Thanks a lot for your inputs!
Sincerely,

Fabien

[frobaux]

Well, after a lot of print and research everywhere, I've found the culprit!!!


The problem was in

psip_
    (
        volPointInterpolation::New(psi.mesh()).interpolate
        (
            psi,
            "volPointInterpolate(" + psi.name() + ')',
            true        // use cache // THIS IS THE CULPRIT
        )
    )

Where the "cache" was used. Digging deeper I saw that my modification of the sourceV didn't modify sourceV.eventNo_.

Now, I do it manually

sourceV.evenNo()=sourceV.eventNo()+1 // this is wrong, see edits

and it works fine


EDIT:

I was wrong, it does not run well. I forgot to recompile the "true".. Inside the volPointInterpolation::New(psi.mesh()).interpolate the eventNo is way to important. +1 does not help matching it...


EDIT2:
Well I was an idiot, I found out that -If this time I'm right- the eventNo is more of a global number counting all occurence. To set a n object as updated, used myObject.setUpToDate() instead (it sets myObject.eventNo_= db().getEvent() )