PABLO_example_00010.cpp

Parallel 3D dynamic adaptive mesh refinement (AMR) with data using PABLO.

Parallel 3D dynamic adaptive mesh refinement (AMR) with data using PABLOIn this example an Adaptive Mesh Refinement (AMR) around a moving sphere is performed.

To run: ./PABLO_example_00010
To see the result visit: PABLO website

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit is distributed in the hope that it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#if BITPIT_ENABLE_MPI==1
#include <mpi.h>
#endif
 
#include "bitpit_PABLO.hpp"
 
#if BITPIT_ENABLE_MPI==1
#include "PABLO_userDataComm.hpp"
#include "PABLO_userDataLB.hpp"
#endif
 
using namespace std;
using namespace bitpit;
 
// =================================================================================== //
// =================================================================================== //

void run()
{
    int iter = 0;

    PabloUniform pablo10(3);

    int idx = 0;
    pablo10.setBalance(idx,false);

    for (iter=1; iter<6; iter++){
        pablo10.adaptGlobalRefine();
    }
 
#if BITPIT_ENABLE_MPI==1
    pablo10.loadBalance();
#endif

    double xc, yc, zc;
    double radius = 0.15;
    double Dt = 0.000025;
    double omega = 2.0*3.14/0.001;

    xc = 0.25*cos(omega* Dt) + 0.5 ;
    yc = 0.25*sin(omega* Dt) + 0.5;
    zc = 100*Dt;

    uint32_t nocts = pablo10.getNumOctants();
    vector<double> oct_data(nocts, 0.0);

    int itstart = 1;
    int itend = 460;
 
    for (iter=itstart; iter<itend; iter++){
        xc = 0.25*cos(omega* Dt* (double) iter) + 0.5 ;
        yc = 0.25*sin(omega* Dt* (double) iter) + 0.5;
        zc = 100*Dt*iter;
 
        for (unsigned int i=0; i<nocts; i++){
            bool inside = false;
            vector<array<double,3> > nodes = pablo10.getNodes(i);
            array<double,3> center = pablo10.getCenter(i);
            oct_data[i] = sqrt((pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0)+pow((center[2]-zc),2.0)));
            for (int j=0; j<8; j++){
                double x = nodes[j][0];
                double y = nodes[j][1];
                double z = nodes[j][2];
                if (pow((x-xc),2.0)+pow((y-yc),2.0)+pow((z-zc),2.0) <= pow(radius,2.0)){
                    if (pablo10.getLevel(i) < 6){
                        pablo10.setMarker(i,1);
                    }
                    else{
                        pablo10.setMarker(i,0);
                    }
                    inside = true;
                }
            }
            if (pablo10.getLevel(i) > 5 && !inside){
                pablo10.setMarker(i,-1);
            }
        }

        pablo10.preadapt();
        pablo10.adapt();
 
#if BITPIT_ENABLE_MPI==1
        pablo10.loadBalance();
#endif

        nocts = pablo10.getNumOctants();
        vector<double> oct_data_new(nocts, 0.0);
        for (unsigned int i=0; i<nocts; i++){
            array<double,3> center = pablo10.getCenter(i);
            oct_data_new[i] = sqrt((pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0)+pow((center[2]-zc),2.0)));
        }
 
        oct_data.resize(nocts);
        oct_data = oct_data_new;
        oct_data_new.clear();

        pablo10.updateConnectivity();
        pablo10.writeTest("pablo00010_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);
    }
}

int main(int argc, char *argv[])
{
#if BITPIT_ENABLE_MPI==1
    MPI_Init(&argc,&argv);
#else
    BITPIT_UNUSED(argc);
    BITPIT_UNUSED(argv);
#endif
 
    int nProcs;
    int rank;
#if BITPIT_ENABLE_MPI==1
    MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#else
    nProcs = 1;
    rank   = 0;
#endif
 
    // Initialize the logger
    log::manager().initialize(log::MODE_SEPARATE, false, nProcs, rank);
    log::cout() << log::fileVerbosity(log::LEVEL_INFO);
    log::cout() << log::disableConsole();
 
    // Run the example
    try {
        run();
    } catch (const std::exception &exception) {
        log::cout() << exception.what();
        exit(1);
    }
 
#if BITPIT_ENABLE_MPI==1
    MPI_Finalize();
#endif
}