add fmmdata.py, exp_plot_fmms.py and new experimental datas

3dpp/exp/exp_plot_fmms.py

@@ -0,0 +1,25 @@
+import fmmdata
+import sys
+
+if __name__ == '__main__':
+    filenames = ['slurm-152683.out']
+    if len(sys.argv) > 1:
+        argname = sys.argv[1]
+        if argname[-4:] == '.out':
+            filenames = [argname]
+        elif argname[-4:] == '.txt':
+            file = open(argname)
+            filenames = file.read()
+            filenames = filenames.splitlines()
+            comments = [filename for filename in filenames if len(filename) > 0 and filename[0] == '#']
+            filenames = [filename for filename in filenames if len(filename) > 0 and filename[0] != '#']
+            print(*comments,sep='\n')
+        else:
+            raise 'file type not supported'
+    print(f"plot from {filenames}")
+    fmm_avg = fmmdata.get_averaged_fmms(filenames)
+    num = len(fmm_avg['filenames'])
+    p_fmms = fmmdata.collect_fmms_by_p(fmm_avg['fmms'])
+    i_fmms = fmmdata.collect_fmms_by_i(fmm_avg['fmms'])
+    fmmdata.plot_fmms(p_fmms, num, 'i-p', 'P')
+    fmmdata.plot_fmms(i_fmms, num, 'p-i', 'i')

3dpp/exp/fmmdata.py

@@ -0,0 +1,145 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.pyplot import figure
+import re
+import sys
+
+class fmmdata:
+    def __init__(self,p,i,l2p,l2f,kt,tt):
+        self.p = p
+        self.i = i
+        self.l2p = l2p
+        self.l2f = l2f
+        self.kt = kt
+        self.tt = tt
+    def __str__(self) -> str:
+        return f'FMM(P={self.p}, i={self.i}, l2p={self.l2p}, l2f={self.l2f}, kernel_time={self.kt}, total_time={self.tt})'
+
+    def __repr__(self):
+        return str(self)
+
+    def same_setting(self,fmm):
+        return self.p == fmm.p and self.i == fmm.i
+
+    def get_setting(self):
+        return (self.p, self.i)
+
+    def __add__(self, other):
+        if self.get_setting() == other.get_setting():
+            return fmmdata(self.p,self.i,(self.l2p+other.l2p),(self.l2f+other.l2f),(self.kt+other.kt),(self.tt+other.tt))
+        raise 'FMMs with different settings cannot be added'
+
+    def __radd__(self, other):
+        if other == 0:
+            return self
+        else:
+            return self.__add__(other)
+
+    def __truediv__(self, n):
+        return fmmdata(self.p,self.i,self.l2p/n,self.l2f/n,self.kt/n,self.tt/n)
+
+
+def get_fmmdatas(filename : str):
+    res = open(filename)
+    txt = res.read()
+    data1 = re.findall(r'L2  \(p\)  : (.*)L2  \(f\)',txt)
+    data2 = re.findall(r'L2  \(f\)  : (.*)L2  \(e\)',txt)
+    data3 = re.findall(r'\[FMM_kernels time measured : (.*) seconds.\]\n',txt)
+    data4 = re.findall(r'\[FMM time measured : (.*) seconds.\]\n',txt)
+    data5 = re.findall(r'images               = (.*)', txt)
+    data6 = re.findall(r'P                    = (.*)', txt)
+    L2ps = [eval(e) for e in data1]
+    L2fs = [eval(e) for e in data2]
+    kts = [eval(e) for e in data3]
+    tts = [eval(e) for e in data4]
+    imgs = [eval(fmm) for fmm in data5]
+    ps = [eval(fmm) for fmm in data6]
+    fmms = [fmmdata(p,i,l2p,l2f,kt,tt) for p,i,l2p,l2f,kt,tt in zip(ps,imgs,L2ps,L2fs,kts,tts)]
+    return {'filename':filename,'data':fmms}
+
+def collect_fmms_by_p(fmms):
+    res = {}
+    all_P = {fmm.p for fmm in fmms}
+    for p in all_P:
+        res[p] = [fmm for fmm in fmms if fmm.p == p]
+    return res
+
+def collect_fmms_by_i(fmms):
+    res = {}
+    all_i = {fmm.i for fmm in fmms}
+    for i in all_i:
+        res[i] = [fmm for fmm in fmms if fmm.i == i]
+    return res
+
+def average_fmms(fmms):
+    res = {'filenames':[],'fmms':[]}
+    fmm_grouped_by_setting = {}
+    for fmm in fmms:
+        datas = fmm['data']
+        res['filenames'].append(fmm['filename'])
+        for data in datas:
+            setting = data.get_setting()
+            if setting not in fmm_grouped_by_setting:
+                fmm_grouped_by_setting[setting] = []
+            fmm_grouped_by_setting[setting].append(data)
+    for setting in fmm_grouped_by_setting:
+        datas = fmm_grouped_by_setting[setting]
+        num = len(datas)
+        data_averaged = sum(datas) / num
+        res['fmms'].append(data_averaged)
+    return res
+
+def plot_fmms(fmms,fmm_file_num,filename,keyname):
+    color_codes = ['#002c53','#ffa510','#0c84c6','#628449','#f74d4d','#2455a4','#41b7ac']
+    markers = ['-o','-*','-x','-d']
+    plt.rcParams["font.family"] = "DejaVu Serif" 
+    fig = plt.figure()
+    fig.set_figheight(4)
+    fig.set_figwidth(11)
+    ax1 = fig.add_subplot(1, 2, 1)
+    ax2 = fig.add_subplot(1, 2, 2)
+    xs = []
+    cnt = 0
+    if keyname == 'P' or keyname == 'p':
+        xname = 'images'
+    else:
+        xname = 'P'
+    for key in fmms:
+        c = color_codes[cnt % len(color_codes)]
+        marker = markers[int(cnt / len(color_codes)) % len(markers)]
+        fmm = fmms[key]
+        if keyname == 'P' or keyname == 'p':
+            x = [fmm.i for fmm in fmm]
+        else:
+            x = [fmm.p for fmm in fmm]
+        if len(x) > len(xs):
+            xs = x
+        l2f = [fmm.l2f for fmm in fmm]
+        kt = [fmm.kt for fmm in fmm]
+        ax1.plot(x,l2f,marker,label=f'{keyname}={key}',color=c)
+        ax2.plot(x,kt,marker,label=f'{keyname}={key}',color=c)
+        cnt+=1
+
+    ax1.set_xticks(xs)
+    ax1.set_yscale('log')
+    ax1.set_xlabel(xname)
+    ax1.set_ylabel('L2-err')
+    ax1.grid()
+    ax1.set_title(f'L2 error of Force(avg_num={fmm_file_num})')
+    ax1.legend()
+    ax2.set_xticks(xs)
+    ax2.set_xlabel(xname)
+    ax2.set_ylabel('time(s)')
+    ax2.legend()
+    ax2.grid()
+    ax2.set_title(f'kernel time(avg_num={fmm_file_num})')
+    ax2.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+
+    fig.tight_layout()
+    plt.savefig('exp(' + filename + ").png", dpi=250)
+
+def get_fmms(filenames : list):
+    return [get_fmmdatas(filename) for filename in filenames]
+
+def get_averaged_fmms(filenames : list):
+    return average_fmms(get_fmms(filenames))

