io.h

Go to the documentation of this file.

//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //
//                                                                                   //
//                       / /   / /    __    \ \  / /                                 //
//                      / /__ / /   / _ \    \ \/ /                                  //
//                     /  ___  /   | |/_/    / /\ \                                  //
//                    / /   / /    \_\      / /  \ \                                 //
//                                                                                   //
//                                                                                   //
//  Copyright (c) 2016, Jakub Benda, Charles University in Prague                    //
//                                                                                   //
// MIT License:                                                                      //
//                                                                                   //
//  Permission is hereby granted, free of charge, to any person obtaining a          //
// copy of this software and associated documentation files (the "Software"),        //
// to deal in the Software without restriction, including without limitation         //
// the rights to use, copy, modify, merge, publish, distribute, sublicense,          //
// and/or sell copies of the Software, and to permit persons to whom the             //
// Software is furnished to do so, subject to the following conditions:              //
//                                                                                   //
//  The above copyright notice and this permission notice shall be included          //
// in all copies or substantial portions of the Software.                            //
//                                                                                   //
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS          //
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE       //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, //
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF         //
// OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.  //
//                                                                                   //
//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //

#ifndef HEX_IO_H
#define HEX_IO_H

#include <cstdio>
#include <fstream>
#include <vector>
#include <string>

#include "hex-arrays.h"

#include "bspline.h"
#include "luft.h"
#include "parallel.h"

class CommandLine
{
    public:
        
        typedef enum {
            StgNone    = 0x00,
            StgRadial  = 0x01,
            StgSolve   = 0x02,
            StgExtract = 0x04
        } HexEcsStg;
        
        // constructor
        CommandLine (int argc, char* argv[])
            : writegrid(false), zipdata(), parallel(false), preconditioner("ILU"),
              droptol(1e-8), itinerary(StgNone), outofcore(false), cont(false), wholematrix(false), cache_all_radint(true), cache_own_radint(true),
              itertol(1e-8), prec_itertol(1e-8), parallel_precondition(false), gpu_large_data(false),
              lightweight_full(false), lightweight_radial_cache(false), shared_scratch(false), reuse_dia_blocks(false),
              kpa_simple_rad(false), ocl_platform(0), ocl_device(0), factorizer("umfpack"), groupsize(1),
              parallel_factorization(false), parallel_extraction(true), ilu_max_iter(10), max_sub_iter(0), fail_on_sub_iter(true),
              carry_initial_guess(false), gpu_multiply(false), extract_extrapolate(false), extract_rho(-1), extract_rho_begin(-1), extract_samples(-1),
              refine_solution(false), map_solution(), map_solution_target(), ssor(-1), noluupdate(false), coupling_limit(1000),
              gpu_host_multiply(false), mumps_outofcore(false), mumps_verbose(0), kpa_drop(-1), write_intermediate_solutions(false),
              fast_bessel(false), hyb_additional_levels(0), multigrid_depth(0), multigrid_coarse_prec(0), dom_x_panels(1), dom_y_panels(1),
              dom_preconditioner("ILU"), dom_sweeps(-1), scratch(std::getenv("SCRATCHDIR") ? std::getenv("SCRATCHDIR") : "."), analytic_eigenstates(false),
              runtime_postprocess(false), sub_prec_verbose(false)
        {
            // get command line options
            parse(argc, argv);
            
            // run the whole sequence if nothing specified
            if (itinerary == StgNone)
                itinerary = StgRadial | StgSolve | StgExtract;
        }
        
        void parse (int argc, char* argv[]);
        
        //
        // public attributes
        //
        
        std::ifstream inputfile;
        
        // Write grid layout to a VTK file.
        bool writegrid;
        
        struct s_zipdata
        {
            // default constructor
            s_zipdata () : file(), Xmin(-1), Ymin(-1), Xmax(-1), Ymax(-1), nX(-1), nY(-1) {}
            
            std::string file;
            Real Xmin, Ymin, Xmax, Ymax;
            int nX, nY;
        }
        zipdata;
        
        bool parallel;
        
        std::string preconditioner;
        
        Real droptol;
        
        int itinerary;
        
        bool outofcore;
        
        bool cont;
        
        bool wholematrix;
        
        bool cache_all_radint;
        bool cache_own_radint;
        
        Real itertol;
        
        Real prec_itertol;
        
        bool parallel_precondition;
        
        bool gpu_large_data;
        
        bool lightweight_full, lightweight_radial_cache;
        
        bool shared_scratch;
        
        bool reuse_dia_blocks;
        
        bool kpa_simple_rad;
        
        unsigned ocl_platform;
        
        unsigned ocl_device;
        
        std::string factorizer;
        
        int groupsize;
        
        int parallel_factorization;
        
        int parallel_extraction;
        
