上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答3.2.2

adf322Foam.C

/*---------------------------------------------------------------------------*\
	Changed from scalarTransportFoam to adv121Foam

Application
    adf322Foam

Description
    Solves training examples 3.2.2 problem.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"

scalar get_xhat(scalar x0, scalar y0, scalar t){
    using std::sin;
    using std::cos;
    scalar xhat = x0*cos(t) - y0*sin(t);
    return xhat;
}

scalar get_yhat(scalar x0, scalar y0, scalar t){
    using std::sin;
    using std::cos;
    scalar yhat = y0*cos(t) - x0*sin(t);
    return yhat;
}

scalar get_my_r2(scalar x, scalar xhat, scalar y, scalar yhat){
    scalar my_r2 = 0;
    my_r2 = (x-xhat)*(x-xhat) + (y-yhat)*(y-yhat);
    return my_r2;
}

scalar get_phiExt(scalar r2, scalar t){
    scalar epsilon = 1e-3;
    scalar phiExt = 0;
    using Foam::exp;
    phiExt = (1/(4*M_PI*epsilon*t)*exp(-(r2/(4*epsilon*t))));
    return phiExt;
}

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

int main(int argc, char *argv[])
{
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    simpleControl simple(mesh);

    #include "createFields.H"

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

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

    #include "CourantNo.H"

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            fvScalarMatrix TEqn
            (
                fvm::ddt(T)
              + fvm::div(phi, T)
              - fvm::laplacian(DT, T)
             ==
                fvOptions(T)
            );

            TEqn.relax();
            fvOptions.constrain(TEqn);
            TEqn.solve();
            fvOptions.correct(T);
        }

        runTime.write();
    }
    
    // t = runTime.value()
    // scalar my_t = 1.0;
    // scalar my_a = 1.0;
    
    volScalarField T_ex(T);
    using std::sqrt;
    
    // Calculate Exact Number
    
    forAll(T_ex,celli)
    {
        scalar xx = mesh.C()[celli].x();
        scalar yy = mesh.C()[celli].y();
        scalar my_r21 = get_my_r2(xx, xhat1, yy, yhat1);
        scalar phiExt1 = get_phiExt(my_r21, my_t1);
        T_ex[celli] = phiExt1;
    }
    
    // Error Analysis
    
    scalar L1=0;
    scalar up = 0;
    scalar low =0;
    
    forAll(T,celli)
    {
        up += mag(T_ex[celli]-T[celli]) * mesh.V()[celli];
        low += mag(T_ex[celli])*mesh.V()[celli];
    }
    L1 = up/low;
    
    Info << "L1 error = " << L1 << endl;
    
    scalar L2 = 0;
    up = 0;
    low = 0;
    forAll(T,celli)
    {
    	up += mag(T_ex[celli]-T[celli])*mag(T_ex[celli]-T[celli])*mesh.V()[celli];
    	low += T_ex[celli]*T_ex[celli]*mesh.V()[celli];
    }
    using std::sqrt;
    L2 = sqrt(up/low);
    
    Info << "L2 error = " << L2 << endl;
    
    scalar L3 = 0;
    scalar up_max = 0;
    scalar low_max = 0;
    
    forAll(T,c