Issues with UnsteadyNSTurb and 21unsteadyNSTurb_RBF tutorial: nut field update and RBF interpolation

[carpemonf]

Hi @giovastabile and contributors,

I've been experimenting with UnsteadyNSTurb and encountered a few issues. While I've managed to resolve some, I wanted to consult with you before diving deeper into the matter.

The 21unsteadyNSTurb_RBF tutorial serves as the reproducible example.

1. The FOM nut field doesn't change over time
I've added the following line for updating nut in L166 ITHACA_FOMPROBLEMS/UnsteadyNSTurb/UnsteadyNSTurb.C:

+ nut = turbulence->nut());

2. Error when reading the parameters matrix
Reading the parameters from mu_samples_mat.txt resulted in missing the last line of the file, leading to a vector size mismatch error.

    auto mu_mat =
        ITHACAstream::readMatrix("./ITHACAoutput/Offline/mu_samples_mat.txt");
    example.velRBF = mu_mat.col(1);

A quick fix commenting out the condition on L139 src/ITHACA_CORE/ITHACAstream/ITHACAstream.C fixed the issue for this test.

-            if (i != (matrix.rows() - 1))
-            {
                ofs << std::endl;
-            }

3. RBF velocity coefficients interpolation

I got some results with the above changes but they were very bad (FOM up, ROM down).

Are the ReducedUnsteadyNSTurb and UnsteadyNSTurb classes designed to solve the eddy viscosity vector as a function of the velocity coefficients (a and a_dot) as described in Data-driven POD-Galerkin reduced order model for turbulent flows?

It looks to me like an issue or bug around the RBF implementation, which seems to apply velocity coefficients online interpolation on the top of weights obtained from the velocity parameters. Lifting was used for that tutorial, and interChoice is 1 by default. So, the tv variable will end up being the scaled velocity with values of e.g. 3.068395784 for a velocity inlet of 1.05, while velRBF is the inlet velocity vector with values of 1.0 and 1.1 for this case. This results in wrong values of newtonObjectSUP.gNut(i) in my view.

ReducedUnsteadyNSTurb.C:

        Eigen::VectorXd tv;
        Eigen::VectorXd aDer;
        aDer = (y.head(Nphi_u) - newtonObjectSUP.y_old.head(Nphi_u)) / dt;
        tv.resize(dimA);

        switch (interChoice)
        {
            case 1:
                tv << y.segment(firstRBFInd, dimA);
                break;

            case 2:
                tv << muStar, y.segment(firstRBFInd, dimA - muStar.size());
                break;

            case 3:
                tv << y.segment(firstRBFInd, dimA / 2), aDer.segment(firstRBFInd, dimA / 2);
                break;

            case 4:
                tv << muStar, y.segment(firstRBFInd, (dimA - muStar.size()) / 2),
                aDer.segment(firstRBFInd, (dimA - muStar.size()) / 2);
                break;

            default:
                tv << y.segment(firstRBFInd, dimA);
                break;
        }

        for (int i = 0; i < nphiNut; i++)
        {
            newtonObjectSUP.gNut(i) = problem->rbfSplines[i]->eval(tv);
        }

UnsteadyNSTurb.C:

                samples[i] = new SPLINTER::DataTable(1, 1);

                for (label j = 0; j < coeffL2.cols(); j++)
                {
                    samples[i]->addSample(velRBF.row(j), coeffL2(i, j));
                }

                rbfSplines[i] = new SPLINTER::RBFSpline(*samples[i],
                                                        SPLINTER::RadialBasisFunctionType::GAUSSIAN, false, radii(i));

I tried the RBF time-parameters implementation and had relatively good results (without extrapolation), but wanted to get this right.

Is there any suggestion or insights you might have?

Thanks!

[carpemonf]

For future reference, I looked further into the issue and got good results by addressing:

• FOM. There are already functions implemented such as velDerivativeCoeff, and others, to compute the matrices needed for the RBF interpolation with velocity coefficients. These functions are not used within the class so I employed them to get the values to "train" the RBF.
• ROM. In the online stage, gNutAve calculations need to be added on solveOnlineSUPAve and solveOnlinePPEAve.

[giovastabile]

Hi @carpemonf. I am sorry for the late response. From my understanding, you managed to solve the issue right? If this is the case could you please do a PR?

[carpemonf]

@giovastabile, the issue has been partially resolved but additional logic and refactoring are needed for a PR. Further work is still needed to make the different combinations of DAEs and RBF interpolation strategies functional.

I also wanted to verify that these changes don't disrupt other tutorials that use this class, such the offline stage in 22datadriven_corrections. I tried to run that tutorial and reproduce the results reported in the Jupyter notebooks but found a number of issues. The offline run breaks at the RBF part of the project function - it expects data to build the mapping for g. I added rbfInterp to false to skip this step.

I found some other issues, such as missing ROM matrices *.npy files loaded in the notebook, and other minor issues. The offline app and dictionary are configured for 50 modes. I understand that for the comparison with the projection
of the full order solution we need coeffs for both 50 and 5 modes together with the matrices for 5 modes, right?

I'm willing to work on these issues and can submit a PR once they're resolved if that would be helpful.

[carpemonf]

Sorry, another note for the 22datadriven_corrections tutorial. Running it with a high number of modes (50) results in a division by zero error at modesEig.col(i) = modesEig.col(i).array() / normFact(i, 0); in ITHACAPOD.C for nut. Shall we set the modes to zero for zero eigenvalues, suggesting that the modes do not capture any dynamics of the system?

Thank you.

[carpemonf]

Never mind my previous message, this was originated by calculating the modes for a uniform field (nut was constant and not updated over time).

I also noticed that the FOM results were not good. It took awhile to figure out but a copy of the U field was used for the PIMPLE solution and then values are not updated back. As a result, the calculation of the turbulence model always uses the initial U and not the calculated one. I fixed that in PR #549 for your consideration.

[carpemonf]

The 22datadriven_corrections tutorial worked for me with the changes in PR #549 and other fixes. I will create another PR and come back to the initial issue with 21unsteadyNSTurb_RBF.

Thanks.