add gaussian explosion code

Renaming host layer functions and introducing acSetVertexBuffer

Approved-by: Miikka Väisälä <mvaisala@asiaa.sinica.edu.tw>

CMakeLists.txt

@@ -11,7 +11,7 @@ project(astaroth C CXX CUDA)
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR})
 
 ## Project-wide compilation flags
-set(COMMON_FLAGS "-mavx -DOMPI_SKIP_MPICXX -Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion -Wshadow") # -DOMPI_SKIP_MPICXX is to force OpenMPI to use the C interface
+set(COMMON_FLAGS "-DOMPI_SKIP_MPICXX -Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion -Wshadow") # -DOMPI_SKIP_MPICXX is to force OpenMPI to use the C interface
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${COMMON_FLAGS}")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COMMON_FLAGS}")
 set(CMAKE_C_STANDARD 11)

acc/src/code_generator.c

@@ -606,7 +606,7 @@ generate_headers(void)
 !  -*-f90-*-  (for emacs)    vim:set filetype=fortran:  (for vim)
 
 ! Utils (see astaroth_fortran.cc for definitions)
-external acupdatebuiltinparams
+external achostupdatebuiltinparams
 external acgetdevicecount
 
 ! Device interface (see astaroth_fortran.cc for definitions)

analysis/python/astar/data/read.py

@@ -21,6 +21,7 @@
 
 import numpy as np
 import os
+import pandas as pd
 
 #Optional YT interface
 try:
@@ -302,9 +303,116 @@ class Mesh:
             self.jj = curl_of_curl(self.aa, self.minfo)
         if trim:
             self.bb = (    self.bb[0][3:-3, 3:-3, 3:-3],self.bb[1][3:-3, 3:-3, 3:-3],self.bb[2][3:-3, 3:-3, 3:-3])
+            self.xx_trim = self.xx[3:-3]
+            self.yy_trim = self.yy[3:-3]
+            self.zz_trim = self.zz[3:-3]
             if get_jj:
                 self.jj = (self.jj[0][3:-3, 3:-3, 3:-3],self.jj[1][3:-3, 3:-3, 3:-3],self.jj[2][3:-3, 3:-3, 3:-3])
 
+
+
+    def Bfieldlines(self, footloc = 'default', vartype = 'B', maxstep = 1000):
+        dx = self.minfo.contents['AC_dsx']
+        dy = self.minfo.contents['AC_dsy']
+        dz = self.minfo.contents['AC_dsz']
+ 
+        if vartype == 'U':
+            #Trim to match
+            self.uu = (self.uu[0][3:-3, 3:-3, 3:-3],self.uu[1][3:-3, 3:-3, 3:-3],self.uu[2][3:-3, 3:-3, 3:-3])
+ 
+        def field_line_step(self, coord, ds): 
+            #TODO assume that grid is at a cell centre
+            ix = np.argmin(np.abs(self.xx_trim - coord[0]))            
+            iy = np.argmin(np.abs(self.yy_trim - coord[1]))            
+            iz = np.argmin(np.abs(self.zz_trim - coord[2]))
+            if vartype == 'U':
+                Bcell_vec = np.array([self.uu[0][ix, iy, iz],
+                                      self.uu[1][ix, iy, iz],
+                                      self.uu[2][ix, iy, iz]])
+            else:
+                Bcell_vec = np.array([self.bb[0][ix, iy, iz],
+                                      self.bb[1][ix, iy, iz],
+                                      self.bb[2][ix, iy, iz]])
+ 
+            Bcell_abs = np.sqrt(Bcell_vec[0]**2.0 + Bcell_vec[1]**2.0 + Bcell_vec[2]**2.0) 
+ 
+            coord_new = coord + (Bcell_vec/Bcell_abs)*ds
+            return coord_new
+ 
+        self.df_lines = pd.DataFrame()
+ 
+        ds = np.amin([self.minfo.contents['AC_dsx'], 
+                      self.minfo.contents['AC_dsy'], 
+                      self.minfo.contents['AC_dsz']])
+        ii = 0
+    
+        if footloc == 'middlez':
+            ixtot = 6
+            iytot = 6
+            iztot = 1
+            xfoots = np.linspace(self.xx_trim.min(), self.xx_trim.max(), num = ixtot)
+            yfoots = np.linspace(self.yy_trim.min(), self.yy_trim.max(), num = iytot)
+            zfoots = np.array([(self.zz_trim.max() - self.zz_trim.min())/2.0 + self.zz_trim.min()])
+        elif footloc == 'cube':
+            ixtot = 5
+            iytot = 5
+            iztot = 5
+            xfoots = np.linspace(self.xx_trim.min()+3.0*dx, self.xx_trim.max()-3.0*dx, num = ixtot)
+            yfoots = np.linspace(self.yy_trim.min()+3.0*dy, self.yy_trim.max()-3.0*dy, num = iytot)
+            zfoots = np.linspace(self.zz_trim.min()+3.0*dz, self.zz_trim.max()-3.0*dz, num = iztot)
+        else:
+            ixtot = 6
+            iytot = 6
+            iztot = 1
+            xfoots = np.linspace(self.xx_trim.min(), self.xx_trim.max(), num = ixtot)
+            yfoots = np.linspace(self.yy_trim.min(), self.yy_trim.max(), num = iytot)
+            zfoots = np.array([self.zz_trim.min()])
+ 
+        imax = ixtot * iytot * iztot
+ 
+        for zfoot in zfoots:
+            for yfoot in yfoots:
+                for xfoot in xfoots:
+                    print(ii, "/", imax-1)
+                    integrate = 1 
+                    counter = 0
+                    dstot = 0.0
+                    coord = np.array([xfoot, yfoot, zfoot])
+                    self.df_lines = self.df_lines.append({"line_num":ii, 
+                                                          "dstot":dstot, 
+                                                          "coordx":coord[0], 
+                                                          "coordy":coord[1], 
+                                                          "coordz":coord[2]}, 
+                                                          ignore_index=True) 
+                    while integrate:
+                        coord = field_line_step(self, coord, ds)
+                        dstot += ds
+                        self.df_lines = self.df_lines.append({"line_num":ii, 
+                                                              "dstot":dstot, 
+                                                              "coordx":coord[0], 
+                                                              "coordy":coord[1], 
+                                                              "coordz":coord[2]}, 
+                                                              ignore_index=True) 
+                         
+                        counter += 1
+                        if counter >= maxstep:
+                            integrate = 0   
+                        if ((coord[0] > self.xx_trim.max()) or 
+                            (coord[1] > self.yy_trim.max()) or 
+                            (coord[2] > self.zz_trim.max()) or 
+                            (coord[0] < self.xx_trim.min()) or 
+                            (coord[1] < self.yy_trim.min()) or 
+                            (coord[2] < self.zz_trim.min())): 
+                            #print("out of bounds")
+                            integrate = 0   
+                        if (np.isnan(coord[0]) or 
+                            np.isnan(coord[1]) or 
+                            np.isnan(coord[2])): 
+                            integrate = 0   
+                    ii += 1
+        #print(self.df_lines)
+
+
     def get_jj(self, trim=False):
         self.jj = curl_of_curl(self.aa, minfo, trim=False)
         if trim:

analysis/python/astar/visual/slices.py

