ExtractCriticalLayerFunctions.h File Reference

#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>

Go to the source code of this file.

Functions
void	Computestreakpositions (MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans=NekConstants::kNekUnsetDouble)

Function Documentation

Computestreakpositions()

void Computestreakpositions	(	MultiRegions::ExpListSharedPtr &	streak,
		Array< OneD, NekDouble > &	xc,
		Array< OneD, NekDouble > &	yc,
		NekDouble	cr,
		NekDouble	trans = NekConstants::kNekUnsetDouble
	)

Definition at line 115 of file ExtractCriticalLayer.cpp.

{
    int i;
    int npts = xc.size();
 
    int nq = streak->GetTotPoints();
    Array<OneD, NekDouble> derstreak(nq);
    Array<OneD, NekDouble> x(nq);
    Array<OneD, NekDouble> y(nq);
    streak->GetCoords(x, y);
 
    streak->BwdTrans(streak->GetCoeffs(), streak->UpdatePhys());
    streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
 
    // set intiial xc to be equispaced over mesh and yc to be zero
    NekDouble x_max = Vmath::Vmax(nq, x, 1);
    NekDouble x_min = Vmath::Vmin(nq, x, 1);
 
    for (i = 0; i < npts; ++i)
    {
        xc[i] = x_min + (x_max - x_min) * i / ((NekDouble)(npts - 1));
        yc[i] = 0.0;
    }
 
    int elmtid, offset, cnt;
    NekDouble U, dU;
    NekDouble F;
    NekDouble ConvTol  = 1e-9;
    NekDouble CoordTol = 1e-5;
    int maxiter        = 100;
    Array<OneD, NekDouble> coord(2);
    NekDouble ytmp;
    // Do Newton iteration on y direction
    cerr << "[";
    for (int e = 0; e < npts; e++)
    {
        coord[0] = xc[e];
        coord[1] = yc[e];
        cout << e << endl;
        if (!(e % 10))
        {
            cerr << ".";
        }
 
        F           = 1000;
        cnt         = 0;
        int its     = 0;
        int attempt = 0;
        // cout<<"start guess:  x="<<xc[e]<<"    y="<<yc[e]<<endl;
        while (abs(F) > 0.000000001)
        {
            ytmp   = coord[1];
            elmtid = streak->GetExpIndex(coord, 0.00001);
            offset = streak->GetPhys_Offset(elmtid);
            U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +
                                                                offset);
            dU =
                streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
            coord[1] = coord[1] - (U - cr) / dU;
            F        = U - cr;
            ASSERTL0(coord[0] == xc[e], " x coordinate must remain the same");
            // stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord,
            // streak->GetPhys() +offset );
            // cout<<"elmtid="<<elmtid<<"  x="<<coord[0]<<"   y="<<coord[1]<<"
            // stvalue="<<U<<"   dU="<<dU<<endl;
            if (abs(coord[1]) > 1)
            {
                coord[1]       = ytmp + 0.01;
                elmtid         = streak->GetExpIndex(coord, 0.00001);
                offset         = streak->GetPhys_Offset(elmtid);
                NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(
                    coord, streak->GetPhys() + offset);
                NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(
                    coord, derstreak + offset);
                coord[1] = coord[1] - (Utmp - cr) / dUtmp;
 
                if ((abs(Utmp - cr) > abs(F)) || (abs(coord[1]) > 1))
                {
                    coord[1] = ytmp - 0.01;
                }
 
                attempt++;
            }
            else
            {
                ASSERTL0(abs(coord[1]) <= 1, " y value out of bound +/-1");
            }
 
            its++;
            if (its > 1000 && abs(F) < 0.0001)
            {
                cout << "warning streak position obtained with precision:" << F
                     << endl;
                break;
            }
            else if (its > 1000)
            {
                ASSERTL0(false, "no convergence after 1000 iterations");
            }
        }
 
        yc[e] = coord[1];
    }
    cerr << "]" << endl;
 
    // output to interface file
    FILE *fp = fopen("interfacedat.geo", "w");
 
    NekDouble y_max = Vmath::Vmax(nq, y, 1);
    NekDouble y_min = Vmath::Vmin(nq, y, 1);
 
    cnt = 1;
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n", cnt++, x_min, y_min);
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n", cnt++, x_max, y_min);
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n", cnt++, x_max, y_max);
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n", cnt++, x_min, y_max);
 
    for (i = 0; i < npts; ++i)
    {
        fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++, xc[i],
                yc[i]);
    }
 
    fclose(fp);
 
    // output to interface_up file as bend of vertical shift and 45 degrees
    // shift
    fp = fopen("interfacedat_up.geo", "w");
 
    NekDouble nx, ny, norm;
 
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++, xc[0],
            yc[0] + trans);
 
    for (i = 1; i < npts - 1; ++i)
    {
        norm = sqrt((xc[i + 1] - xc[i - 1]) * (xc[i + 1] - xc[i - 1]) +
                    (yc[i + 1] - yc[i - 1]) * (yc[i + 1] - yc[i - 1]));
        nx   = (yc[i - 1] - yc[i + 1]) / norm;
        ny   = (xc[i + 1] - xc[i - 1]) / norm;
 
        fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++,
                xc[i] + nx * trans, yc[i] + ny * trans);
    }
 
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++, xc[npts - 1],
            yc[npts - 1] + trans);
 
    // output to interface_up file as bend of vertical shift and 45 degrees
    // shift
    fp = fopen("interfacedat_dn.geo", "w");
 
    trans = -trans;
 
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++, xc[0],
            yc[0] + trans);
 
    for (i = 1; i < npts - 1; ++i)
    {
        norm = sqrt((xc[i + 1] - xc[i - 1]) * (xc[i + 1] - xc[i - 1]) +
                    (yc[i + 1] - yc[i - 1]) * (yc[i + 1] - yc[i - 1]));
        nx   = (yc[i - 1] - yc[i + 1]) / norm;
        ny   = (xc[i + 1] - xc[i - 1]) / norm;
 
        fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++,
                xc[i] + nx * trans, yc[i] + ny * trans);
    }
 
    fprintf(fp, "Point(%d)={%12.10lf,%12.10lf,0,1.0};  \n", cnt++, xc[npts - 1],
            yc[npts - 1] + trans);
}

References tinysimd::abs(), ASSERTL0, Nektar::MultiRegions::eY, tinysimd::sqrt(), Vmath::Vmax(), and Vmath::Vmin().

Referenced by GetStreakLocation(), and main().