        int ilu_max_iter;
        
        int max_sub_iter;
        
        bool fail_on_sub_iter;
        
        bool carry_initial_guess;
        
        bool gpu_multiply;
        
        bool extract_extrapolate;
        
        Real extract_rho;
        
        Real extract_rho_begin;
        
        int extract_samples;
        
        bool refine_solution;
        
        std::vector<std::string> map_solution;
        
        std::string map_solution_target;
        
        Real ssor;
        
        bool noluupdate;
        
        int coupling_limit;
        
        bool gpu_host_multiply;
        
        bool mumps_outofcore;
        
        int mumps_verbose;
        
        Real kpa_drop;
        
        bool write_intermediate_solutions;
        
        bool fast_bessel;
        
        int hyb_additional_levels;
        
        int multigrid_depth;
        
        int multigrid_coarse_prec;
        
        int dom_x_panels, dom_y_panels;
        
        std::string dom_preconditioner;
        
        int dom_sweeps;
        
        std::string scratch;
        
        bool analytic_eigenstates;
        
        bool runtime_postprocess;
        
        bool sub_prec_verbose;
};

class InputFile
{
    public:
        
        // constructor
        InputFile (std::ifstream & inputfile)
        {
            // read inputfile
            read(inputfile);
            
            // maximal angular momentum (the same for both parities)
            maxell = levels + L;
        }
        
        // read data from file
        void read (std::ifstream & inputfile);
        
        //
        // public attributes
        //
        
            // B-spline order
            int order;
            
            // initial atomic principal quantum number
            int ni;
            
            // total angular momentum
            int L;
            
            // total parity
            int Pi;
            
            // 'nL', the limit on number of coupled angular states;
            // there will be 'nL * (L + 1 - Pi)' coupled angular state pairs
            int levels;
            
            // upper limit on the smaller of the one-electron angular momenta
            int limit;
            
            // maximal one-electron orbital quantum number
            int maxell;
            
            // total spins to calculate
            iArray Spin;
            
            // ECS rotation angle
            Real ecstheta;
            
            // weak magnetic field in atomic units (involved only perturbatively)
            Real B;
            
            // maximal energy (Ry) of states included in the asymptotic (outer) region
            Real channel_max_E;
            
            // real B-spline knots
            rArray rknots;
            
            // real B-spline knot projectile extention
            rArray rknots_ext;
            
            // complex-to-become knots (after rotation)
            rArray cknots;
            
            // total energies for which to solve the system
            rArray Etot;
            
            // maximal total energy
            Real max_Etot;
            
            // initial and final atomic states
            std::vector<std::tuple<int,int,int>> instates, outstates;
            
            // atom charge (must be positive integer)
            Real Za;
            
            // projectile charge (only sign)
            Real Zp;
            
            // whether to calculate just the inner problem (decided from the knot sequence)
            bool inner_only;
            
            // keep only l1 <= l2; this is useful for large total angular momenta for reduction of the angular basis
            bool exchange;
};

class SolutionIO
{
    public:
        
        SolutionIO () {}
        
        SolutionIO
        (
            int L, int S, int Pi,
            int ni, int li, int mi,
            Real E,
            std::vector<std::pair<int,int>> const & ang,
            std::vector<std::pair<int,int>> const & chann = std::vector<std::pair<int,int>>(),
            std::string prefix = "psi"
        ) : L_(L), S_(S), Pi_(Pi),
            ni_(ni), li_(li), mi_(mi),
            E_(E), ang_(ang), prefix_(prefix)
        {
            if (chann.empty())
            {
                // initialize channel counts to zero
                chann_.resize(ang.size(), std::make_pair(0,0));
            }
            else
            {
                // check consistency of the parameters
                if (ang.size() != chann.size())
                    HexException("Incompatible size of angular momentum pairs list and channel count list.");
                
                // use the supplied channel counts
                chann_ = chann;
            }
        }
        
        static const int All = -1;
        
        std::string name (int ill = SolutionIO::All) const
        {
            if (ill == -1)
                return format("%s-%g-%d-%d-%d-%d-%d-%d.hdf", prefix_.c_str(), E_ + 1, L_, S_, Pi_, ni_, li_, mi_);
            else
                return format("%s-%g-%d-%d-%d-%d-%d-%d-(%d,%d).hdf", prefix_.c_str(), E_ + 1, L_, S_, Pi_, ni_, li_, mi_, ang_[ill].first, ang_[ill].second);
        }
        