@@ -142,4 +142,103 @@ def plot_3(mesh, input_grid, title = '', fname = 'default', bitmap=False,
         print('Saved %s_%s.png' % (fname, mesh.framenum))
         plt.close(fig)
          
+def volume_render(mesh, val1 = {"variable": None, "min": None, "max":None, "opacity":1.0}):
+        
+    if val1["variable"] == "btot":
+        plt.figure()
+        bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
+        array = bb_tot
+        varname = "btot"
+        meshxx = mesh.xx[3:-3]
+        meshyy = mesh.yy[3:-3]
+        meshzz = mesh.zz[3:-3]
+
+    if val1["variable"] == "utot":
+        plt.figure()
+        uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
+        array = uu_tot
+        varname = "utot"
+        meshxx = mesh.xx
+        meshyy = mesh.yy
+        meshzz = mesh.zz
+
+    if val1["variable"] == "rho":
+        plt.figure()
+        array = np.exp(mesh.lnrho)
+        varname = "rho"
+        meshxx = mesh.xx
+        meshyy = mesh.yy
+        meshzz = mesh.zz
+
+    if val1["variable"] == "aa":
+        plt.figure()
+        aa_tot = np.sqrt(mesh.aa[0]**2.0 + mesh.aa[1]**2.0 + mesh.aa[2]**2.0)
+        array = aa_tot
+        varname = "aa"
+        meshxx = mesh.xx
+        meshyy = mesh.yy
+        meshzz = mesh.zz
+
+    #Histogram plot to find value ranges. 
+    hist, bedges = np.histogram(array, bins=mesh.xx.size)
+    plt.plot(bedges[:-1], hist)
+    plt.yscale('log')
+    if val1["min"] != None or val1["max"] != None:
+        plt.plot([val1["min"],val1["min"]], [1,hist.max()], label=varname+" min")
+        plt.plot([val1["max"],val1["max"]], [1,hist.max()], label=varname+" max")
+        plt.legend()
+
+    plt.savefig('volrend_hist_%s_%s.png' % (varname, mesh.framenum))
+    plt.close()
+
+    if val1["min"] != None or val1["max"] != None:
+
+        #print(np.where(bb_tot < val1["min"]))
+    
+        array[np.where(array < val1["min"])] = 0.0
+        array[np.where(array > val1["max"])] = 0.0
+        array[np.where(array > 0.0)] = val1["opacity"]
+    
+        #plt.figure()
+        #plt.plot(bb_tot[:,64,64])
+
+        mapyz = array.sum(axis=0)
+        mapxz = array.sum(axis=1)
+        mapxy = array.sum(axis=2)
+
+        yy_yz, zz_yz = np.meshgrid(meshyy, meshzz, indexing='ij')
+        xx_xz, zz_xz = np.meshgrid(meshxx, meshzz, indexing='ij')
+        xx_xy, yy_xy = np.meshgrid(meshxx, meshyy, indexing='ij')
+
+        fig, ax = plt.subplots()
+        #plt.imshow(mapyz, vmin=0.0, vmax=1.0)
+        plt.pcolormesh(yy_yz, zz_yz, mapyz, vmin=0.0, vmax=1.0, shading='auto')
+        ax.set_aspect('equal')
+        ax.set_title(varname)
+        ax.set_xlabel('y')
+        ax.set_ylabel('z')
+        plt.savefig('volrend_%s_%s_%s.png' % (varname, "yz", mesh.framenum))
+        plt.close()
+    
+        fig, ax = plt.subplots()
+        #plt.imshow(mapxz, vmin=0.0, vmax=1.0)
+        plt.pcolormesh(xx_xz, zz_xz, mapxz, vmin=0.0, vmax=1.0, shading='auto')
+        ax.set_aspect('equal')
+        ax.set_title(varname)
+        ax.set_xlabel('x')
+        ax.set_ylabel('z')
+        plt.savefig('volrend_%s_%s_%s.png' % (varname, "xz", mesh.framenum))
+        plt.close()
+    
+        fig, ax = plt.subplots()
+        #plt.imshow(mapxy, vmin=0.0, vmax=1.0)
+        plt.pcolormesh(xx_xy, yy_xy, mapxy, vmin=0.0, vmax=1.0, shading='auto')
+        ax.set_aspect('equal')
+        ax.set_title(varname)
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        plt.savefig('volrend_%s_%s_%s.png' % (varname, "xy", mesh.framenum))
+        plt.close()
+
+    #plt.show()
  

analysis/python/samples/readtest.py

@@ -25,6 +25,16 @@ import sys
 import os
 import pandas as pd
 
+from mpl_toolkits.mplot3d import Axes3D
+
+#Optional YT interface
+try:
+    import yt
+    yt_present = True 
+except ImportError:
+    yt_present = False
+
+
 ##mesh = ad.read.Mesh(500, fdir="/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/")
 ##
 ##print(np.shape(mesh.uu))
@@ -122,28 +132,110 @@ if "single" in sys.argv:
     print( mesh.uu[2][100, 101, 00], "periodic")
 
 if 'xline' in sys.argv:
-    mesh = ad.read.Mesh(0, fdir=meshdir)
+    mesh_file_numbers = ad.read.parse_directory(meshdir)
-    plt.figure()
+    print(mesh_file_numbers)
-    plt.plot(mesh.uu[0][100, 50, :] , label="z")
+    maxfiles = np.amax(mesh_file_numbers)
-    plt.plot(mesh.uu[0][100, :, 100], label="x")
-    plt.plot(mesh.uu[0][:, 50, 100] , label="y")
-    plt.legend()
 
-    plt.figure()
+    for i in mesh_file_numbers[-3:]:
-    plt.plot(mesh.uu[0][197, 50, :] , label="z edge")
+        mesh = ad.read.Mesh(i, fdir=meshdir)
+        mesh.Bfield(trim=True)
 
-    plt.figure()
+        xhalf = int(mesh.uu[0].shape[0]/2.0)
-    plt.plot(mesh.uu[1][100, 50, :] , label="z")
+        yhalf = int(mesh.uu[0].shape[1]/2.0)
-    plt.plot(mesh.uu[1][100, :, 100], label="x")
+        zhalf = int(mesh.uu[0].shape[2]/2.0)
-    plt.plot(mesh.uu[1][:, 50, 100] , label="y")
+  
-    plt.legend()
+        print(xhalf, yhalf, zhalf)
+
+        plt.figure()
+        plt.plot(mesh.uu[0][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.uu[0][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.uu[0][    :, yhalf, zhalf], label="x")
+        plt.title("UUX")
+        plt.legend()
+
+        #plt.figure()
+        #plt.plot(mesh.uu[0][197, 50, :] , label="z edge")
+
+        plt.figure()
+        plt.plot(mesh.uu[1][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.uu[1][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.uu[1][    :, yhalf, zhalf], label="x")
+        plt.title("UUY")
+        plt.legend()
+
+        plt.figure()
+        plt.plot(mesh.uu[2][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.uu[2][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.uu[2][    :, yhalf, zhalf], label="x")
+        plt.legend()
+        plt.title("UUZ")
+
+        plt.figure()
+        plt.plot(mesh.bb[0][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.bb[0][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.bb[0][    :, yhalf, zhalf], label="x")
+        plt.title("BBX")
+        plt.legend()
+
+        plt.figure()
+        plt.plot(mesh.bb[1][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.bb[1][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.bb[1][    :, yhalf, zhalf], label="x")
+        plt.title("BBY")
+        plt.legend()
+
+        plt.figure()
+        plt.plot(mesh.bb[2][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.bb[2][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.bb[2][    :, yhalf, zhalf], label="x")
+        plt.legend()
+        plt.title("BBZ")
+
+        plt.figure()
+        plt.plot(mesh.aa[0][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.aa[0][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.aa[0][    :, yhalf, zhalf], label="x")
+        plt.title("AX")
+        plt.legend()
+
+        plt.figure()
+        plt.plot(mesh.aa[1][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.aa[1][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.aa[1][    :, yhalf, zhalf], label="x")
+        plt.title("AY")
+        plt.legend()
+
+        plt.figure()
+        plt.plot(mesh.aa[2][xhalf, yhalf,     :], label="z")
+        plt.plot(mesh.aa[2][xhalf,     :, zhalf], label="y")
+        plt.plot(mesh.aa[2][    :, yhalf, zhalf], label="x")
+        plt.legend()
+        plt.title("AZ")
+
+        uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
+        bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
+
+        plt.figure()
+        plt.plot(uu_tot[xhalf, yhalf,     :], label="z")
+        plt.plot(uu_tot[xhalf,     :, zhalf], label="y")
+        plt.plot(uu_tot[    :, yhalf, zhalf], label="x")
+        plt.legend()
+        plt.title("UTOT")
+
+        plt.figure()
+        plt.plot(bb_tot[xhalf, yhalf,     :], label="z")
+        plt.plot(bb_tot[xhalf,     :, zhalf], label="y")
+        plt.plot(bb_tot[    :, yhalf, zhalf], label="x")
+        plt.legend()
+        plt.title("BTOT")
+
+        plt.figure()
+        plt.plot(np.exp(mesh.lnrho[xhalf, yhalf,     :]), label="z")
+        plt.plot(np.exp(mesh.lnrho[xhalf,     :, zhalf]), label="y")
+        plt.plot(np.exp(mesh.lnrho[    :, yhalf, zhalf]), label="x")
+        plt.legend()
+        plt.title("RHO")
 
-    plt.figure()
-    plt.plot(mesh.uu[2][100, 50, :] , label="z")
-    plt.plot(mesh.uu[2][100, :, 100], label="x")
-    plt.plot(mesh.uu[2][:, 50, 100] , label="y")
-    plt.legend()
-    plt.show()
 
 if 'check' in sys.argv:
     mesh = ad.read.Mesh(0, fdir=meshdir)
@@ -180,22 +272,43 @@ if '1d' in sys.argv:
  
             plt.legend()
 
+
     plt.show()
 
+if 'csv' in sys.argv:
+    filenum = sys.argv[1]
+    mesh = ad.read.Mesh(filenum, fdir=meshdir)
+    mesh.Bfield()
+    mesh.export_csv()
+
+if 'raw' in sys.argv:
+    filenum = sys.argv[1]
+    mesh = ad.read.Mesh(filenum, fdir=meshdir)
+    mesh.Bfield()
+    mesh.export_raw()
+
+if 'findnan' in sys.argv:
+    filenum = sys.argv[1]
+    mesh = ad.read.Mesh(filenum, fdir=meshdir)
+    print("nan uu", np.where(np.isnan(mesh.uu)))
+    print("nan aa", np.where(np.isnan(mesh.aa)))
+    print("nan lnrho", np.where(np.isnan(mesh.lnrho)))
+    print("inf uu", np.where(np.isinf(mesh.uu)))
+    print("inf aa", np.where(np.isinf(mesh.aa)))
+    print("inf lnrho", np.where(np.isinf(mesh.lnrho)))
+
 
 if 'sl' in sys.argv:
-    #maxfiles = 200002
+    mesh_file_numbers = ad.read.parse_directory(meshdir)
-    #stride = 10000
+    print(mesh_file_numbers)
-    maxfiles = 500000001
+    maxfiles = np.amax(mesh_file_numbers)
-    stride = 1
+    for i in mesh_file_numbers:
-    for i in range(0, maxfiles, stride):
-        #mesh = ad.read.Mesh(i, fdir=meshdir) 
         mesh = ad.read.Mesh(i, fdir=meshdir) 
         print(" %i / %i" % (i, maxfiles))
         if mesh.ok:
             uu_tot = np.sqrt(mesh.uu[0]**2.0 + mesh.uu[1]**2.0 + mesh.uu[2]**2.0)
             aa_tot = np.sqrt(mesh.aa[0]**2.0 + mesh.aa[1]**2.0 + mesh.aa[2]**2.0)
-            mesh.Bfield()
+            mesh.Bfield(trim=True)
             bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0)
 
             if 'sym' in sys.argv:
@@ -248,6 +361,133 @@ if 'sl' in sys.argv:
                 vis.slices.plot_3(mesh, mesh.bb[1],         title = r'$B_y$', bitmap = True, fname = 'bby',     trimghost=3)#, bfieldlines=True)
                 vis.slices.plot_3(mesh, mesh.bb[2],         title = r'$B_z$', bitmap = True, fname = 'bbz',     trimghost=3)#, bfieldlines=True)
                  
+            if 'yt' in sys.argv:
+                mesh.yt_conversion()
+                from mpl_toolkits.axes_grid1 import AxesGrid
+
+                coords = ['x', 'y','z']
+                for coord in coords:
+                    fields = ['density', 'uux', 'uuy', 'uuz']
+                    fig = plt.figure()
+                    grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
+                                    nrows_ncols = (2, 2),
+                                    axes_pad = 1.0,
+                                    label_mode = "1",
+                                    share_all = True,
+                                    cbar_location="right",
+                                    cbar_mode="each",
+                                    cbar_size="3%",
+                                    cbar_pad="0%")
+
+                    p = yt.SlicePlot(mesh.ytdata, coord, fields)
+                    p.set_log('uux', False)
+                    p.set_log('uuy', False)
+                    p.set_log('uuz', False)
+
+                    for i, field in enumerate(fields):
+                        plot = p.plots[field]
+                        plot.figure = fig
+                        plot.axes = grid[i].axes
+                        plot.cax = grid.cbar_axes[i]
+
+                    p._setup_plots()
+                    plt.savefig('yt_rho_uu_%s_%s.png' % (coord, mesh.framenum))
+
+                    plt.close(fig=fig)
+
+                    ###
+
+                    fields = ['density', 'bbx', 'bby', 'bbz']
+                    fig = plt.figure()
+                    grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
+                                    nrows_ncols = (2, 2),
+                                    axes_pad = 1.0,
+                                    label_mode = "1",
+                                    share_all = True,
+                                    cbar_location="right",
+                                    cbar_mode="each",
+                                    cbar_size="3%",
+                                    cbar_pad="0%")
+
+                    p = yt.SlicePlot(mesh.ytdata, coord, fields)
+                    p.set_log('bbx', False)
+                    p.set_log('bby', False)
+                    p.set_log('bbz', False)
+
+                    for i, field in enumerate(fields):
+                        plot = p.plots[field]
+                        plot.figure = fig
+                        plot.axes = grid[i].axes
+                        plot.cax = grid.cbar_axes[i]
+
+                    p._setup_plots()
+                    plt.savefig('yt_rho_bb_%s_%s.png' % (coord, mesh.framenum))
+
+                    plt.close(fig=fig)
+
+            elif 'csvall' in sys.argv:
+                mesh.export_csv()
+
+            elif 'rawall' in sys.argv:
+                mesh.export_raw()
+   
+if "vol" in sys.argv: 
+   print("VOLUME RENDERING") 
+   mesh_file_numbers = ad.read.parse_directory(meshdir)
+   print(mesh_file_numbers)
+   maxfiles = np.amax(mesh_file_numbers)
+
+   for i in mesh_file_numbers:
+       mesh = ad.read.Mesh(i, fdir=meshdir) 
+       mesh.Bfield(trim=True)
+       print(" %i / %i" % (i, maxfiles))
+       if mesh.ok:
+           vis.slices.volume_render(mesh, val1 = {"variable": "btot", "min":0.5, "max":2.0, "opacity":0.05})
+           vis.slices.volume_render(mesh, val1 = {"variable": "utot", "min":0.5, "max":2.0, "opacity":0.05})
+           vis.slices.volume_render(mesh, val1 = {"variable": "rho", "min":10.0, "max":300.0, "opacity":0.05})
+           vis.slices.volume_render(mesh, val1 = {"variable": "aa", "min":0.1, "max":0.25, "opacity":0.05})
+
+if (("bline" in sys.argv) or ("uline" in sys.argv)): 
+   print("Field line computation") 
+   mesh_file_numbers = ad.read.parse_directory(meshdir)
+   print(mesh_file_numbers)
+   maxfiles = np.amax(mesh_file_numbers)
+
+   for i in mesh_file_numbers:
+       mesh = ad.read.Mesh(i, fdir=meshdir) 
+       mesh.Bfield(trim=True)
+       print(" %i / %i" % (i, maxfiles))
+       if mesh.ok:
+           if "uline" in sys.argv:
+               mesh.Bfieldlines(footloc = 'cube', vartype='U', maxstep = 200)
+           else: 
+               mesh.Bfieldlines(footloc = 'default')
+           print(mesh.df_lines)
+    
+       fig = plt.figure(figsize=(5.0,5.0))
+       ax = fig.gca(projection='3d')
+       for line_num in range(int(mesh.df_lines['line_num'].max()+1)):
+           df_myline = mesh.df_lines.loc[mesh.df_lines['line_num'] == line_num]
+           print(df_myline)
+           my_xscale = [mesh.xx_trim.min(), mesh.xx_trim.max()]
+           my_yscale = [mesh.yy_trim.min(), mesh.yy_trim.max()]
+           my_zscale = [mesh.zz_trim.min(), mesh.zz_trim.max()]
+           ax.plot(df_myline["coordx"], df_myline["coordy"], df_myline["coordz"], color="red")
+           ax.set_xlim3d(my_xscale)
+           ax.set_ylim3d(my_yscale)
+           ax.set_zlim3d(my_zscale)
+
+       if "uline" in sys.argv:
+           filename = 'Ugeometry_%s.png' % (mesh.framenum)
+       else:
+           filename = 'Bgeometry_%s.png' % (mesh.framenum)  
+
+       plt.savefig(filename)
+       plt.close()
+
+       
+   #plt.show()
+
 
 if 'ts' in sys.argv:
    ts = ad.read.TimeSeries(fdir=meshdir)

config/astaroth.conf

@@ -5,9 +5,9 @@
  * "Compile-time" params
  * =============================================================================
  */
-AC_nx = 128
+AC_nx = 64
-AC_ny = 128
+AC_ny = 64
-AC_nz = 128
+AC_nz = 64
 
 AC_dsx = 0.04908738521
 AC_dsy = 0.04908738521

include/astaroth.h

@@ -46,13 +46,16 @@ typedef struct {
 } AcMatrix;
 
 #include "user_defines.h" // Autogenerated defines from the DSL
+//#include "../src/core/kernels/kernels.h"
 
 typedef enum { AC_SUCCESS = 0, AC_FAILURE = 1 } AcResult;
 
-// Neming the associated number of the boundary condition types
+// Naming the associated number of the boundary condition types
-typedef enum { AC_BOUNDCOND_PERIODIC = 0, 
+typedef enum {
-               AC_BOUNDCOND_SYMMETRIC = 1, 
+    AC_BOUNDCOND_PERIODIC      = 0,
-               AC_BOUNDCOND_ANTISYMMETRIC = 2 } AcBoundcond;
+    AC_BOUNDCOND_SYMMETRIC     = 1,
+    AC_BOUNDCOND_ANTISYMMETRIC = 2
+} AcBoundcond;
 
 #define AC_GEN_ID(X) X,
 typedef enum {
@@ -225,6 +228,9 @@ AcResult acLoadDeviceConstant(const AcRealParam param, const AcReal value);
 /** Loads an AcMesh to the devices visible to the caller */
 AcResult acLoad(const AcMesh host_mesh);
 
+/** Sets the whole mesh to some value */
+AcResult acSetVertexBuffer(const VertexBufferHandle handle, const AcReal value);
+
 /** Stores the AcMesh distributed among the devices visible to the caller back to the host*/
 AcResult acStore(AcMesh* host_mesh);
 
@@ -252,10 +258,11 @@ AcReal acReduceScal(const ReductionType rtype, const VertexBufferHandle vtxbuf_h
 AcReal acReduceVec(const ReductionType rtype, const VertexBufferHandle a,
                    const VertexBufferHandle b, const VertexBufferHandle c);
 
-/** Does a reduction for an operation which requires a vector and a scalar with vertex buffers 
+/** Does a reduction for an operation which requires a vector and a scalar with vertex buffers
  *  * where the vector components are (a, b, c) and scalr is (d) */
 AcReal acReduceVecScal(const ReductionType rtype, const VertexBufferHandle a,
-                       const VertexBufferHandle b, const VertexBufferHandle c, const VertexBufferHandle d);
+                       const VertexBufferHandle b, const VertexBufferHandle c,
+                       const VertexBufferHandle d);
 
 /** Stores a subset of the mesh stored across the devices visible to the caller back to host memory.
  */
@@ -322,7 +329,7 @@ AcResult acGridIntegrate(const Stream stream, const AcReal dt);
 
 /** */
 /*   MV: Commented out for a while, but save for the future when standalone_MPI
-         works with periodic boundary conditions. 
+         works with periodic boundary conditions.
 AcResult
 acGridIntegrateNonperiodic(const Stream stream, const AcReal dt)
 
@@ -432,6 +439,10 @@ AcResult acNodeLoadVertexBuffer(const Node node, const Stream stream, const AcMe
 /** */
 AcResult acNodeLoadMesh(const Node node, const Stream stream, const AcMesh host_mesh);
 
+/** */
+AcResult acNodeSetVertexBuffer(const Node node, const Stream stream,
+                               const VertexBufferHandle handle, const AcReal value);
+
 /** Deprecated ? */
 AcResult acNodeStoreVertexBufferWithOffset(const Node node, const Stream stream,
                                            const VertexBufferHandle vtxbuf_handle, const int3 src,
@@ -467,8 +478,9 @@ AcResult acNodePeriodicBoundcondStep(const Node node, const Stream stream,
 AcResult acNodePeriodicBoundconds(const Node node, const Stream stream);
 
 /** */
-AcResult acNodeGeneralBoundcondStep(const Node node, const Stream stream,   
+AcResult acNodeGeneralBoundcondStep(const Node node, const Stream stream,
-                                    const VertexBufferHandle vtxbuf_handle, const AcMeshInfo config);
+                                    const VertexBufferHandle vtxbuf_handle,
+                                    const AcMeshInfo config);
 
 /** */
 AcResult acNodeGeneralBoundconds(const Node node, const Stream stream, const AcMeshInfo config);
@@ -483,8 +495,8 @@ AcResult acNodeReduceVec(const Node node, const Stream stream_type, const Reduct
 /** */
 AcResult acNodeReduceVecScal(const Node node, const Stream stream_type, const ReductionType rtype,
                              const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1,
-                             const VertexBufferHandle vtxbuf2, const VertexBufferHandle vtxbuf3, AcReal* result);
+                             const VertexBufferHandle vtxbuf2, const VertexBufferHandle vtxbuf3,
-
+                             AcReal* result);
 
 /*
  * =============================================================================
@@ -554,6 +566,10 @@ AcResult acDeviceLoadVertexBuffer(const Device device, const Stream stream, cons
 /** */
 AcResult acDeviceLoadMesh(const Device device, const Stream stream, const AcMesh host_mesh);
 
+/** */
+AcResult acDeviceSetVertexBuffer(const Device device, const Stream stream,
+                                 const VertexBufferHandle handle, const AcReal value);
+
 /** */
 AcResult acDeviceStoreVertexBufferWithOffset(const Device device, const Stream stream,
                                              const VertexBufferHandle vtxbuf_handle, const int3 src,
@@ -610,7 +626,6 @@ AcResult acDeviceGeneralBoundcondStep(const Device device, const Stream stream,
 AcResult acDeviceGeneralBoundconds(const Device device, const Stream stream, const int3 start,
                                    const int3 end, const AcMeshInfo config, const int3 bindex);
 
-
 /** */
 AcResult acDeviceReduceScal(const Device device, const Stream stream, const ReductionType rtype,
                             const VertexBufferHandle vtxbuf_handle, AcReal* result);
@@ -619,9 +634,10 @@ AcResult acDeviceReduceVec(const Device device, const Stream stream_type, const
                            const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1,
                            const VertexBufferHandle vtxbuf2, AcReal* result);
 /** */
-AcResult acDeviceReduceVecScal(const Device device, const Stream stream_type, const ReductionType rtype,
+AcResult acDeviceReduceVecScal(const Device device, const Stream stream_type,
-                               const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1,
+                               const ReductionType rtype, const VertexBufferHandle vtxbuf0,
-                               const VertexBufferHandle vtxbuf2, const VertexBufferHandle vtxbuf3, AcReal* result);
+                               const VertexBufferHandle vtxbuf1, const VertexBufferHandle vtxbuf2,
+                               const VertexBufferHandle vtxbuf3, AcReal* result);
 /** */
 AcResult acDeviceRunMPITest(void);
 
@@ -631,16 +647,16 @@ AcResult acDeviceRunMPITest(void);
  * =============================================================================
  */
 /** Updates the built-in parameters based on nx, ny and nz */
-AcResult acUpdateBuiltinParams(AcMeshInfo* config);
+AcResult acHostUpdateBuiltinParams(AcMeshInfo* config);
 
 /** Creates a mesh stored in host memory */
-AcResult acMeshCreate(const AcMeshInfo mesh_info, AcMesh* mesh);
+AcResult acHostMeshCreate(const AcMeshInfo mesh_info, AcMesh* mesh);
 
 /** Randomizes a host mesh */
-AcResult acMeshRandomize(AcMesh* mesh);
+AcResult acHostMeshRandomize(AcMesh* mesh);
 
 /** Destroys a mesh stored in host memory */
-AcResult acMeshDestroy(AcMesh* mesh);
+AcResult acHostMeshDestroy(AcMesh* mesh);
 
 #ifdef __cplusplus
 } // extern "C"

include/astaroth_utils.h

@@ -45,25 +45,25 @@ typedef struct {
 AcResult acLoadConfig(const char* config_path, AcMeshInfo* config);
 
 /** */
-AcResult acVertexBufferSet(const VertexBufferHandle handle, const AcReal value, AcMesh* mesh);
+AcResult acHostVertexBufferSet(const VertexBufferHandle handle, const AcReal value, AcMesh* mesh);
 
 /** */
-AcResult acMeshSet(const AcReal value, AcMesh* mesh);
+AcResult acHostMeshSet(const AcReal value, AcMesh* mesh);
 
 /** */
-AcResult acMeshApplyPeriodicBounds(AcMesh* mesh);
+AcResult acHostMeshApplyPeriodicBounds(AcMesh* mesh);
 
 /** */
-AcResult acMeshClear(AcMesh* mesh);
+AcResult acHostMeshClear(AcMesh* mesh);
 
 /** */
-AcResult acModelIntegrateStep(AcMesh mesh, const AcReal dt);
+AcResult acHostIntegrateStep(AcMesh mesh, const AcReal dt);
 
 /** */
-AcReal acModelReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a);
+AcReal acHostReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a);
 
 /** */
-AcReal acModelReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a,
+AcReal acHostReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a,
                         const VertexBufferHandle b, const VertexBufferHandle c);
 
 Error acGetError(const AcReal model, const AcReal candidate);

samples/benchmark/main.cc

@@ -98,7 +98,7 @@ main(int argc, char** argv)
             info.int_params[AC_nx] = nx;
             info.int_params[AC_ny] = ny;
             info.int_params[AC_nz] = nz;
-            acUpdateBuiltinParams(&info);
+            acHostUpdateBuiltinParams(&info);
             printf("Benchmark mesh dimensions: (%d, %d, %d)\n", nx, ny, nz);
         }
         else {
@@ -107,7 +107,7 @@ main(int argc, char** argv)
         }
     }
 
-    const TestType test = TEST_WEAK_SCALING;
+    const TestType test = TEST_STRONG_SCALING;
     if (test == TEST_WEAK_SCALING) {
         uint3_64 decomp = decompose(nprocs);
         info.int_params[AC_nx] *= decomp.x;
@@ -118,10 +118,10 @@ main(int argc, char** argv)
     /*
     AcMesh model, candidate;
     if (pid == 0) {
-        acMeshCreate(info, &model);
+        acHostMeshCreate(info, &model);
-        acMeshCreate(info, &candidate);
+        acHostMeshCreate(info, &candidate);
-        acMeshRandomize(&model);
+        acHostMeshRandomize(&model);
-        acMeshRandomize(&candidate);
+        acHostMeshRandomize(&candidate);
     }*/
 
     // GPU alloc & compute
@@ -130,8 +130,8 @@ main(int argc, char** argv)
 
     /*
     AcMesh model;
-    acMeshCreate(info, &model);
+    acHostMeshCreate(info, &model);
-    acMeshRandomize(&model);
+    acHostMeshRandomize(&model);
     acGridLoadMesh(STREAM_DEFAULT, model);
     */
 
@@ -145,12 +145,12 @@ main(int argc, char** argv)
 
     // Verify
     if (pid == 0) {
-        acModelIntegrateStep(model, FLT_EPSILON);
+        acHostIntegrateStep(model, FLT_EPSILON);
-        acMeshApplyPeriodicBounds(&model);
+        acHostMeshApplyPeriodicBounds(&model);
 
         AcResult retval = acVerifyMesh(model, candidate);
-        acMeshDestroy(&model);
+        acHostMeshDestroy(&model);
-        acMeshDestroy(&candidate);
+        acHostMeshDestroy(&candidate);
 
         if (retval != AC_SUCCESS) {
             fprintf(stderr, "Failures found, benchmark invalid. Skipping\n");

samples/cpptest/main.cc

@@ -31,12 +31,12 @@ main(void)
 
     // Alloc
     AcMesh model, candidate;
-    acMeshCreate(info, &model);
+    acHostMeshCreate(info, &model);
-    acMeshCreate(info, &candidate);
+    acHostMeshCreate(info, &candidate);
 
     // Init
-    acMeshRandomize(&model);
+    acHostMeshRandomize(&model);
-    acMeshApplyPeriodicBounds(&model);
+    acHostMeshApplyPeriodicBounds(&model);
 
     // Verify that the mesh was loaded and stored correctly
     acInit(info);
@@ -55,8 +55,8 @@ main(void)
 
     // Destroy
     acQuit();
-    acMeshDestroy(&model);
+    acHostMeshDestroy(&model);
-    acMeshDestroy(&candidate);
+    acHostMeshDestroy(&candidate);
 
     puts("cpptest complete.");
     return EXIT_SUCCESS;

samples/ctest/main.c

@@ -30,12 +30,12 @@ main(void)
 
     // Alloc
     AcMesh model, candidate;
-    acMeshCreate(info, &model);
+    acHostMeshCreate(info, &model);
-    acMeshCreate(info, &candidate);
+    acHostMeshCreate(info, &candidate);
 
     // Init
-    acMeshRandomize(&model);
+    acHostMeshRandomize(&model);
-    acMeshApplyPeriodicBounds(&model);
+    acHostMeshApplyPeriodicBounds(&model);
 
     // Verify that the mesh was loaded and stored correctly
     acInit(info);
@@ -46,6 +46,7 @@ main(void)
     // Attempt to integrate and check max and min
     printf("Integrating... ");
     acIntegrate(FLT_EPSILON);
+
     printf("Done.\nVTXBUF ranges after one integration step:\n");
     for (size_t i = 0; i < NUM_VTXBUF_HANDLES; ++i)
         printf("\t%-15s... [%.3g, %.3g]\n", vtxbuf_names[i], //
@@ -54,8 +55,8 @@ main(void)
 
     // Destroy
     acQuit();
-    acMeshDestroy(&model);
+    acHostMeshDestroy(&model);
-    acMeshDestroy(&candidate);
+    acHostMeshDestroy(&candidate);
 
     puts("ctest complete.");
     return EXIT_SUCCESS;

samples/fortrantest/main.f90

@@ -14,7 +14,7 @@ program pc
   info%int_params(AC_nx + 1) = 128
   info%int_params(AC_ny + 1) = 128
   info%int_params(AC_nz + 1) = 128
-  call acupdatebuiltinparams(info)
+  call achostupdatebuiltinparams(info)
 
   call acdevicecreate(0, info, device)
   call acdeviceprintinfo(device)

samples/mpitest/main.cc

@@ -23,6 +23,271 @@
 #include "astaroth_utils.h"
 #include "errchk.h"
 
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+using namespace std;
+
+void
+gaussian_radial_explosion(AcMesh* mesh)
+{
+    AcReal* uu_x     = mesh->vertex_buffer[VTXBUF_UUX];
+    AcReal* uu_y     = mesh->vertex_buffer[VTXBUF_UUY];
+    AcReal* uu_z     = mesh->vertex_buffer[VTXBUF_UUZ];
+    const int mx     = mesh->info.int_params[AC_mx];
+    const int my     = mesh->info.int_params[AC_my];
+    const int nx_min = mesh->info.int_params[AC_nx_min];
+    const int nx_max = mesh->info.int_params[AC_nx_max];
+    const int ny_min = mesh->info.int_params[AC_nx_min];
+    const int ny_max = mesh->info.int_params[AC_nx_max];
+    const int nz_min = mesh->info.int_params[AC_nx_min];
+    const int nz_max = mesh->info.int_params[AC_nx_max];
+    const AcReal DX  = mesh->info.real_params[AC_dsx];
+    const AcReal DY  = mesh->info.real_params[AC_dsy];
+    const AcReal DZ  = mesh->info.real_params[AC_dsz];
+
+
+    const AcReal xorig = 0.00001;
+    const AcReal yorig = 0.00001;
+    const AcReal zorig = 0.00001;
+
+    const AcReal INIT_LOC_UU_X = 0.01;
+    const AcReal INIT_LOC_UU_Y = 32 * DY;
+    const AcReal INIT_LOC_UU_Z = 32 * DZ;
+    const AcReal AMPL_UU       = 1;
+    const AcReal UU_SHELL_R    = 0.8;
+    const AcReal WIDTH_UU      = 0.2;
+    // Outward explosion with gaussian initial velocity profile.
+    int idx;
+    AcReal xx, yy, zz, rr2, rr, theta = 0.0, phi = 0.0;
+    AcReal uu_radial;
+    // real theta_old = 0.0;
+    for (int k = nz_min; k < nz_max; k++) {
+        for (int j = ny_min; j < ny_max; j++) {
+            for (int i = nx_min; i < nx_max; i++) {
+                // Calculate the value of velocity in a particular radius.
+                idx = i + j * mx + k * mx * my;
+                // Determine the coordinates
+                xx  = DX * (i - nx_min) - xorig;
+                xx  = xx - INIT_LOC_UU_X;
+                yy  = DY * (j - ny_min) - yorig;
+                yy  = yy - INIT_LOC_UU_Y;
+                zz  = DZ * (k - nz_min) - zorig;
+                zz  = zz - INIT_LOC_UU_Z;
+                rr2 = pow(xx, 2.0) + pow(yy, 2.0) + pow(zz, 2.0);
+                rr  = sqrt(rr2);
+
+                // printf("[%d %d %d] %e %e %e\n", i, j, k, DX, DY, DZ);
+
+                // Origin is different!
+                AcReal xx_abs, yy_abs, zz_abs;
+                if (rr > 0.0) {
+                    // theta range [0, PI]
+                    if (zz >= 0.0) {
+                        theta = acos(min(1.0, zz / rr));
+                        if (theta > M_PI / 2.0 || theta < 0.0) {
+                            printf("Explosion THETA WRONG: zz = %.3f, rr = %.3f, theta = %.3e/PI, "
+                                   "M_PI = %.3e\n",
+                                   zz, rr, theta / M_PI, M_PI);
+                        }
+                    }
+                    else {
+                        zz_abs = -zz; // Needs a posite value for acos
+                        theta  = M_PI - acos(zz_abs / rr);
+                        if (theta < M_PI / 2.0 || theta > 2 * M_PI) {
+                            printf("Explosion THETA WRONG: zz = %.3f, rr = %.3f, theta = %.3e/PI, "
+                                   "M_PI = %.3e\n",
+                                   zz, rr, theta / M_PI, M_PI);
+                        }
+                    }
+                    // phi range [0, 2*PI]i
+                    if (xx != 0.0) {
+                        if (xx < 0.0 && yy >= 0.0) {
+                            //-+
+                            xx_abs = -xx; // Needs a posite value for atan
+                            phi    = M_PI - atan(yy / xx_abs);
+                            if (phi < (M_PI / 2.0) || phi > M_PI) {
+                                printf("Explosion PHI WRONG -+: xx = %.3f, yy = %.3f, phi = "
+                                       "%.3e/PI, M_PI = %.3e\n",
+                                       xx, yy, phi / M_PI, M_PI);
+                            }
+                        }
+                        else if (xx > 0.0 && yy < 0.0) {
+                            //+-
+                            yy_abs = -yy;
+                            phi    = 2.0 * M_PI - atan(yy_abs / xx);
+                            if (phi < (3.0 * M_PI) / 2.0 || phi > (2.0 * M_PI + 1e-6)) {
+                                printf("Explosion PHI WRONG +-: xx = %.3f, yy = %.3f, phi = "
+                                       "%.3e/PI, M_PI = %.3e\n",
+                                       xx, yy, phi / M_PI, M_PI);
+                            }
+                        }
+                        else if (xx < 0.0 && yy < 0.0) {
+                            //--
+                            yy_abs = -yy;
+                            xx_abs = -xx;
+                            phi    = M_PI + atan(yy_abs / xx_abs);
+                            if (phi < M_PI || phi > ((3.0 * M_PI) / 2.0 + 1e-6)) {
+                                printf("Explosion PHI WRONG --: xx = %.3f, yy = %.3f, xx_abs = "
+                                       "%.3f, yy_abs = %.3f, phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
+                                       xx, yy, xx_abs, yy_abs, phi, (3.0 * M_PI) / 2.0);
+                            }
+                        }
+                        else {
+                            //++
+                            phi = atan(yy / xx);
+                            if (phi < 0 || phi > M_PI / 2.0) {
+                                printf("Explosion PHI WRONG --: xx = %.3f, yy = %.3f, phi = %.3e, "
+                                       "(3.0*M_PI)/2.0 = %.3e\n",
+                                       xx, yy, phi, (3.0 * M_PI) / 2.0);
+                            }
+                        }
+                    }
+                    else { // To avoid div by zero with atan
+                        if (yy > 0.0) {
+                            phi = M_PI / 2.0;
+                        }
+                        else if (yy < 0.0) {
+                            phi = (3.0 * M_PI) / 2.0;
+                        }
+                        else {
+                            phi = 0.0;
+                        }
+                    }
+                    // Set zero for explicit safekeeping
+                    if (xx == 0.0 && yy == 0.0) {
+                        phi = 0.0;
+                    }
+                    // Gaussian velocity
+                    // uu_radial = AMPL_UU*exp( -rr2 / (2.0*pow(WIDTH_UU, 2.0)) );
+                    // New distribution, where that gaussion wave is not in the exact centre
+                    // coordinates uu_radial = AMPL_UU*exp( -pow((rr - 4.0*WIDTH_UU),2.0) /
+                    // (2.0*pow(WIDTH_UU, 2.0)) ); //TODO: Parametrize the peak location.
+                    uu_radial = AMPL_UU *
+                                exp(-pow((rr - UU_SHELL_R), 2.0) / (2.0 * pow(WIDTH_UU, 2.0)));
+                }
+                else {
+                    uu_radial = 0.0; // TODO: There will be a discontinuity in the origin... Should
+                                     // the shape of the distribution be different?
+                }
+
+                // if (rr > 0.2 && rr < 1.0) {
+                //     printf("%e\n", uu_radial);
+                // }
+
+                // Determine the carthesian velocity components and lnrho
+                uu_x[idx] = uu_radial * sin(theta) * cos(phi);
+                uu_y[idx] = uu_radial * sin(theta) * sin(phi);
+                uu_z[idx] = uu_radial * cos(theta);
+                // Temporary diagnosticv output (TODO: Remove after not needed)
+                // if (theta > theta_old) {
+                // if (theta > M_PI || theta < 0.0 || phi < 0.0 || phi > 2*M_PI) {
+                /*  printf("Explosion: xx = %.3f, yy = %.3f, zz = %.3f, rr = %.3f, phi = %.3e/PI,
+                   theta = %.3e/PI\n, M_PI = %.3e", xx, yy, zz, rr, phi/M_PI, theta/M_PI, M_PI);
+                    printf("           uu_radial = %.3e, uu_x[%i] = %.3e, uu_y[%i] = %.3e, uu_z[%i]
+                   = %.3e \n", uu_radial, idx, uu_x[idx], idx, uu_y[idx], idx, uu_z[idx]); theta_old
+                   = theta;
+                */
+            }
+        }
+    }
+}
+
+static AcReal
+randf(void)
+{
+    return (AcReal)rand() / (AcReal)RAND_MAX;
+}
+
+AcResult
+meshRadial(AcMesh* mesh)
+{
+
+    const int n = acVertexBufferSize(mesh->info);
+
+    // lnrho to constant
+    for (int i = 0; i < n; ++i)
+        mesh->vertex_buffer[VTXBUF_LNRHO][i] = 0.5;
+
+    // A to random
+    for (int i = 0; i < n; ++i)
+        mesh->vertex_buffer[VTXBUF_AX][i] = randf();
+    for (int i = 0; i < n; ++i)
+        mesh->vertex_buffer[VTXBUF_AY][i] = randf();
+    for (int i = 0; i < n; ++i)
+        mesh->vertex_buffer[VTXBUF_AZ][i] = randf();
+
+    // entropy random
+    for (int i = 0; i < n; ++i)
+        mesh->vertex_buffer[VTXBUF_ENTROPY][i] = randf();
+
+    // velocity is radial explosion
+    gaussian_radial_explosion(mesh);
+
+    return AC_SUCCESS;
+}
+
+void
+dumpMesh(AcMesh* mesh, const std::string& path)
+{
+    const char delim[] = ",";
+
+    FILE* outf = fopen(path.c_str(), "w");
+    if (!outf) {
+        std::cerr << "unable to open \"" << path << "\" for writing\n";
+        exit(1);
+    }
+
+    // column headers
+    fprintf(outf, "Z%sY%sX", delim, delim);
+    for (int64_t qi = 0; qi < NUM_VTXBUF_HANDLES; ++qi) {
+        std::string colName = "data" + std::to_string(qi);
+        fprintf(outf, "%s%s", delim, colName.c_str());
+    }
+    fprintf(outf, "\n");
+
+    // is this right?
+    int3 origin = mesh->info.int3_params[AC_multigpu_offset];
+    origin.x += 1;
+    origin.y += 1;
+    origin.z += 1;
+
+    const int nz = mesh->info.int_params[AC_nz];
+    const int ny = mesh->info.int_params[AC_ny];
+    const int nx = mesh->info.int_params[AC_nx];
+    // const int mz = mesh->info.int_params[AC_mz];
+    const int my = mesh->info.int_params[AC_my];
+    const int mx = mesh->info.int_params[AC_mx];
+
+    // rows
+    for (int lz = 0; lz < nz; ++lz) {
+        for (int ly = 0; ly < ny; ++ly) {
+            for (int lx = 0; lx < nx; ++lx) {
+
+                const int3 pos{origin.x + lx, origin.y + ly, origin.z + lz};
+
+                fprintf(outf, "%d%s%d%s%d", pos.z, delim, pos.y, delim, pos.x);
+
+                for (int64_t qi = 0; qi < NUM_VTXBUF_HANDLES; ++qi) {
+                    AcReal val = mesh->vertex_buffer[qi][lz * (my * mx) + ly * mx + lx];
+                    if (8 == sizeof(AcReal)) {
+                        fprintf(outf, "%s%.17f", delim, val);
+                    }
+                    else if (4 == sizeof(AcReal)) {
+                        fprintf(outf, "%s%.9f", delim, val);
+                    }
+                }
+
+                fprintf(outf, "\n");
+            }
+        }
+    }
+
+    fclose(outf);
+}
+
 #if AC_MPI_ENABLED
 
 #include <mpi.h>
@@ -47,10 +312,18 @@ main(void)
 
     AcMesh model, candidate;
     if (pid == 0) {
-        acMeshCreate(info, &model);
+        acHostMeshCreate(info, &model);
-        acMeshCreate(info, &candidate);
+        acHostMeshCreate(info, &candidate);
-        acMeshRandomize(&model);
+        meshRadial(&model);
-        acMeshRandomize(&candidate);
+        meshRadial(&candidate);
+    }
+
+    {
+        std::stringstream ss;
+        ss << "candidate_init_";
+        ss << pid << ".txt";
+        std::cerr << "dump to " << ss.str() << "\n";
+        dumpMesh(&candidate, ss.str());
     }
 
     // GPU alloc & compute
@@ -61,10 +334,10 @@ main(void)
     acGridPeriodicBoundconds(STREAM_DEFAULT);
     acGridStoreMesh(STREAM_DEFAULT, &candidate);
     if (pid == 0) {
-        acMeshApplyPeriodicBounds(&model);
+        acHostMeshApplyPeriodicBounds(&model);
         const AcResult res = acVerifyMesh("Boundconds", model, candidate);
         ERRCHK_ALWAYS(res == AC_SUCCESS);
-        acMeshRandomize(&model);
+        meshRadial(&model);
     }
 
     // Integration
@@ -72,12 +345,21 @@ main(void)
     acGridIntegrate(STREAM_DEFAULT, FLT_EPSILON);
     acGridPeriodicBoundconds(STREAM_DEFAULT);
     acGridStoreMesh(STREAM_DEFAULT, &candidate);
+
+    {
+        std::stringstream ss;
+        ss << "candidate_final_";
+        ss << pid << ".txt";
+        std::cerr << "dump to " << ss.str() << "\n";
+        dumpMesh(&candidate, ss.str());
+    }
+
     if (pid == 0) {
-        acModelIntegrateStep(model, FLT_EPSILON);
+        acHostIntegrateStep(model, FLT_EPSILON);
-        acMeshApplyPeriodicBounds(&model);
+        acHostMeshApplyPeriodicBounds(&model);
         const AcResult res = acVerifyMesh("Integration", model, candidate);
         ERRCHK_ALWAYS(res == AC_SUCCESS);
-        acMeshRandomize(&model);
+        meshRadial(&model);
     }
 
     // Scalar reductions
@@ -93,10 +375,10 @@ main(void)
         AcReal candval;
         acGridReduceScal(STREAM_DEFAULT, (ReductionType)i, v0, &candval);
         if (pid == 0) {
-            const AcReal modelval   = acModelReduceScal(model, (ReductionType)i, v0);
+            const AcReal modelval   = acHostReduceScal(model, (ReductionType)i, v0);
             Error error             = acGetError(modelval, candval);
-            error.maximum_magnitude = acModelReduceScal(model, RTYPE_MAX, v0);
+            error.maximum_magnitude = acHostReduceScal(model, RTYPE_MAX, v0);
-            error.minimum_magnitude = acModelReduceScal(model, RTYPE_MIN, v0);
+            error.minimum_magnitude = acHostReduceScal(model, RTYPE_MIN, v0);
             ERRCHK_ALWAYS(acEvalError(rtype_names[i], error));
         }
     }
@@ -114,17 +396,17 @@ main(void)
         AcReal candval;
         acGridReduceVec(STREAM_DEFAULT, (ReductionType)i, v0, v1, v2, &candval);
         if (pid == 0) {
-            const AcReal modelval   = acModelReduceVec(model, (ReductionType)i, v0, v1, v2);
+            const AcReal modelval   = acHostReduceVec(model, (ReductionType)i, v0, v1, v2);
             Error error             = acGetError(modelval, candval);
-            error.maximum_magnitude = acModelReduceVec(model, RTYPE_MAX, v0, v1, v2);
+            error.maximum_magnitude = acHostReduceVec(model, RTYPE_MAX, v0, v1, v2);
-            error.minimum_magnitude = acModelReduceVec(model, RTYPE_MIN, v0, v1, v1);
+            error.minimum_magnitude = acHostReduceVec(model, RTYPE_MIN, v0, v1, v1);
             ERRCHK_ALWAYS(acEvalError(rtype_names[i], error));
         }
     }
 
     if (pid == 0) {
-        acMeshDestroy(&model);
+        acHostMeshDestroy(&model);
-        acMeshDestroy(&candidate);
+        acHostMeshDestroy(&candidate);
     }
 
     acGridQuit();

samples/standalone_mpi/host_memory.cc

@@ -576,7 +576,7 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
 
     switch (init_type) {
     case INIT_TYPE_RANDOM: {
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         const AcReal range = AcReal(0.01);
         for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
             for (int i = 0; i < n; ++i)
@@ -585,14 +585,14 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
         break;
     }
     case INIT_TYPE_GAUSSIAN_RADIAL_EXPL:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
-        acVertexBufferSet(VTXBUF_LNRHO, mesh->info.real_params[AC_ampl_lnrho], mesh);
+        acHostVertexBufferSet(VTXBUF_LNRHO, mesh->info.real_params[AC_ampl_lnrho], mesh);
         // acmesh_init_to(INIT_TYPE_RANDOM, mesh);
         gaussian_radial_explosion(mesh);
 
         break;
     case INIT_TYPE_XWAVE:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         acmesh_init_to(INIT_TYPE_RANDOM, mesh);
         for (int k = 0; k < mz; k++) {
             for (int j = 0; j < my; j++) {
@@ -605,24 +605,24 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
         }
         break;
     case INIT_TYPE_SIMPLE_CORE:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         simple_uniform_core(mesh);
         break;
     case INIT_TYPE_VEDGE:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         inflow_vedge_freefall(mesh);
         break;
     case INIT_TYPE_VEDGEX:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         inflow_freefall_x(mesh);
         break;
     case INIT_TYPE_RAYLEIGH_TAYLOR:
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         inflow_freefall_x(mesh);
         lnrho_step(mesh);
         break;
     case INIT_TYPE_ABC_FLOW: {
-        acMeshClear(mesh);
+        acHostMeshClear(mesh);
         acmesh_init_to(INIT_TYPE_RANDOM, mesh);
         for (int k = nz_min; k < nz_max; k++) {
             for (int j = ny_min; j < ny_max; j++) {

samples/standalone_mpi/main.cc

@@ -199,9 +199,9 @@ print_diagnostics_host(const AcMesh mesh, const int step, const AcReal dt, const
     const int max_name_width = 16;
 
     // Calculate rms, min and max from the velocity vector field
-    buf_max = acModelReduceVec(mesh, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+    buf_max = acHostReduceVec(mesh, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-    buf_min = acModelReduceVec(mesh, RTYPE_MIN, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+    buf_min = acHostReduceVec(mesh, RTYPE_MIN, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
-    buf_rms = acModelReduceVec(mesh, RTYPE_RMS, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+    buf_rms = acHostReduceVec(mesh, RTYPE_RMS, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
 
     // MV: The ordering in the earlier version was wrong in terms of variable
     // MV: name and its diagnostics.
@@ -213,9 +213,9 @@ print_diagnostics_host(const AcMesh mesh, const int step, const AcReal dt, const
 
     // Calculate rms, min and max from the variables as scalars
     for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
-        buf_max = acModelReduceScal(mesh, RTYPE_MAX, VertexBufferHandle(i));
+        buf_max = acHostReduceScal(mesh, RTYPE_MAX, VertexBufferHandle(i));
-        buf_min = acModelReduceScal(mesh, RTYPE_MIN, VertexBufferHandle(i));
+        buf_min = acHostReduceScal(mesh, RTYPE_MIN, VertexBufferHandle(i));
-        buf_rms = acModelReduceScal(mesh, RTYPE_RMS, VertexBufferHandle(i));
+        buf_rms = acHostReduceScal(mesh, RTYPE_RMS, VertexBufferHandle(i));
 
         printf("  %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, vtxbuf_names[i],
                double(buf_min), double(buf_rms), double(buf_max));
@@ -330,7 +330,7 @@ main(int argc, char** argv)
 
     AcMesh mesh;
     if (pid == 0) {
-        acMeshCreate(info, &mesh);
+        acHostMeshCreate(info, &mesh);
         acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, &mesh);
     }
     acGridInit(info);
@@ -362,7 +362,7 @@ main(int argc, char** argv)
         }
     }
     if (pid == 0)
-        acMeshDestroy(&mesh);
+        acHostMeshDestroy(&mesh);
 
     acGridQuit();
     /////////////// Simple example END
@@ -376,7 +376,7 @@ main(int argc, char** argv)
         if (argc == 3 && (!strcmp(argv[1], "-c") || !strcmp(argv[1], "--config"))) {
             acLoadConfig(argv[2], &info);
             load_config(argv[2], &info);
-            acUpdateBuiltinParams(&info);
+            acHostUpdateBuiltinParams(&info);
         }
         else {
             printf("Usage: ./ac_run\n");
@@ -388,7 +388,7 @@ main(int argc, char** argv)
     else {
         acLoadConfig(AC_DEFAULT_CONFIG, &info);
         load_config(AC_DEFAULT_CONFIG, &info);
-        acUpdateBuiltinParams(&info);
+        acHostUpdateBuiltinParams(&info);
     }
 
     const int start_step     = info.int_params[AC_start_step];
@@ -406,7 +406,7 @@ main(int argc, char** argv)
     AcMesh mesh;
     ///////////////////////////////// PROC 0 BLOCK START ///////////////////////////////////////////
     if (pid == 0) {
-        acMeshCreate(info, &mesh);
+        acHostMeshCreate(info, &mesh);
         // TODO: This need to be possible to define in astaroth.conf
         acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, &mesh);
         // acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); //Initial condition for a collapse test
@@ -445,7 +445,7 @@ main(int argc, char** argv)
 #endif
         }
 
-        acMeshApplyPeriodicBounds(&mesh);
+        acHostMeshApplyPeriodicBounds(&mesh);
         if (start_step == 0) {
             save_mesh(mesh, 0, t_step);
         }
@@ -564,7 +564,7 @@ main(int argc, char** argv)
 
     acGridQuit();
     if (pid == 0)
-        acMeshDestroy(&mesh);
+        acHostMeshDestroy(&mesh);
     fclose(diag_file);
 #endif
 

src/core/astaroth.cc

@@ -74,6 +74,12 @@ acLoad(const AcMesh host_mesh)
     return acNodeLoadMesh(nodes[0], STREAM_DEFAULT, host_mesh);
 }
 
+AcResult
+acSetVertexBuffer(const VertexBufferHandle handle, const AcReal value)
+{
+    return acNodeSetVertexBuffer(nodes[0], STREAM_DEFAULT, handle, value);
+}
+
 AcResult
 acStore(AcMesh* host_mesh)
 {
@@ -182,7 +188,7 @@ acGetNode(void)
 }
 
 AcResult
-acUpdateBuiltinParams(AcMeshInfo* config)
+acHostUpdateBuiltinParams(AcMeshInfo* config)
 {
     config->int_params[AC_mx] = config->int_params[AC_nx] + STENCIL_ORDER;
     ///////////// PAD TEST
@@ -221,7 +227,7 @@ acUpdateBuiltinParams(AcMeshInfo* config)
 }
 
 AcResult
-acMeshCreate(const AcMeshInfo info, AcMesh* mesh)
+acHostMeshCreate(const AcMeshInfo info, AcMesh* mesh)
 {
     mesh->info = info;
 
@@ -241,7 +247,7 @@ randf(void)
 }
 
 AcResult
-acMeshRandomize(AcMesh* mesh)
+acHostMeshRandomize(AcMesh* mesh)
 {
     const int n = acVertexBufferSize(mesh->info);
     for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
@@ -252,7 +258,7 @@ acMeshRandomize(AcMesh* mesh)
 }
 
 AcResult
-acMeshDestroy(AcMesh* mesh)
+acHostMeshDestroy(AcMesh* mesh)
 {
     for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
         free(mesh->vertex_buffer[w]);

src/core/astaroth_fortran.cc

@@ -7,9 +7,9 @@
  * Utils
  */
 void
-acupdatebuiltinparams_(AcMeshInfo* info)
+achostupdatebuiltinparams_(AcMeshInfo* info)
 {
-    acUpdateBuiltinParams(info);
+    acHostUpdateBuiltinParams(info);
 }
 
 void

src/core/astaroth_fortran.h

@@ -8,7 +8,7 @@ extern "C" {
 /**
  * Utils
  */
-void acupdatebuiltinparams_(AcMeshInfo* info);
+void achostupdatebuiltinparams_(AcMeshInfo* info);
 
 void acgetdevicecount_(int* count);
 

src/core/device.cc

@@ -294,6 +294,26 @@ acDeviceLoadMesh(const Device device, const Stream stream, const AcMesh host_mes
     return AC_SUCCESS;
 }
 
+AcResult
+acDeviceSetVertexBuffer(const Device device, const Stream stream, const VertexBufferHandle handle, const AcReal value)
+{
+    acDeviceSynchronizeStream(device, stream);
+
+    const size_t count = acVertexBufferSize(device->local_config);
+    AcReal* data = (AcReal*) malloc(sizeof(AcReal) * count);
+    ERRCHK_ALWAYS(data);
+
+    for (size_t i = 0; i < count; ++i)
+        data[i] = value;
+
+    // Set both in and out for safety (not strictly needed)
+    ERRCHK_CUDA_ALWAYS(cudaMemcpyAsync(device->vba.in[handle], data, sizeof(data[0]) * count, cudaMemcpyHostToDevice, device->streams[stream]));
+    ERRCHK_CUDA_ALWAYS(cudaMemcpyAsync(device->vba.out[handle], data, sizeof(data[0]) * count, cudaMemcpyHostToDevice, device->streams[stream]));
+
+    free(data);
+    return AC_SUCCESS;
+}
+
 AcResult
 acDeviceStoreVertexBufferWithOffset(const Device device, const Stream stream,
                                     const VertexBufferHandle vtxbuf_handle, const int3 src,
@@ -1269,10 +1289,10 @@ acGridRandomize(void)
     ERRCHK(grid.initialized);
 
     AcMesh host;
-    acMeshCreate(grid.submesh.info, &host);
+    acHostMeshCreate(grid.submesh.info, &host);
-    acMeshRandomize(&host);
+    acHostMeshRandomize(&host);
     acDeviceLoadMesh(grid.device, STREAM_DEFAULT, host);
-    acMeshDestroy(&host);
+    acHostMeshDestroy(&host);
 
     return AC_SUCCESS;
 }
@@ -1320,7 +1340,7 @@ acGridInit(const AcMeshInfo info)
     };
     submesh_info.int3_params[AC_multigpu_offset] = pid3d *
                                                    (int3){submesh_nx, submesh_ny, submesh_nz};
-    acUpdateBuiltinParams(&submesh_info);
+    acHostUpdateBuiltinParams(&submesh_info);
 
     // GPU alloc
     int devices_per_node = -1;
@@ -1331,7 +1351,7 @@ acGridInit(const AcMeshInfo info)
 
     // CPU alloc
     AcMesh submesh;
-    acMeshCreate(submesh_info, &submesh);
+    acHostMeshCreate(submesh_info, &submesh);
 
     // Setup the global grid structure
     grid.device        = device;
@@ -1380,7 +1400,7 @@ acGridQuit(void)
 
     grid.initialized   = false;
     grid.decomposition = (uint3_64){0, 0, 0};
-    acMeshDestroy(&grid.submesh);
+    acHostMeshDestroy(&grid.submesh);
     acDeviceDestroy(grid.device);
 
     acGridSynchronizeStream(STREAM_ALL);

src/core/node.cc

@@ -536,6 +536,18 @@ acNodeLoadMesh(const Node node, const Stream stream, const AcMesh host_mesh)
     return AC_SUCCESS;
 }
 
+AcResult
+acNodeSetVertexBuffer(const Node node, const Stream stream, const VertexBufferHandle handle, const AcReal value)
+{
+    acNodeSynchronizeStream(node, stream);
+
+    for (int i = 0; i < node->num_devices; ++i)
+        acDeviceSetVertexBuffer(node->devices[i], stream, handle, value);
+
+    acNodeSynchronizeStream(node, stream); // For safety
+    return AC_SUCCESS;
+}
+
 AcResult
 acNodeStoreVertexBufferWithOffset(const Node node, const Stream stream,
                                   const VertexBufferHandle vtxbuf_handle, const int3 src,

src/utils/CMakeLists.txt

@@ -1,3 +1,4 @@
 ## Astaroth Utils
 add_library(astaroth_utils STATIC config_loader.c memory.c verification.c modelsolver.c modelreduce.c)
 add_dependencies(astaroth_utils dsl_headers)
+target_compile_options(astaroth_utils PRIVATE "-mavx")

src/utils/config_loader.c

@@ -89,7 +89,7 @@ acLoadConfig(const char* config_path, AcMeshInfo* config)
     memset(config, (uint8_t)0xFF, sizeof(*config));
 
     parse_config(config_path, config);
-    acUpdateBuiltinParams(config);
+    acHostUpdateBuiltinParams(config);
 #if AC_VERBOSE
     printf("###############################################################\n");
     printf("Config dimensions loaded:\n");

src/utils/memory.c

@@ -21,7 +21,7 @@
 #include "errchk.h"
 
 AcResult
-acVertexBufferSet(const VertexBufferHandle handle, const AcReal value, AcMesh* mesh)
+acHostVertexBufferSet(const VertexBufferHandle handle, const AcReal value, AcMesh* mesh)
 {
     const int n = acVertexBufferSize(mesh->info);
     for (int i = 0; i < n; ++i)
@@ -30,16 +30,16 @@ acVertexBufferSet(const VertexBufferHandle handle, const AcReal value, AcMesh* m
     return AC_SUCCESS;
 }
 AcResult
-acMeshSet(const AcReal value, AcMesh* mesh)
+acHostMeshSet(const AcReal value, AcMesh* mesh)
 {
     for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
-        acVertexBufferSet(w, value, mesh);
+        acHostVertexBufferSet(w, value, mesh);
 
     return AC_SUCCESS;
 }
 
 AcResult
-acMeshApplyPeriodicBounds(AcMesh* mesh)
+acHostMeshApplyPeriodicBounds(AcMesh* mesh)
 {
     const AcMeshInfo info = mesh->info;
     for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
@@ -105,7 +105,7 @@ acMeshApplyPeriodicBounds(AcMesh* mesh)
 }
 
 AcResult
-acMeshClear(AcMesh* mesh)
+acHostMeshClear(AcMesh* mesh)
 {
-    return acMeshSet(0, mesh);
+    return acHostMeshSet(0, mesh);
 }

src/utils/modelreduce.c

@@ -74,7 +74,7 @@ exp_squared_vec(const AcReal a, const AcReal b, const AcReal c) { return exp_squ
 // clang-format on
 
 AcReal
-acModelReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a)
+acHostReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a)
 {
     ReduceInitialScalFunc reduce_initial;
     ReduceFunc reduce;
@@ -139,7 +139,7 @@ acModelReduceScal(const AcMesh mesh, const ReductionType rtype, const VertexBuff
 }
 
 AcReal
-acModelReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a,
+acHostReduceVec(const AcMesh mesh, const ReductionType rtype, const VertexBufferHandle a,
                  const VertexBufferHandle b, const VertexBufferHandle c)
 {
     // AcReal (*reduce_initial)(AcReal, AcReal, AcReal);

src/utils/modelsolver.c

@@ -30,7 +30,7 @@
 #include <stdbool.h>
 
 #include "errchk.h"
-#include "memory.h" // acMeshCreate, acMeshDestroy, acMeshApplyPeriodicBounds
+#include "memory.h" // acHostMeshCreate, acHostMeshDestroy, acHostMeshApplyPeriodicBounds
 
 // Standalone flags
 #define LDENSITY (1)
@@ -985,7 +985,7 @@ checkConfiguration(const AcMeshInfo info)
 }
 
 AcResult
-acModelIntegrateStep(AcMesh mesh, const AcReal dt)
+acHostIntegrateStep(AcMesh mesh, const AcReal dt)
 {
     mesh_info = &(mesh.info);
 
@@ -998,7 +998,7 @@ acModelIntegrateStep(AcMesh mesh, const AcReal dt)
     checkConfiguration(*mesh_info);
 
     AcMesh intermediate_mesh;
-    acMeshCreate(mesh.info, &intermediate_mesh);
+    acHostMeshCreate(mesh.info, &intermediate_mesh);
 
     const int nx_min = getInt(AC_nx_min);
     const int nx_max = getInt(AC_nx_max);
@@ -1012,7 +1012,7 @@ acModelIntegrateStep(AcMesh mesh, const AcReal dt)
     for (int step_number = 0; step_number < 3; ++step_number) {
 
         // Boundconds
-        acMeshApplyPeriodicBounds(&mesh);
+        acHostMeshApplyPeriodicBounds(&mesh);
 
         // Alpha step
         // #pragma omp parallel for
@@ -1035,7 +1035,7 @@ acModelIntegrateStep(AcMesh mesh, const AcReal dt)
         }
     }
 
-    acMeshDestroy(&intermediate_mesh);
+    acHostMeshDestroy(&intermediate_mesh);
     mesh_info = NULL;
     return AC_SUCCESS;
 }