simulationsetup2d.hh

来自「open lattice boltzmann project www.open」· HH 代码 · 共 488 行 · 第 1/2 页

                    (lambda/(T)4) * (sumPressure + poissonTerm.get(iX,iY) );
            }
        }

        maxResidue = std::numeric_limits<T>::min();
        for (int iX=1; iX<lx-1; ++iX) {
            for (int iY=1; iY<ly-1; ++iY) {
                T sumPressure =
                      pressure.get(iX+1,iY) +
                      pressure.get(iX,iY+1) +
                      pressure.get(iX-1,iY) +
                      pressure.get(iX,iY-1);
                T residue = fabs(sumPressure -(T)4*pressure.get(iX,iY)
                                            + poissonTerm.get(iX,iY));
                if (residue > maxResidue) maxResidue = residue;
            }
        }

        if (iter%20==0) {
            std::cout << "SOR iteration " << iter
                      << ": max residue= "
                      << maxResidue/averagePoisson << std::endl;
        }
        ++iter;
    }
    while (maxResidue/averagePoisson>epsilon);

    iter = 0;
    T difference = (T)1;
    do {
        for (int iX=1; iX<lx-1; ++iX) {
            pressure2.get(iX,0) = pressure2.get(iX,1);
            pressure2.get(iX,ly-1) = pressure2.get(iX,ly-2);
        }
        for (int iY=1; iY<ly-1; ++iY) {
            pressure2.get(0,iY) = pressure2.get(1,iY);
            pressure2.get(lx-1,iY) = pressure2.get(lx-2,iY);
        }
        pressure2.get(0,0)       = pressure2.get(1,1);
        pressure2.get(lx-1,0   ) = pressure2.get(lx-2,2);
        pressure2.get(0,ly-1)    = pressure2.get(1,ly-2);
        pressure2.get(lx-1,ly-1) = pressure2.get(lx-2,ly-2);

        for (int iX=1; iX<lx-1; ++iX) {
            for (int iY=1; iY<ly-1; ++iY) {

                T sumPressure =
                      pressure2.get(iX+1,iY) +
                      pressure2.get(iX,iY+1) +
                      pressure2.get(iX-1,iY) +
                      pressure2.get(iX,iY-1);

                pressure2.get(iX,iY) = 
                    ((T)1-lambda) * pressure2.get(iX,iY) +
                    (lambda/(T)4) * (sumPressure + poissonTerm.get(iX,iY) );
            }
        }

        difference = T();
        for (int iX=0; iX<lx; ++iX) {
            for (int iY=0; iY<ly; ++iY) {
                difference
                    += util::sqr(pressure.get(iX,iY)-pressure2.get(iX,iY));
            }
        }
        difference = sqrt(difference/(T)(lx*ly))/4e-4;

        if (iter%20==0) {
            std::cout << difference << std::endl;
        }
        ++iter;
    }
    while (difference>epsilon);

    for (int iX=0; iX<lx; ++iX) {
        for (int iY=0; iY<ly; ++iY) {
            lattice.get(iX,iY).defineRho (
                    pressure.get(iX,iY)*Lattice<T>::invCs2 + (T)1 );
        }
    }
}


template <typename T, template<typename U> class Lattice>
void convergeFixedVelocity(BlockLattice2D<T,Lattice>& lattice,
                           T epsilon, int step)
{
    BlockLatticeView2D<T,Lattice> latticeView(lattice);
    TensorFieldBase2D<T,2> const& velocity = latticeView.getDataAnalysis().getVelocity();

    T maxF = T();
    for (int iX=0; iX<latticeView.getNx(); iX+=step) {
        for (int iY=0; iY<latticeView.getNy(); iY+=step) {
            for (int iF=0; iF<Lattice<T>::q; ++iF) {
                T f = latticeView.get(iX,iY)[iF];
                if (f>maxF) {
                    maxF = f;
                }
            }
        }
    }

    std::vector<T> oldValues, currentValues;
    T diff=(T)1;
    do {
        lattice.staticCollide(velocity);
        lattice.stream(true);

        for (int iX=0; iX<latticeView.getNx(); iX+=step) {
            for (int iY=0; iY<latticeView.getNy(); iY+=step) {
                currentValues.push_back (
                    latticeView.get(iX,iY)[(iX+iY) % (Lattice<T>::q)] );
            }
        }

        if (oldValues.size() == currentValues.size()) {
            diff = T();
            for (unsigned iVal=0; iVal<oldValues.size(); ++iVal) {
                T candidate = std::abs(currentValues[iVal]-oldValues[iVal]) 
                                / maxF;
                if (candidate>diff) {
                    diff=candidate;
                }
            }
        }
        oldValues.swap(currentValues);
        currentValues.clear();
    }
    while (diff > epsilon);
}