3dpp/exp/results1.txt

@@ -0,0 +1,5 @@
+#exp on PBC
+#machine = rtx6000-ada_1
+
+slurm-152683.out
+slurm-152683.out

3dpp/src/argument.cpp

@@ -13,7 +13,7 @@ rtfmm::Argument::Argument(int argc, char* argv[])
     cmd.add<int>("ncrit", 'm', "minimum number of bodies per leaf box", false, 16);
     cmd.add<int>("timing", 't', "if measure execution time of computational steps", false, 1);
     cmd.add<int>("images", 'i', "periodic images depth", false, 0, cmdline::range(0, 20));
-    cmd.add<int>("verbose", 'v', "is verbose", false, 1);
+    cmd.add<int>("verbose", 'v', "is verbose", false, 0);
     cmd.add<real>("cycle", 'c', "cycle of images(or box_r)", false, 2 * M_PI);
     cmd.add<int>("num_compare", 0, "compare number of target body(also direct calculation number)", false, -1);
     cmd.add<int>("ewald_ksize", 0, "ksize of ewald DFT", false, 11);

3dpp/src/fmm.cpp

@@ -13,20 +13,24 @@ rtfmm::LaplaceFMM::LaplaceFMM(const Bodies3& bs_, const Argument& args_)
 rtfmm::Bodies3 rtfmm::LaplaceFMM::solve()
 {
     /* build tree */
-    tbegin(tree_traverse);
+    tbegin(build_and_traverse);
+    tbegin(build_tree);
     Tree tree;
     tree.build(bs, args.x, args.r, args.ncrit, Tree::TreeType::nonuniform);
     cs = tree.get_cells();
     if(verbose && args.check_tree) check_tree(cs);
+    tend(build_tree);
 
     /* traverse to get interaction list */
+    tbegin(traverse);
     traverser.traverse(tree, args.cycle, args.images);
     cs = traverser.get_cells();
     if(verbose && args.check_tree) check_traverser(traverser);
     if(verbose && args.check_tree) check_cells(cs);
     tree_depth_range = get_min_max_depth(cs);
     if(verbose) std::cout<<"tree_depth_range = "<<tree_depth_range<<std::endl;
-    tend(tree_traverse);
+    tend(traverse);
+    tend(build_and_traverse);
 
     tbegin(init_cell_matrix);
     init_cell_matrix(cs);

3dpp/src/traverser.cpp

@@ -27,12 +27,20 @@ void rtfmm::Traverser::traverse(Tree& tree, real cycle, int images)
     }
     else if(images >= 1)
     {
+        tbegin(horizontal_periodic_near);
         horizontal_periodic_near(cycle);
+        tend(horizontal_periodic_near);
+        tbegin(make_M2L_parent_map_i1);
         make_M2L_parent_map_i1(cycle); // for t
+        tend(make_M2L_parent_map_i1);
         if(images >= 2)
         {
+            tbegin(horizontal_periodic_far);
             horizontal_periodic_far(cycle,images);
+            tend(horizontal_periodic_far);
+
             // for t
+            tbegin(t_image2);
             int idx = cells.size() - 1 - images;
             for(int m = 1; m < images; m++)
             {
@@ -54,6 +62,7 @@ void rtfmm::Traverser::traverse(Tree& tree, real cycle, int images)
                     }
                 }
             }
+            tend(t_image2);
         }
     }
 }
@@ -242,6 +251,7 @@ void rtfmm::Traverser::make_M2L_parent_map()
 void rtfmm::Traverser::make_M2L_parent_map_i1(real cycle)
 {
     make_M2L_parent_map();
+    #pragma omp parallel for
     for(int j = 0; j < cells.size(); j++)
     {
         Cell3& ctar = cells[j];