<commit_before>#include <iostream>
#include <iomanip>
#include <sstream>
#include <thrust/host_vector.h>

#include "dg/algorithm.h"
#include "dg/file/json_utilities.h"

#ifndef WITHOUT_GLFW
#include "draw/host_window.h"
#endif
#include "dg/file/file.h"

#include "init.h"
#include "diag.h"
#include "shu.cuh"


int main( int argc, char* argv[])
{
    ////Parameter initialisation ////////////////////////////////////////////
    Json::Value js;
    enum dg::file::error mode = dg::file::error::is_throw;
    if( argc == 1)
        dg::file::file2Json( "input/default.json", js, dg::file::comments::are_discarded);
    else
        dg::file::file2Json( argv[1], js);
    std::cout << js <<std::endl;

    /////////////////////////////////////////////////////////////////
    dg::CartesianGrid2d grid = shu::createGrid( js, mode);
    dg::DVec w2d( dg::create::weights(grid));
    /////////////////////////////////////////////////////////////////

    std::string initial = dg::file::get( mode, js, "init", "type", "lamb").asString();
    dg::HVec omega = shu::initial_conditions.at( initial)(js, mode, grid);

    dg::DVec y0( omega ), y1( y0);
    //subtract mean mass
    if( grid.bcx() == dg::PER && grid.bcy() == dg::PER)
    {
        double meanMass = dg::blas1::dot( y0, w2d)/(double)(grid.lx()*grid.ly());
        dg::blas1::axpby( -meanMass, 1., 1., y0);
    }
    //make solver and stepper
    shu::Shu<dg::CartesianGrid2d, dg::DMatrix, dg::DVec>
        shu( grid, js, mode);
    shu::Diffusion<dg::CartesianGrid2d, dg::DMatrix, dg::DVec> diffusion( grid, js, mode);
    if( "mms" == initial)
    {
        double sigma = dg::file::get( mode, js, "init", "sigma", 0.2).asDouble();
        double velocity = dg::file::get( mode, js, "init", "velocity", 0.1).asDouble();
        shu.set_mms_source( sigma, velocity, grid.ly());
    }

    double time = 0;
    shu::Variables var = {shu, grid, y0, time, w2d, 0., mode, js};
    dg::Timer t;
    t.tic();
    shu( 0., y0, y1);
    t.toc();
    var.duration = t.diff();
    for( auto record : shu::diagnostics1d_list)
    {
        double result = record.function( var);
        std::cout  << "Diagnostics "<<record.name<<" "<<result<<"\n";
    }

    /// ////////////// Initialize timestepper ///////////////////////
    std::string stepper = dg::file::get( mode, js, "timestepper", "type", "FilteredExplicitMultistep").asString();
    std::string regularization = dg::file::get( mode, js, "regularization", "type", "moddal").asString();
    dg::ModalFilter<dg::DMatrix, dg::DVec> filter;
    dg::IdentityFilter id;
    bool apply_filter = true;
    dg::ImExMultistep<dg::DVec> imex;
    dg::ShuOsher<dg::DVec> shu_osher;
    dg::FilteredExplicitMultistep<dg::DVec> multistep;
    if( regularization == "modal")
    {
        double alpha = dg::file::get( mode, js, "regularization", "alpha", 36).asDouble();
        double order = dg::file::get( mode, js, "regularization", "order", 8).asDouble();
        double eta_c = dg::file::get( mode, js, "regularization", "eta_c", 0.5).asDouble();
        filter.construct( dg::ExponentialFilter(alpha, eta_c, order, grid.n()), grid);
    }
    else
        apply_filter = false;

    double dt = dg::file::get( mode, js, "timestepper", "dt", 2e-3).asDouble();
    if( "ImExMultistep" == stepper)
    {
        if( regularization != "viscosity")
        {
            throw dg::Error(dg::Message(_ping_)<<"Error: ImExMultistep only works with viscosity regularization! Exit now!");

            return -1;
        }
        double eps_time = dg::file::get( mode, js, "timestepper", "eps_time", 1e-10).asDouble();
        std::string tableau = dg::file::get( mode, js, "timestepper", "tableau", "ImEx-BDF-3-3").asString();
        imex.construct( tableau, y0, y0.size(), eps_time);
        imex.init( shu, diffusion, time, y0, dt);
    }
    else if( "Shu-Osher" == stepper)
    {
        shu_osher.construct( "SSPRK-3-3", y0);
    }
    else if( "FilteredExplicitMultistep" == stepper)
    {
        multistep.construct( "eBDF-3-3", y0);
        if( apply_filter)
            multistep.init( shu, filter, time, y0, dt);
        else
            multistep.init( shu, id, time, y0, dt);
    }
    else
    {
        throw dg::Error(dg::Message(_ping_)<<"Error! Timestepper not recognized!\n");

        return -1;
    }
    unsigned maxout = dg::file::get( mode, js, "output", "maxout", 100).asUInt();
    unsigned itstp = dg::file::get( mode, js, "output", "itstp", 5).asUInt();
    std::string output = dg::file::get( mode, js, "output", "type", "glfw").asString();
#ifndef WITHOUT_GLFW
    if( "glfw" == output)
    {
        ////////////////////////////////glfw//////////////////////////////
        //create visualisation vectors
        dg::DVec visual( grid.size());
        dg::HVec hvisual( grid.size());
        //transform vector to an equidistant grid
        std::stringstream title;
        GLFWwindow* w = draw::glfwInitAndCreateWindow(600, 600, "");
        draw::RenderHostData render( 1,1);
        draw::ColorMapRedBlueExt colors( 1.);
        dg::IDMatrix equidistant = dg::create::backscatter( grid );
        while (!glfwWindowShouldClose(w) && time < maxout*itstp*dt)
        {
            dg::blas2::symv( equidistant, y0, visual);
            colors.scale() =  (float)thrust::reduce( visual.begin(), visual.end(), -1., dg::AbsMax<double>() );
            //draw and swap buffers
            dg::blas1::transfer( visual, hvisual);
            render.renderQuad( hvisual, grid.n()*grid.Nx(), grid.n()*grid.Ny(), colors);
            title << "Time "<<time<< " \ttook "<<t.diff()/(double)itstp<<"\t per step";
            glfwSetWindowTitle(w, title.str().c_str());
            title.str("");
            glfwPollEvents();
            glfwSwapBuffers(w);
            //step
            t.tic();
            try{
                for( unsigned j=0; j<itstp; j++)
                {
                    if( "ImExMultistep" == stepper)
                        imex.step( shu, diffusion, time, y0);
                    else if ( "FilteredExplicitMultistep" == stepper)
                    {
                        if( apply_filter)
                            multistep.step( shu, filter, time, y0);
                        else
                            multistep.step( shu, id, time, y0);
                    }
                    else if ( "Shu-Osher" == stepper)
                    {
                        if( apply_filter)
                            shu_osher.step( shu, filter, time, y0, time, y0, dt);
                        else
                            shu_osher.step( shu, id, time, y0, time, y0, dt);
                    }
                }
            } catch( dg::Fail& fail) {
                std::cerr << "CG failed to converge to "<<fail.epsilon()<<"\n";
                std::cerr << "Does Simulation respect CFL condition?\n";
                return -1;
            }
            t.toc();
            var.duration = t.diff()/(double)itstp;
        }
        glfwTerminate();
    }
    else if( "netcdf" == output)
#endif //WITHOUT_GLFW
    {
        std::string inputfile = js.toStyledString(); //save input without comments, which is important if netcdf file is later read by another parser
        std::string outputfile;
        if( argc == 1 || argc == 2)
            outputfile = "shu.nc";
        else
            outputfile = argv[2];
        /// //////////////////////set up netcdf/////////////////////////////////////
        dg::file::NC_Error_Handle err;
        int ncid=-1;
        try{
            err = nc_create( outputfile.c_str(),NC_NETCDF4|NC_CLOBBER, &ncid);
        }catch( std::exception& e)
        {
            std::cerr << "ERROR creating file "<<outputfile<<std::endl;<commit_msg>Fix missing if bug in shu_hpc<commit_after>#include <iostream>
#include <iomanip>
#include <sstream>
#include <thrust/host_vector.h>

#include "dg/algorithm.h"
#include "dg/file/json_utilities.h"

#ifndef WITHOUT_GLFW
#include "draw/host_window.h"
#endif
#include "dg/file/file.h"

#include "init.h"
#include "diag.h"
#include "shu.cuh"


int main( int argc, char* argv[])
{
    ////Parameter initialisation ////////////////////////////////////////////
    Json::Value js;
    enum dg::file::error mode = dg::file::error::is_throw;
    if( argc == 1)
        dg::file::file2Json( "input/default.json", js, dg::file::comments::are_discarded);
    else
        dg::file::file2Json( argv[1], js);
    std::cout << js <<std::endl;

    /////////////////////////////////////////////////////////////////
    dg::CartesianGrid2d grid = shu::createGrid( js, mode);
    dg::DVec w2d( dg::create::weights(grid));
    /////////////////////////////////////////////////////////////////

    std::string initial = dg::file::get( mode, js, "init", "type", "lamb").asString();
    dg::HVec omega = shu::initial_conditions.at( initial)(js, mode, grid);

    dg::DVec y0( omega ), y1( y0);
    //subtract mean mass
    if( grid.bcx() == dg::PER && grid.bcy() == dg::PER)
    {
        double meanMass = dg::blas1::dot( y0, w2d)/(double)(grid.lx()*grid.ly());
        dg::blas1::axpby( -meanMass, 1., 1., y0);
    }
    //make solver and stepper
    shu::Shu<dg::CartesianGrid2d, dg::DMatrix, dg::DVec>
        shu( grid, js, mode);
    shu::Diffusion<dg::CartesianGrid2d, dg::DMatrix, dg::DVec> diffusion( grid, js, mode);
    if( "mms" == initial)
    {
        double sigma = dg::file::get( mode, js, "init", "sigma", 0.2).asDouble();
        double velocity = dg::file::get( mode, js, "init", "velocity", 0.1).asDouble();
        shu.set_mms_source( sigma, velocity, grid.ly());
    }

    double time = 0;
    shu::Variables var = {shu, grid, y0, time, w2d, 0., mode, js};
    dg::Timer t;
    t.tic();
    shu( 0., y0, y1);
    t.toc();
    var.duration = t.diff();
    for( auto record : shu::diagnostics1d_list)
    {
        double result = record.function( var);
        std::cout  << "Diagnostics "<<record.name<<" "<<result<<"\n";
    }

    /// ////////////// Initialize timestepper ///////////////////////
    std::string stepper = dg::file::get( mode, js, "timestepper", "type", "FilteredExplicitMultistep").asString();
    std::string regularization = dg::file::get( mode, js, "regularization", "type", "moddal").asString();
    dg::ModalFilter<dg::DMatrix, dg::DVec> filter;
    dg::IdentityFilter id;
    bool apply_filter = true;
    dg::ImExMultistep<dg::DVec> imex;
    dg::ShuOsher<dg::DVec> shu_osher;
    dg::FilteredExplicitMultistep<dg::DVec> multistep;
    if( regularization == "modal")
    {
        double alpha = dg::file::get( mode, js, "regularization", "alpha", 36).asDouble();
        double order = dg::file::get( mode, js, "regularization", "order", 8).asDouble();
        double eta_c = dg::file::get( mode, js, "regularization", "eta_c", 0.5).asDouble();
        filter.construct( dg::ExponentialFilter(alpha, eta_c, order, grid.n()), grid);
    }
    else
        apply_filter = false;

    double dt = dg::file::get( mode, js, "timestepper", "dt", 2e-3).asDouble();
    if( "ImExMultistep" == stepper)
    {
        if( regularization != "viscosity")
        {
            throw dg::Error(dg::Message(_ping_)<<"Error: ImExMultistep only works with viscosity regularization! Exit now!");

            return -1;
        }
        double eps_time = dg::file::get( mode, js, "timestepper", "eps_time", 1e-10).asDouble();
        std::string tableau = dg::file::get( mode, js, "timestepper", "tableau", "ImEx-BDF-3-3").asString();
        imex.construct( tableau, y0, y0.size(), eps_time);
        imex.init( shu, diffusion, time, y0, dt);
    }
    else if( "Shu-Osher" == stepper)
    {
        shu_osher.construct( "SSPRK-3-3", y0);
    }
    else if( "FilteredExplicitMultistep" == stepper)
    {
        multistep.construct( "eBDF-3-3", y0);
        if( apply_filter)
            multistep.init( shu, filter, time, y0, dt);
        else
            multistep.init( shu, id, time, y0, dt);
    }
    else
    {
        throw dg::Error(dg::Message(_ping_)<<"Error! Timestepper not recognized!\n");

        return -1;
    }
    unsigned maxout = dg::file::get( mode, js, "output", "maxout", 100).asUInt();
    unsigned itstp = dg::file::get( mode, js, "output", "itstp", 5).asUInt();
    std::string output = dg::file::get( mode, js, "output", "type", "glfw").asString();
#ifndef WITHOUT_GLFW
    if( "glfw" == output)
    {
        ////////////////////////////////glfw//////////////////////////////
        //create visualisation vectors
        dg::DVec visual( grid.size());
        dg::HVec hvisual( grid.size());
        //transform vector to an equidistant grid
        std::stringstream title;
        GLFWwindow* w = draw::glfwInitAndCreateWindow(600, 600, "");
        draw::RenderHostData render( 1,1);
        draw::ColorMapRedBlueExt colors( 1.);
        dg::IDMatrix equidistant = dg::create::backscatter( grid );
        while (!glfwWindowShouldClose(w) && time < maxout*itstp*dt)
        {
            dg::blas2::symv( equidistant, y0, visual);
            colors.scale() =  (float)thrust::reduce( visual.begin(), visual.end(), -1., dg::AbsMax<double>() );
            //draw and swap buffers
            dg::blas1::transfer( visual, hvisual);
            render.renderQuad( hvisual, grid.n()*grid.Nx(), grid.n()*grid.Ny(), colors);
            title << "Time "<<time<< " \ttook "<<t.diff()/(double)itstp<<"\t per step";
            glfwSetWindowTitle(w, title.str().c_str());
            title.str("");
            glfwPollEvents();
            glfwSwapBuffers(w);
            //step
            t.tic();
            try{
                for( unsigned j=0; j<itstp; j++)
                {
                    if( "ImExMultistep" == stepper)
                        imex.step( shu, diffusion, time, y0);
                    else if ( "FilteredExplicitMultistep" == stepper)
                    {
                        if( apply_filter)
                            multistep.step( shu, filter, time, y0);
                        else
                            multistep.step( shu, id, time, y0);
                    }
                    else if ( "Shu-Osher" == stepper)
                    {
                        if( apply_filter)
                            shu_osher.step( shu, filter, time, y0, time, y0, dt);
                        else
                            shu_osher.step( shu, id, time, y0, time, y0, dt);
                    }
                }
            } catch( dg::Fail& fail) {
                std::cerr << "CG failed to converge to "<<fail.epsilon()<<"\n";
                std::cerr << "Does Simulation respect CFL condition?\n";
                return -1;
            }
            t.toc();
            var.duration = t.diff()/(double)itstp;
        }
        glfwTerminate();
    }
#endif //WITHOUT_GLFW
    if( "netcdf" == output)
    {
        std::string inputfile = js.toStyledString(); //save input without comments, which is important if netcdf file is later read by another parser
        std::string outputfile;
        if( argc == 1 || argc == 2)
            outputfile = "shu.nc";
        else
            outputfile = argv[2];
        /// //////////////////////set up netcdf/////////////////////////////////////
        dg::file::NC_Error_Handle err;
        int ncid=-1;
        try{
            err = nc_create( outputfile.c_str(),NC_NETCDF4|NC_CLOBBER, &ncid);
        }catch( std::exception& e)
        {
            std::cerr << "ERROR creating file "<<outputfile<<std::endl;
{'fraction_non_alphanumeric': 0.20248919842651705, 'fraction_numerical': 0.011091765009350616, 'mean_word_length': 1.9657678332377126, 'pattern_counts': {'":': 0, '<': 113, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4406725', 'n_tokens_mistral': 5223, 'n_tokens_neox': 5003, 'n_words': 1135}