PABLO_example_00008.cpp

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

Parallel 2D adaptive mesh refinement (AMR) with data using PABLOExample 00007 is here enriched by an additional feature.

The load-balance is performed by keep a family of desired level together on the same process.

In particular in this example the families, containing elements from the maximum level reached in the quadtree to that level minus 3, are kept on the same partition.

To run: ./PABLO_example_00008
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 pablo8(2);
 
    int idx = 0;
    pablo8.setBalance(idx,false);
 
    for (iter=1; iter<6; iter++){
        pablo8.adaptGlobalRefine();
    }
 
    double xc, yc;
    xc = yc = 0.5;
    double radius = 0.25;
 
    uint32_t nocts = pablo8.getNumOctants();
    uint32_t nghosts = pablo8.getNumGhosts();
    vector<double> oct_data(nocts, 0.0), ghost_data(nghosts, 0.0);
 
    for (unsigned int i=0; i<nocts; i++){
        vector<array<double,3> > nodes = pablo8.getNodes(i);
        array<double,3> center = pablo8.getCenter(i);
        for (int j=0; j<4; j++){
            double x = nodes[j][0];
            double y = nodes[j][1];
            if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                oct_data[i] = (pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0));
            }
        }
    }
 
    iter = 0;
    pablo8.updateConnectivity();
    pablo8.writeTest("pablo00008_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);
 
    int start = 1;
    for (iter=start; iter<start+2; iter++){
        for (unsigned int i=0; i<nocts; i++){
            vector<array<double,3> > nodes = pablo8.getNodes(i);
            array<double,3> center = pablo8.getCenter(i);
            for (int j=0; j<4; j++){
                double x = nodes[j][0];
                double y = nodes[j][1];
                if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                    if (center[0]<=xc){
 
                        pablo8.setMarker(i,1);
                    }
                    else{
 
                        pablo8.setMarker(i,-1);
                    }
                }
            }
        }
 
        vector<double> oct_data_new;
        vector<uint32_t> mapper;
        vector<bool> isghost;
        pablo8.adapt(true);
        nocts = pablo8.getNumOctants();
        oct_data_new.resize(nocts, 0.0);
 
        for (uint32_t i=0; i<nocts; i++){
            pablo8.getMapping(i, mapper, isghost);
            if (pablo8.getIsNewC(i)){
                for (int j=0; j<4; j++){
                    if (isghost[j]){
                        oct_data_new[i] += ghost_data[mapper[j]]/4;
                    }
                    else{
                        oct_data_new[i] += oct_data[mapper[j]]/4;
                    }
                }
            }
            else{
                oct_data_new[i] += oct_data[mapper[0]];
            }
        }
 
        pablo8.updateConnectivity();
        pablo8.writeTest("pablo00008_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data_new);
 
        oct_data = oct_data_new;
    }
 
#if BITPIT_ENABLE_MPI==1
    uint8_t levels = 4;
    UserDataLB<vector<double> > data_lb(oct_data,ghost_data);
    pablo8.loadBalance(data_lb, levels);
#endif
 
    pablo8.updateConnectivity();
    pablo8.writeTest("pablo00008_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
}