        bool check (int ill, std::size_t & total_size) const
        {
            // look for specific solution segment file
            if (ill != SolutionIO::All)
            {
                HDFFile fsingle (name(ill), HDFFile::readonly);
                return fsingle.valid() ? fsingle.size("array")/2 : 0;
            }
            
            // look for all solution segment files
            std::vector<std::size_t> size (ang_.size());
            for (unsigned illp = 0; illp < ang_.size(); illp++)
            {
                HDFFile fsingle (name(illp), HDFFile::readonly);
                size[illp] = (fsingle.valid() ? fsingle.size("array")/2 : 0);
            }
            
            // calculate total size
            total_size = std::accumulate(size.begin(), size.end(), 0);
            
            // check that all blocks existed
            return std::find(size.begin(), size.end(), 0) == size.end();
        }
        bool check (int ill = SolutionIO::All) const
        {
            std::size_t total_size = 0;
            return check(ill, total_size);
        }
        
        bool load (BlockArray<Complex> & solution, int ill = SolutionIO::All)
        {
            // check that the requested files are present
            if (not check(ill))
                return false;
            
            // select segments
            iArray segments_to_load = { ill };
            if (ill == SolutionIO::All)
                segments_to_load = linspace<int>(0, solution.size() - 1, solution.size());
            
            // for all blocks to load
            for (int i : segments_to_load)
            {
                //
                // Either read from the whole monolithic solution file ...
                //
                
                    HDFFile hdf (name(), HDFFile::readonly);
                    if (hdf.valid())
                    {
                        // get segment size
                        std::size_t segsize = hdf.size("array") / ang_.size();
                        
                        // allocate memory
                        solution[i].resize(segsize);
                        
                        // read data from the file
                        if (not hdf.read("array", solution[i].data(), segsize, i * segsize))
                            return false;
                        
                        // dump to disk
                        if (not solution.inmemory())
                        {
                            solution.hdfsave(i);
                            solution[i].drop();
                        }
                    }
                
                //
                // ... or read from a solution segment.
                //
                
                    else
                    {
                        // simply load the requested solution segment file
                        if (not solution[i].hdfload(name(i)))
                            return false;
                        
                        // dump to disk
                        if (not solution.inmemory())
                        {
                            solution.hdfsave(i);
                            solution[i].drop();
                        }
                        
                        // read numbers of channels
                        HDFFile hdf (name(i), HDFFile::readonly);
                        if (hdf.valid())
                        {
                            hdf.read("Nchan1", &chann_[i].first,  1);
                            hdf.read("Nchan2", &chann_[i].second, 1);
                        }
                    }
            }
            
            return true;
        }
        
        bool save (BlockArray<Complex> const & solution, int ill = SolutionIO::All) const
        {
            // select segments
            iArray segments_to_save = { ill };
            if (ill == SolutionIO::All)
                segments_to_save = linspace<int>(0, solution.size() - 1, solution.size());
            
            // for all segments
            for (int iseg : segments_to_save)
            {
                // load data
                if (not solution.inmemory())
                    const_cast<BlockArray<Complex>&>(solution).hdfload(iseg);
                
                // write data
                if (not save_segment(solution[iseg], iseg))
                    return false;
                
                // unload data
                if (not solution.inmemory())
                    const_cast<BlockArray<Complex>&>(solution)[iseg].drop();
            }
            
            // success
            return true;
        }
        
        bool save_segment (const cArrayView solution, int ill) const
        {
            // open file
            HDFFile hdf (name(ill), HDFFile::overwrite);
            
            // check file status
            if (not hdf.valid())
                return false;
            
            // write data
            bool success = true;
            success = (success and hdf.write("L",  &L_,  1));
            success = (success and hdf.write("S",  &S_,  1));
            success = (success and hdf.write("Pi", &Pi_, 1));
            success = (success and hdf.write("ni", &ni_, 1));
            success = (success and hdf.write("li", &li_, 1));
            success = (success and hdf.write("mi", &mi_, 1));
            success = (success and hdf.write("E",  &E_,  1));
            success = (success and hdf.write("l1",  &ang_[ill].first,  1));
            success = (success and hdf.write("l2",  &ang_[ill].second, 1));
            success = (success and hdf.write("Nchan1", &chann_[ill].first, 1));
            success = (success and hdf.write("Nchan2", &chann_[ill].second, 1));
            success = (success and hdf.write("array", solution.data(), solution.size()));
            return success;
        }
        
        std::vector<std::pair<int,int>> const & channels () const { return chann_; }
    
    private:
        
        int L_, S_, Pi_, ni_, li_, mi_;
        Real E_;
        std::vector<std::pair<int,int>> ang_, chann_;
        std::string prefix_;
};

void zip_solution
(
    CommandLine const & cmd,
    InputFile const & inp,
    Parallel const & par, 
    Bspline const & bspline_inner,
    Bspline const & bspline_full,
    std::vector<std::pair<int,int>> const & ll
);

void write_grid
(
    Bspline const & bspline,
    std::string const & basename
);

#endif /* HEX_IO_H */