template <typename T, template<typename U> class Lattice>
void testLiShi(BlockLattice2D<T,Lattice>& lattice, T epsilon, T lambda)
{
    BlockLatticeView2D<T,Lattice> latticeView(lattice);
    DataAnalysisBase2D<T,Lattice> const& analysis = latticeView.getDataAnalysis();
    TensorFieldBase2D<T,2> const& velocity = analysis.getVelocity();

    int lx = lattice.getNx();
    int ly = lattice.getNy();

    ScalarField2D<T> const& poissonTerm = analysis.getPoissonTerm();

    ScalarField2D<T> pressure(lx, ly); pressure.construct();

    T averagePoisson = T();
    for (int iX=0; iX<lx; ++iX) {
        for (int iY=0; iY<ly; ++iY) {
            pressure.get(iX,iY) = T();
            averagePoisson += util::sqr(poissonTerm.get(iX,iY));
        }
    }
    averagePoisson /= (lx*ly);
    averagePoisson = sqrt(averagePoisson);

    int iter=0;
    T maxResidue = (T)1;
    do {
        for (int iX=1; iX<lx-1; ++iX) {
            pressure.get(iX,0) = pressure.get(iX,1);
            pressure.get(iX,ly-1) = pressure.get(iX,ly-2);
        }
        for (int iY=1; iY<ly-1; ++iY) {
            pressure.get(0,iY) = pressure.get(1,iY);
            pressure.get(lx-1,iY) = pressure.get(lx-2,iY);
        }
        pressure.get(0,0)       = pressure.get(1,1);
        pressure.get(lx-1,0   ) = pressure.get(lx-2,2);
        pressure.get(0,ly-1)    = pressure.get(1,ly-2);
        pressure.get(lx-1,ly-1) = pressure.get(lx-2,ly-2);

        for (int iX=1; iX<lx-1; ++iX) {
            for (int iY=1; iY<ly-1; ++iY) {

                T sumPressure =
                      pressure.get(iX+1,iY) +
                      pressure.get(iX,iY+1) +
                      pressure.get(iX-1,iY) +
                      pressure.get(iX,iY-1);

                pressure.get(iX,iY) = 
                    ((T)1-lambda) * pressure.get(iX,iY) +
                    (lambda/(T)4) * (sumPressure + poissonTerm.get(iX,iY) );
            }
        }

        maxResidue = std::numeric_limits<T>::min();
        for (int iX=1; iX<lx-1; ++iX) {
            for (int iY=1; iY<ly-1; ++iY) {
                T sumPressure =
                      pressure.get(iX+1,iY) +
                      pressure.get(iX,iY+1) +
                      pressure.get(iX-1,iY) +
                      pressure.get(iX,iY-1);
                T residue = fabs(sumPressure -(T)4*pressure.get(iX,iY)
                                            + poissonTerm.get(iX,iY));
                if (residue > maxResidue) maxResidue = residue;
            }
        }

        if (iter%20==0) {
            std::cout << "SOR iteration " << iter
                      << ": max residue= "
                      << maxResidue/averagePoisson << std::endl;
        }
        ++iter;
    }
    while (maxResidue/averagePoisson>epsilon);

    T diff=(T)1;
    for(int i=0; i<50000; ++i) {
        lattice.staticCollide(velocity);
        lattice.stream();

        T averageRho = T();
        for (int iX=0; iX<lx; ++iX) {
            for (int iY=0; iY<ly; ++iY) {
                averageRho += lattice.get(iX,iY).computeRho();
            }
        }
        
        averageRho /= (T)(lx*ly);

        T difference = T();
        for (int iX=0; iX<lx; ++iX) {
            for (int iY=0; iY<ly; ++iY) {
                difference += util::sqr(
                        pressure.get(iX,iY)-
                        (lattice.get(iX,iY).computeRho()-averageRho)/3.);
            }
        }
        difference = sqrt(difference/(T)(lx*ly))/4e-4;

        if (i%20==0) {
            std::cout << difference << std::endl;
        }
    }
    while (diff > epsilon);
}

}  // namespace olb

#endif