diff --git a/acc/mhd_solver/stencil_kernel.ac b/acc/mhd_solver/stencil_kernel.ac index 0f37ea5..e0efa94 100644 --- a/acc/mhd_solver/stencil_kernel.ac +++ b/acc/mhd_solver/stencil_kernel.ac @@ -2,12 +2,15 @@ #define LDENSITY (1) #define LHYDRO (1) +// MV: Currenly only magnetic with entropy. Support for isothermal MHD required +// MV: (matter of switch combination). #define LMAGNETIC (1) #define LENTROPY (1) #define LTEMPERATURE (0) #define LFORCING (0) #define LUPWD (0) #define LSINK (0) +#define LBFIELD (0) #define AC_THERMAL_CONDUCTIVITY (0.001) // TODO: make an actual config parameter #define H_CONST (0) // TODO: make an actual config parameter @@ -74,6 +77,9 @@ uniform Scalar AC_ampl_uu; uniform Scalar AC_angl_uu; uniform Scalar AC_lnrho_edge; uniform Scalar AC_lnrho_out; +uniform Scalar AC_ampl_aa; +uniform Scalar AC_init_k_wave; +uniform Scalar AC_init_sigma_hel; // Forcing parameters. User configured. uniform Scalar AC_forcing_magnitude; uniform Scalar AC_relhel; @@ -133,6 +139,12 @@ uniform ScalarField VTXBUF_LNRHO; uniform ScalarField VTXBUF_ACCRETION; #endif +#if LBFIELD +uniform ScalarField BFIELDX; +uniform ScalarField BFIELDY; +uniform ScalarField BFIELDZ; +#endif + #if LUPWD Preprocessed Scalar @@ -449,11 +461,14 @@ induction(in VectorField uu, in VectorField aa) // yes this actually works. See pg.28 in arXiv:astro-ph/0109497) // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ]) const Vector B = curl(aa); - const Vector grad_div = gradient_of_divergence(aa); + //MV: Due to gauge freedom we can reduce the gradient of scalar (divergence) from the equation + //const Vector grad_div = gradient_of_divergence(aa); const Vector lap = laplace_vec(aa); // Note, AC_mu0 is cancelled out - const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap); + //MV: Due to gauge freedom we can reduce the gradient of scalar (divergence) from the equation + //const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap); + const Vector ind = cross(value(uu), B) + AC_eta * lap; return ind; } @@ -626,6 +641,15 @@ forcing(int3 globalVertexIdx, Scalar dt) } #endif // LFORCING +#if LBFIELD +// Calculate the B-field for VTXBUFF +Device Vector +get_bfield(in VectorField aa) +{ + return curl(aa); +} +#endif + // Declare input and output arrays using locations specified in the // array enum in astaroth.h in ScalarField lnrho(VTXBUF_LNRHO); @@ -654,6 +678,11 @@ in ScalarField accretion(VTXBUF_ACCRETION); out ScalarField out_accretion(VTXBUF_ACCRETION); #endif +#if LBFIELD +out VectorField b_out(BFIELDX, BFIELDY, BFIELDZ); +#endif + + Kernel void solve() { @@ -688,4 +717,11 @@ solve() out_accretion = out_accretion * AC_dsx * AC_dsy * AC_dsz; // unit is now mass! } #endif + +#if LBFIELD + if (step_number == 2) { + b_out = get_bfield(aa); + } +#endif + } diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c index b0dcb13..017591a 100644 --- a/acc/src/code_generator.c +++ b/acc/src/code_generator.c @@ -736,6 +736,9 @@ external acdevicesynchronizestream fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_MIN"); fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_RMS"); fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_RMS_EXP"); + fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_ALFVEN_MAX"); + fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_ALFVEN_MIN"); + fprintf(DSLHEADER, "FUNC(%s)\\\n", "RTYPE_ALFVEN_RMS"); fprintf(DSLHEADER, "FUNC(%s)\n", "RTYPE_SUM"); size_t counter = 0; @@ -747,6 +750,12 @@ external acdevicesynchronizestream ++counter; fprintf(FHEADER, "integer(c_int), parameter :: RTYPE_RMS_EXP = %lu\n", counter); ++counter; + fprintf(FHEADER, "integer(c_int), parameter :: RTYPE_ALFVEN_MAX = %lu\n", counter); + ++counter; + fprintf(FHEADER, "integer(c_int), parameter :: RTYPE_ALFVEN_MIN = %lu\n", counter); + ++counter; + fprintf(FHEADER, "integer(c_int), parameter :: RTYPE_ALFVEN_RMS = %lu\n", counter); + ++counter; fprintf(FHEADER, "integer(c_int), parameter :: RTYPE_SUM = %lu\n", counter); ++counter; fprintf(FHEADER, "integer(c_int), parameter :: NUM_REDUCTION_TYPES = %lu\n", counter); diff --git a/analysis/python/astar/data/read.py b/analysis/python/astar/data/read.py index 9b2b820..3ae6ff8 100644 --- a/analysis/python/astar/data/read.py +++ b/analysis/python/astar/data/read.py @@ -20,6 +20,7 @@ # This module is for reading data. import numpy as np +import os #Optional YT interface try: @@ -28,26 +29,28 @@ try: except ImportError: yt_present = False -def set_dtype(endian, AcRealSize): +def set_dtype(endian, AcRealSize, print_type = True): if endian == 0: en = '>' elif endian == 1: en = '<' type_instruction = en + 'f' + str(AcRealSize) - print("type_instruction", type_instruction) + if print_type: + print("type_instruction", type_instruction) my_dtype = np.dtype(type_instruction) return my_dtype -def read_bin(fname, fdir, fnum, minfo, numtype=np.longdouble): +def read_bin(fname, fdir, fnum, minfo, numtype=np.longdouble, getfilename=True): '''Read in a floating point array''' filename = fdir + fname + '_' + fnum + '.mesh' datas = np.DataSource() read_ok = datas.exists(filename) - my_dtype = set_dtype(minfo.contents['endian'], minfo.contents['AcRealSize']) + my_dtype = set_dtype(minfo.contents['endian'], minfo.contents['AcRealSize'], print_type=getfilename) if read_ok: - print(filename) + if getfilename: + print(filename) array = np.fromfile(filename, dtype=my_dtype) timestamp = array[0] @@ -92,10 +95,58 @@ def read_meshtxt(fdir, fname, dbg_output): print(line[1], contents[line[1]]) else: print(line) - print('ERROR: ' + line[0] +' not recognized!') + if dbg_output: + print('ERROR: ' + line[0] +' not recognized!') return contents +def parse_directory(meshdir): + dirlist = os.listdir(meshdir) + dirlist = [k for k in dirlist if 'LNRHO' in k] + for i, item in enumerate(dirlist): + tmp = item.strip('.mesh') + tmp = tmp.strip('VTXBUF_LNRHO') + dirlist[i] = int(tmp) + dirlist.sort() + return dirlist + +def apply_boundcond(array, btype): + + if btype == "p": + be = 3 + bi = 6 + # Edges + # xx + array[ : 3, : , : ] = array[-6:-3, : , : ] + array[-3: , : , : ] = array[ 3: 6, : , : ] + # yy + array[ : , : 3, : ] = array[ : ,-6:-3, : ] + array[ : ,-3: , : ] = array[ : , 3: 6, : ] + # zz + array[ : , : , : 3] = array[ : , : ,-6:-3] + array[ : , : ,-3: ] = array[ : , : , 3: 6] + # Corner parts + # xy + array[ : 3, : 3, : ] = array[-6:-3,-6:-3, : ] + array[-3: ,-3: , : ] = array[ 3: 6, 3: 6, : ] + array[-3: , : 3, : ] = array[ 3: 6,-6:-3, : ] + array[ : 3,-3: , : ] = array[-6:-3, 3: 6, : ] + # xz + array[ : 3, : , : 3] = array[-6:-3, : ,-6:-3] + array[-3: , : ,-3: ] = array[ 3: 6, : , 3: 6] + array[-3: , : , : 3] = array[ 3: 6, : ,-6:-3] + array[ : 3, : ,-3: ] = array[-6:-3, : , 3: 6] + # yz + array[ : , : 3, : 3] = array[ : ,-6:-3,-6:-3] + array[ : ,-3: ,-3: ] = array[ : , 3: 6, 3: 6] + array[ : ,-3: , : 3] = array[ : , 3: 6,-6:-3] + array[ : , : 3,-3: ] = array[ : ,-6:-3, 3: 6] + else: + print("Unknown btype", btype) + + return array + + def DERX(array, dx): output = np.zeros_like(array) for i in range(3, array.shape[0]-3): #Keep boundary poits as 0 @@ -120,6 +171,34 @@ def DERZ(array, dz): - array[:,:,i-3] + array[:,:,i+3] )/(60.0*dz) return output +def DER2X(array, dx): + output = np.zeros_like(array) + for i in range(3, array.shape[0]-3): #Keep boundary poits as 0 + output[i,:,:] =( 2.0*array[i-1,:,:] + 2.0*array[i+1,:,:] + - 27.0*array[i-2,:,:] - 27.0*array[i+2,:,:] + +270.0*array[i-3,:,:] + 270.0*array[i+3,:,:] + -490.0*array[i ,:,:] )/(180.0*dx*dx) + return output + +def DER2Y(array, dy): + output = np.zeros_like(array) + for i in range(3,array.shape[1]-3): + output[:,i,:] =( 2.0*array[:,i-1,:] + 2.0*array[:,i+1,:] + - 27.0*array[:,i-2,:] - 27.0*array[:,i+2,:] + +270.0*array[:,i-3,:] + 270.0*array[:,i+3,:] + -490.0*array[:,i ,:] )/(180.0*dy*dy) + return output + +def DER2Z(array, dz): + output = np.zeros_like(array) + for i in range(3, array.shape[2]-3): + output[:,:,i] =( 2.0*array[:,:,i-1] + 2.0*array[:,:,i+1] + - 27.0*array[:,:,i-2] - 27.0*array[:,:,i+2] + +270.0*array[:,:,i-3] + 270.0*array[:,:,i+3] + -490.0*array[:,:,i ] )/(180.0*dz*dz) + return output + + def curl(aa, minfo): dx = minfo.contents['AC_dsx'] dy = minfo.contents['AC_dsy'] @@ -128,6 +207,44 @@ def curl(aa, minfo): DERZ(aa[0], dz)-DERX(aa[2], dx), DERX(aa[1], dx)-DERY(aa[0], dy)) +def div(array, minfo): + dx = minfo.contents['AC_dsx'] + dy = minfo.contents['AC_dsy'] + dz = minfo.contents['AC_dsz'] + return ( DERX(array[0], dx) + + DERY(array[1], dy) + + DERZ(array[2], dz)) + +def grad(array, minfo): + dx = minfo.contents['AC_dsx'] + dy = minfo.contents['AC_dsy'] + dz = minfo.contents['AC_dsz'] + return (DERX(array, dx), + DERY(array, dy), + DERZ(array, dz)) + +def grad_div(array, minfo): + scalar = div(array, minfo) + scalar = apply_boundcond(scalar, "p") + vec = grad(scalar, minfo) + return vec + +def laplace_scal(array, minfo): + dx = minfo.contents['AC_dsx'] + dy = minfo.contents['AC_dsy'] + dz = minfo.contents['AC_dsz'] + return (DER2X(array, dx) + DER2Y(array, dy) + DER2Z(array, dz)) + +def laplace_vec(array, minfo): + return (laplace_scal(array[0], minfo), + laplace_scal(array[1], minfo), + laplace_scal(array[2], minfo)) + +def curl_of_curl(array, minfo): + array1 = curl(array, minfo) + array2 = (apply_boundcond(array1[0], "p"), apply_boundcond(array1[1], "p"), apply_boundcond(array1[2], "p")) + return curl(array2, minfo) + class MeshInfo(): '''Object that contains all mesh info''' @@ -138,45 +255,40 @@ class MeshInfo(): class Mesh: '''Class tha contains all 3d mesh data''' - def __init__(self, fnum, fdir=""): + def __init__(self, fnum, fdir="", only_info = False, pdiag = True): fnum = str(fnum) self.framenum = fnum.zfill(10) self.minfo = MeshInfo(fdir) - self.lnrho, self.timestamp, self.ok = read_bin('VTXBUF_LNRHO', fdir, fnum, self.minfo) + if only_info == False: + self.lnrho, self.timestamp, self.ok = read_bin('VTXBUF_LNRHO', fdir, fnum, self.minfo, getfilename=pdiag) + else: + self.ok = False if self.ok: - self.ss, timestamp, ok = read_bin('VTXBUF_ENTROPY', fdir, fnum, self.minfo) + self.ss, timestamp, ok = read_bin('VTXBUF_ENTROPY', fdir, fnum, self.minfo, getfilename=pdiag) - self.accretion, timestamp, ok = read_bin('VTXBUF_ACCRETION', fdir, fnum, self.minfo) + self.accretion, timestamp, ok = read_bin('VTXBUF_ACCRETION', fdir, fnum, self.minfo, getfilename=pdiag) #TODO Generalize is a dict. Do not hardcode! - uux, timestamp, ok = read_bin('VTXBUF_UUX', fdir, fnum, self.minfo) - uuy, timestamp, ok = read_bin('VTXBUF_UUY', fdir, fnum, self.minfo) - uuz, timestamp, ok = read_bin('VTXBUF_UUZ', fdir, fnum, self.minfo) + uux, timestamp, ok = read_bin('VTXBUF_UUX', fdir, fnum, self.minfo, getfilename=pdiag) + uuy, timestamp, ok = read_bin('VTXBUF_UUY', fdir, fnum, self.minfo, getfilename=pdiag) + uuz, timestamp, ok = read_bin('VTXBUF_UUZ', fdir, fnum, self.minfo, getfilename=pdiag) self.uu = (uux, uuy, uuz) uux = [] uuy = [] uuz = [] - aax, timestamp, ok = read_bin('VTXBUF_AX', fdir, fnum, self.minfo) - aay, timestamp, ok = read_bin('VTXBUF_AY', fdir, fnum, self.minfo) - aaz, timestamp, ok = read_bin('VTXBUF_AZ', fdir, fnum, self.minfo) + aax, timestamp, ok = read_bin('VTXBUF_AX', fdir, fnum, self.minfo, getfilename=pdiag) + aay, timestamp, ok = read_bin('VTXBUF_AY', fdir, fnum, self.minfo, getfilename=pdiag) + aaz, timestamp, ok = read_bin('VTXBUF_AZ', fdir, fnum, self.minfo, getfilename=pdiag) self.aa = (aax, aay, aaz) aax = [] aay = [] aaz = [] - - #self.aa[0][:,:,:] = 0.0 - #self.aa[1][:,:,:] = 0.0 - #self.aa[2][:,:,:] = 0.0 - #for i in range(0, self.aa[0].shape[0]): - # self.aa[0][:,i,:] = float(i) - - self.xx = np.arange(self.minfo.contents['AC_mx']) * self.minfo.contents['AC_dsx'] self.yy = np.arange(self.minfo.contents['AC_my']) * self.minfo.contents['AC_dsy'] self.zz = np.arange(self.minfo.contents['AC_mz']) * self.minfo.contents['AC_dsz'] @@ -184,10 +296,44 @@ class Mesh: self.xmid = int(self.minfo.contents['AC_mx']/2) self.ymid = int(self.minfo.contents['AC_my']/2) self.zmid = int(self.minfo.contents['AC_mz']/2) - def Bfield(self, get_jj = False): + def Bfield(self, get_jj = False, trim=False): self.bb = curl(self.aa, self.minfo) if get_jj: - self.jj = curl(self.bb, self.minfo) + 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]) + 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 get_jj(self, trim=False): + self.jj = curl_of_curl(self.aa, minfo, trim=False) + if trim: + 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 vorticity(self, trim=False): + self.oo = curl(self.uu, self.minfo) + if trim: + self.oo = (self.oo[0][3:-3, 3:-3, 3:-3],self.oo[1][3:-3, 3:-3, 3:-3],self.oo[2][3:-3, 3:-3, 3:-3]) + + def rad_vel(self): + print("Calculating spherical velocity components") + self.uu_pherical = np.zeros_like(self.uu) + xx, yy, zz = np.meshgrid(self.xx - self.xmid, self.yy - self.ymid, self.zz - self.zmid) + rr = np.sqrt(xx**2.0 + yy**2.0 + zz**2.0) + theta = np.arccos(zz/rr) + phi = np.arctan2(yy,xx) + sin_theta_sin_phi = np.sin(theta)*np.sin(phi) + cos_theta_cos_phi = np.cos(theta)*np.cos(phi) + sin_theta_cos_phi = np.sin(theta)*np.cos(phi) + cos_theta_sin_phi = np.cos(theta)*np.sin(phi) + ux = self.uu[0]; uy = self.uu[1]; uz = self.uu[2]; + vr = sin_theta_cos_phi*ux + sin_theta_sin_phi*uy + np.cos(theta)*uz + vtheta = cos_theta_cos_phi*ux + cos_theta_sin_phi*uy - np.sin(theta)*uz + vphi = -np.sin(phi)*ux + np.cos(phi)*uy + self.uu_pherical[0] = vr + self.uu_pherical[1] = vtheta + self.uu_pherical[2] = vphi + def yt_conversion(self): if yt_present: @@ -266,42 +412,143 @@ class Mesh: f.write(binary_format) f.close() + def export_vtk_ascii(self, Beq = 1.0): + #BASED ON https://lorensen.github.io/VTKExamples/site/VTKFileFormats/#dataset-attribute-format + self.Bfield() + + f = open("GRID%s.vtk" % self.framenum, 'w') + + Ntot = self.minfo.contents['AC_mx']*self.minfo.contents['AC_my']*self.minfo.contents['AC_mz'] + + mx = self.minfo.contents['AC_mx'] + my = self.minfo.contents['AC_my'] + mz = self.minfo.contents['AC_mz'] + + print("Writing GRID%s.vtk" % self.framenum) + + f.write("# vtk DataFile Version 2.0\n") + f.write("Astaroth grid for visualization\n") + f.write("ASCII\n") + #f.write("DATASET STRUCTURED_GRID\n") + #f.write("DIMENSIONS %i %i %i \n" % (mx, my, mz)) + #f.write("POINTS %i float \n" % (Ntot)) + #for i in range(mx): + # for j in range(my): + # for k in range(mz): + # f.write("%e %e %e \n" % (i, j, k)) + f.write("DATASET RECTILINEAR_GRID\n") + f.write("DIMENSIONS %i %i %i \n" % (mx, my, mz)) + f.write("X_COORDINATES %i float \n" % mx) + for i in range(mx): + f.write("%e " % (i)) + f.write("\n") + f.write("Y_COORDINATES %i float \n" % my) + for j in range(my): + f.write("%e " % (j)) + f.write("\n") + f.write("Z_COORDINATES %i float \n" % mz) + for k in range(mz): + f.write("%e " % (k)) + f.write("\n") + f.write("POINT_DATA %i \n" % (Ntot)) + f.write("VECTORS velocity float \n" ) + for i in range(mx): + if (i % 8) == 0: + print("i = %i / %i" %(i, mx-1)) + for j in range(my): + for k in range(mz): + f.write("%e %e %e \n" % ( self.uu[0][i, j, k], self.uu[1][i, j, k], self.uu[2][i, j, k])) + #f.write("POINT_DATA %i \n" % (Ntot)) + f.write("VECTORS bfield float \n") + eqprint = True + for i in range(mx): + if (i % 8) == 0: + print("i = %i / %i" %(i, mx-1)) + eqprint = True + for j in range(my): + for k in range(mz): + #Beq is the equipartition magnetic field. + while(eqprint): + print("normal B %e %e %e \n" % (self.bb[0][i, j, k], self.bb[1][i, j, k], self.bb[2][i, j, k] )) + print("equipartition B %e %e %e \n" % (self.bb[0][i, j, k]/Beq, self.bb[1][i, j, k]/Beq, self.bb[2][i, j, k]/Beq)) + eqprint = False + f.write("%e %e %e \n" % ( self.bb[0][i, j, k]/Beq, self.bb[1][i, j, k]/Beq, self.bb[2][i, j, k]/Beq)) + #ADD DENSITY SCALAR + f.write("\n") + + print("Done.") + + f.close() + +def find_explosion_index(array, criterion = 1e5): + for i, ar in enumerate(array): + if (np.abs(array[i])-np.abs(array[i-1])) > criterion: + return i + return -1 + +def mask_bad_values_ts(ts, criterion = 1e5): + indexmask = np.zeros_like(ts) + + ts = np.ma.masked_invalid(ts) + + index = find_explosion_index(ts, criterion = criterion) + + if index >= 0: + indexmask[index:] = 1 + + ts = np.ma.array(ts, mask=indexmask) + + return ts -def parse_ts(fdir, fname): - with open(fdir+fname) as f: - filetext = f.read().splitlines() - +def parse_ts(fdir, fname, debug = False): var = {} - line = filetext[0].split() - for i in range(len(line)): - line[i] = line[i].replace('VTXBUF_', "") - line[i] = line[i].replace('UU', "uu") - line[i] = line[i].replace('_total', "tot") - line[i] = line[i].replace('ACCRETION', "acc") - line[i] = line[i].replace('A', "aa") - line[i] = line[i].replace('LNRHO', "lnrho") - line[i] = line[i].replace('ENTROPY', "ss") - line[i] = line[i].replace('X', "x") - line[i] = line[i].replace('Y', "y") - line[i] = line[i].replace('Z', "z") + tsfile = fdir+fname - tsdata = np.loadtxt(fdir+fname,skiprows=1) + if os.path.exists(tsfile): - for i in range(len(line)): - var[line[i]] = tsdata[:,i] + with open(tsfile) as f: + filetext = f.read().splitlines() - var['step'] = np.int64(var['step']) + line = filetext[0].split() + for i in range(len(line)): + line[i] = line[i].replace('VTXBUF_', "") + line[i] = line[i].replace('UU', "uu") + line[i] = line[i].replace('_total', "tot") + line[i] = line[i].replace('ACCRETION', "acc") + line[i] = line[i].replace('A', "aa") + line[i] = line[i].replace('LNRHO', "lnrho") + line[i] = line[i].replace('ENTROPY', "ss") + line[i] = line[i].replace('BFIELD', "bb") + line[i] = line[i].replace('X', "x") + line[i] = line[i].replace('Y', "y") + line[i] = line[i].replace('Z', "z") + line[i] = line[i].replace('vaa', "vA") - print("HERE ARE ALL KEYS FOR TS DATA:") - print(var.keys()) + #tsdata = np.loadtxt(fdir+fname,skiprows=1) + tsdata = np.genfromtxt(fdir+fname,skip_header=1, skip_footer=1) + + for i in range(len(line)): + var[line[i]] = tsdata[:,i] + + var['step'] = np.int64(var['step']) + + var['exist'] = True + + else: + var['exist'] = False + + if debug: + print("HERE ARE ALL KEYS FOR TS DATA:") + print(var.keys()) return var + class TimeSeries: '''Class for time series data''' - def __init__(self, fdir="", fname="timeseries.ts"): + def __init__(self, fdir="", fname="timeseries.ts", debug = False): - self.var = parse_ts(fdir, fname) + self.var = parse_ts(fdir, fname, debug = debug) diff --git a/analysis/python/astar/visual/lineplot.py b/analysis/python/astar/visual/lineplot.py index e793a6f..b76674a 100644 --- a/analysis/python/astar/visual/lineplot.py +++ b/analysis/python/astar/visual/lineplot.py @@ -39,7 +39,7 @@ def plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3): def plot_ts(ts, isotherm=False, show_all=False, lnrho=False, uutot=False, uux=False, uuy=False, uuz=False, aax=False, aay=False, aaz=False, - ss=False, acc=False, sink=False, rho=False): + ss=False, acc=False, sink=False, rho=False, bb=False, alfven=False): if show_all: lnrho=True @@ -51,6 +51,8 @@ def plot_ts(ts, isotherm=False, show_all=False, lnrho=False, uutot=False, ss=True acc=True sink=True + bb=True + alfven=True if isotherm: lnrho=True @@ -143,6 +145,44 @@ def plot_ts(ts, isotherm=False, show_all=False, lnrho=False, uutot=False, yaxis3 = 'ss_max' plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + if bb: + plt.figure() + xaxis = 't_step' + yaxis1 = 'bbtot_rms' + yaxis2 = 'bbtot_min' + yaxis3 = 'bbtot_max' + plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + + plt.figure() + xaxis = 't_step' + yaxis1 = 'bbx_rms' + yaxis2 = 'bbx_min' + yaxis3 = 'bbx_max' + plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + + plt.figure() + xaxis = 't_step' + yaxis1 = 'bby_rms' + yaxis2 = 'bby_min' + yaxis3 = 'bby_max' + plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + + plt.figure() + xaxis = 't_step' + yaxis1 = 'bbz_rms' + yaxis2 = 'bbz_min' + yaxis3 = 'bbz_max' + plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + + if alfven: + plt.figure() + xaxis = 't_step' + yaxis1 = 'vAtot_rms' + yaxis2 = 'vAtot_min' + yaxis3 = 'vAtot_max' + plot_min_man_rms(ts, xaxis, yaxis1, yaxis2, yaxis3) + + if acc: plt.figure() xaxis = 't_step' diff --git a/analysis/python/samples/readtest.py b/analysis/python/samples/readtest.py index 72e5107..ddede1a 100644 --- a/analysis/python/samples/readtest.py +++ b/analysis/python/samples/readtest.py @@ -22,6 +22,9 @@ import pylab as plt import numpy as np import sys +import os +import pandas as pd + ##mesh = ad.read.Mesh(500, fdir="/tiara/home/mvaisala/astaroth-code/astaroth_2.0/build/") ## ##print(np.shape(mesh.uu)) @@ -44,6 +47,23 @@ print("sys.argv", sys.argv) meshdir = "$HOME/astaroth/build/" + +#Example fixed scaling template +if (meshdir == "$HOME/astaroth/build/"): + rlnrho = [- 0.08, 0.08] + rrho = [ 0.93, 1.08] + rNcol = [ 500.0, 530.0] + + ruu_tot = [ 0.0, 0.3] + ruu_xyz = [-0.3, 0.3] + + raa_tot = [ 0.0, 0.14] + raa_xyz = [-0.14, 0.14] + + rbb_tot = [ 0.0, 0.3] + rbb_xyz = [-0.3, 0.3] + + if "xtopbound" in sys.argv: for i in range(0, 171): mesh = ad.read.Mesh(i, fdir=meshdir) @@ -178,23 +198,37 @@ if 'sl' in sys.argv: mesh.Bfield() bb_tot = np.sqrt(mesh.bb[0]**2.0 + mesh.bb[1]**2.0 + mesh.bb[2]**2.0) - if 'lim' in sys.argv: - vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho', colrange=[-0.0, 3.5]) - vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho', colrange=[1.0, 5.0 ]) - vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux', colrange=[-5.0, 5.0]) - vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$', bitmap = True, fname = 'uuy', colrange=[-5.0, 5.0]) - vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$', bitmap = True, fname = 'uuz', colrange=[-5.0, 5.0]) - vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\mathrm{col}$', bitmap = True, fname = 'colden', slicetype = 'sum', colrange=[0.0, 200.0]) - vis.slices.plot_3(mesh, uu_tot, title = r'$|u|$', bitmap = True, fname = 'uutot', colrange=[0.00, 5.0]) - vis.slices.plot_3(mesh, aa_tot, title = r'$\|A\|$', bitmap = True, fname = 'aatot', colrange=[0.0,0.01]) - vis.slices.plot_3(mesh, mesh.aa[0], title = r'$A_x$', bitmap = True, fname = 'aax', colrange=[-0.01,0.01]) - vis.slices.plot_3(mesh, mesh.aa[1], title = r'$A_y$', bitmap = True, fname = 'aay', colrange=[-0.01,0.01]) - vis.slices.plot_3(mesh, mesh.aa[2], title = r'$A_z$', bitmap = True, fname = 'aaz', colrange=[-0.01,0.01]) - vis.slices.plot_3(mesh, mesh.accretion, title = r'$Accretion$', bitmap = True, fname = 'accretion', colrange=[0.0,0.000001]) - vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot', colrange=[0.0,1.0e-4]) - vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True) - vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True) - vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz', colrange=[-1.0e-4, 1.0e-4])#, bfieldlines=True) + if 'sym' in sys.argv: + rlnrho = [np.amin(mesh.lnrho), np.amax(mesh.lnrho)] + rrho = [ np.exp(rlnrho[0]), np.exp(rlnrho[1])] + rNcol = [ 0.0, 1.0] + ruu_tot = [ np.amin(uu_tot), np.amax(uu_tot)] + maxucomp = np.amax([np.amax(np.abs(mesh.uu[0])), np.amax(np.abs(mesh.uu[1])), np.amax(np.abs(mesh.uu[2]))]) + maxacomp = np.amax([np.amax(np.abs(mesh.aa[0])), np.amax(np.abs(mesh.aa[1])), np.amax(np.abs(mesh.aa[2]))]) + maxbcomp = np.amax([np.amax(np.abs(mesh.bb[0])), np.amax(np.abs(mesh.bb[1])), np.amax(np.abs(mesh.bb[2]))]) + ruu_xyz = [-maxucomp, maxucomp] + raa_tot = [ np.amin(aa_tot), np.amax(aa_tot)] + raa_xyz = [-maxacomp, maxacomp] + rbb_tot = [ np.amin(bb_tot), np.amax(bb_tot)] + rbb_xyz = [-maxbcomp, maxbcomp] + + if ('lim' in sys.argv) or ('sym' in sys.argv): + vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho', colrange=rlnrho) + vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho', colrange=rrho) + vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\mathrm{col}$', bitmap = True, fname = 'colden', slicetype = 'sum', colrange=rNcol) + vis.slices.plot_3(mesh, uu_tot, title = r'$|u|$', bitmap = True, fname = 'uutot', colrange=ruu_tot) + vis.slices.plot_3(mesh, mesh.uu[0], title = r'$u_x$', bitmap = True, fname = 'uux', colrange=ruu_xyz) + vis.slices.plot_3(mesh, mesh.uu[1], title = r'$u_y$', bitmap = True, fname = 'uuy', colrange=ruu_xyz) + vis.slices.plot_3(mesh, mesh.uu[2], title = r'$u_z$', bitmap = True, fname = 'uuz', colrange=ruu_xyz) + vis.slices.plot_3(mesh, aa_tot, title = r'$\|A\|$', bitmap = True, fname = 'aatot', colrange=raa_tot) + vis.slices.plot_3(mesh, mesh.aa[0], title = r'$A_x$', bitmap = True, fname = 'aax', colrange=raa_xyz) + vis.slices.plot_3(mesh, mesh.aa[1], title = r'$A_y$', bitmap = True, fname = 'aay', colrange=raa_xyz) + vis.slices.plot_3(mesh, mesh.aa[2], title = r'$A_z$', bitmap = True, fname = 'aaz', colrange=raa_xyz) + #vis.slices.plot_3(mesh, mesh.accretion, title = r'$Accretion$', bitmap = True, fname = 'accretion', colrange=[0.0,0.000001]) + vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot', colrange=rbb_tot, trimghost=3) + vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx', colrange=rbb_xyz, trimghost=3)#, bfieldlines=True) + vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby', colrange=rbb_xyz, trimghost=3)#, bfieldlines=True) + vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz', colrange=rbb_xyz, trimghost=3)#, bfieldlines=True) else: vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\ln \rho$', bitmap = True, fname = 'lnrho') vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\rho$', bitmap = True, fname = 'rho') @@ -209,10 +243,10 @@ if 'sl' in sys.argv: vis.slices.plot_3(mesh, mesh.aa[0], title = r'$A_x$', bitmap = True, fname = 'aax') vis.slices.plot_3(mesh, mesh.aa[1], title = r'$A_y$', bitmap = True, fname = 'aay') vis.slices.plot_3(mesh, mesh.aa[2], title = r'$A_z$', bitmap = True, fname = 'aaz') - vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot') - vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx')#, bfieldlines=True) - vis.slices.plot_3(mesh, mesh.bb[1], title = r'$B_y$', bitmap = True, fname = 'bby')#, bfieldlines=True) - vis.slices.plot_3(mesh, mesh.bb[2], title = r'$B_z$', bitmap = True, fname = 'bbz')#, bfieldlines=True) + vis.slices.plot_3(mesh, bb_tot, title = r'$\|B\|$', bitmap = True, fname = 'bbtot', trimghost=3) + vis.slices.plot_3(mesh, mesh.bb[0], title = r'$B_x$', bitmap = True, fname = 'bbx', trimghost=3)#, bfieldlines=True) + 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 'ts' in sys.argv: @@ -221,3 +255,46 @@ if 'ts' in sys.argv: #vis.lineplot.plot_ts(ts, isotherm=True) +if 'getvtk' in sys.argv: + mesh_file_numbers = ad.read.parse_directory(meshdir) + print(mesh_file_numbers) + maxfiles = np.amax(mesh_file_numbers) + + if os.path.exists("grouped.csv"): + df_archive = pd.read_csv("grouped.csv") + print(df_archive) + useBeq = True + else: + print("reduced.csv missing!") + useBeq = False + + + #for i in mesh_file_numbers[-1:]: + for i in mesh_file_numbers: + mesh = ad.read.Mesh(i, fdir=meshdir) + resolution = mesh.minfo.contents['AC_nx' ] + eta = mesh.minfo.contents['AC_eta'] + relhel = mesh.minfo.contents['AC_relhel'] + kk = (mesh.minfo.contents['AC_kmax']+mesh.minfo.contents['AC_kmin'])/2.0 + + if i == mesh_file_numbers[0]: + if useBeq: + #MV: Do not use unless you know what you are doing. + df_archive = df_archive.loc[df_archive['relhel'] == relhel] + df_archive = df_archive.loc[df_archive['eta'] == eta] + df_archive = df_archive.loc[df_archive['resolution'] == resolution] + df_archive = df_archive.loc[df_archive['kk'] == kk] + df_archive = df_archive.reset_index() + print(df_archive) + uu_eq = df_archive['urms_growth'].values[0] + print(uu_eq) + myBeq = np.sqrt(1.0*1.0)*uu_eq + print(myBeq) + else: + myBeq = 1.0 + + + print(" %i / %i" % (i, maxfiles)) + if mesh.ok: + #mesh.Bfield() + mesh.export_vtk_ascii(Beq = myBeq) diff --git a/include/astaroth.h b/include/astaroth.h index c73879d..d32b367 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -238,6 +238,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 + * * 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); + /** Stores a subset of the mesh stored across the devices visible to the caller back to host memory. */ AcResult acStoreWithOffset(const int3 dst, const size_t num_vertices, AcMesh* host_mesh); @@ -436,6 +441,11 @@ AcResult acNodeReduceScal(const Node node, const Stream stream, const ReductionT AcResult acNodeReduceVec(const Node node, const Stream stream_type, const ReductionType rtype, const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1, const VertexBufferHandle vtxbuf2, AcReal* result); +/** */ +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); + /* * ============================================================================= @@ -560,6 +570,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, + const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1, + const VertexBufferHandle vtxbuf2, const VertexBufferHandle vtxbuf3, AcReal* result); +/** */ AcResult acDeviceRunMPITest(void); /* diff --git a/samples/standalone/model/host_timestep.cc b/samples/standalone/model/host_timestep.cc index fb1b52c..897adb2 100644 --- a/samples/standalone/model/host_timestep.cc +++ b/samples/standalone/model/host_timestep.cc @@ -31,7 +31,7 @@ static AcReal timescale = AcReal(1.0); AcReal -host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info) +host_timestep(const AcReal& umax, const AcReal& vAmax, const AcMeshInfo& mesh_info) { const long double cdt = mesh_info.real_params[AC_cdt]; const long double cdtv = mesh_info.real_params[AC_cdtv]; @@ -40,7 +40,7 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info) const long double nu_visc = mesh_info.real_params[AC_nu_visc]; const long double eta = mesh_info.real_params[AC_eta]; const long double chi = 0; // mesh_info.real_params[AC_chi]; // TODO not calculated - const long double gamma = mesh_info.real_params[AC_gamma]; + const long double gamma = mesh_info.real_params[AC_gamma]; //TODO this does not make sense here at all. const long double dsmin = mesh_info.real_params[AC_dsmin]; // Old ones from legacy Astaroth @@ -49,12 +49,11 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info) // New, closer to the actual Courant timestep // See Pencil Code user manual p. 38 (timestep section) - const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l)); + //const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l)); + const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + vAmax*vAmax)); const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), - max(gamma, chi)); // + 1; // TODO NOTE: comment the +1 out to - // get scientifically accurate results - // MV: White the +1? It was messing up my computations! + max(gamma, chi)); const long double dt = min(uu_dt, visc_dt); return AcReal(timescale) * AcReal(dt); diff --git a/samples/standalone/model/host_timestep.h b/samples/standalone/model/host_timestep.h index dd5cdf7..9074b26 100644 --- a/samples/standalone/model/host_timestep.h +++ b/samples/standalone/model/host_timestep.h @@ -27,6 +27,6 @@ #pragma once #include "astaroth.h" -AcReal host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info); +AcReal host_timestep(const AcReal& umax, const AcReal& vAmax, const AcMeshInfo& mesh_info); void set_timescale(const AcReal scale); diff --git a/samples/standalone/simulation.cc b/samples/standalone/simulation.cc index 01106af..0bf193e 100644 --- a/samples/standalone/simulation.cc +++ b/samples/standalone/simulation.cc @@ -39,10 +39,16 @@ #include #include #include +#include // NEED TO BE DEFINED HERE. IS NOT NOTICED BY compile_acc call. -#define LFORCING (1) +#define LFORCING (0) #define LSINK (0) +#ifdef BFIELDX +#define LBFIELD (1) +#else +#define LBFIELD (0) +#endif // Write all setting info into a separate ascii file. This is done to guarantee // that we have the data specifi information in the thing, even though in @@ -95,6 +101,7 @@ write_mesh_info(const AcMeshInfo* config) fclose(infotxt); } + // This funtion writes a run state into a set of C binaries. static inline void save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step) @@ -134,6 +141,61 @@ save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step) } } +/* +// This funtion writes a run state into a set of C binaries. +static inline void +save_slice_cut(const AcMesh& save_mesh, const int step, const AcReal t_step) +{ + FILE* save_ptr; + + for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) { + const size_t n = acVertexBufferSize(save_mesh.info); + + const char* buffername = vtxbuf_names[w]; + char cstep[11]; + char bin_filename_xy[80] = "\0"; + char bin_filename_xz[80] = "\0"; + char bin_filename_yz[80] = "\0"; + + // sprintf(bin_filename, ""); + + sprintf(cstep, "%d", step); + + strcat(bin_filename_xy, buffername); + strcat(bin_filename_xy, "_"); + strcat(bin_filename_xy, cstep); + strcat(bin_filename_xy, ".mxy"); + + strcat(bin_filename_xz, buffername); + strcat(bin_filename_xz, "_"); + strcat(bin_filename_xz, cstep); + strcat(bin_filename_xz, ".mxz"); + + strcat(bin_filename_yz, buffername); + strcat(bin_filename_yz, "_"); + strcat(bin_filename_yz, cstep); + strcat(bin_filename_yz, ".myz"); + + printf("Slice files %s, %s, %s, \n", + bin_filename_xy, bin_filename_xz, bin_filename_yz); + + save_ptr = fopen(bin_filename, "wb"); + + // Start file with time stamp + AcReal write_long_buf = (AcReal)t_step; + fwrite(&write_long_buf, sizeof(AcReal), 1, save_ptr); + // Grid data + for (size_t i = 0; i < n; ++i) { + const AcReal point_val = save_mesh.vertex_buffer[VertexBufferHandle(w)][i]; + AcReal write_long_buf2 = (AcReal)point_val; + fwrite(&write_long_buf2, sizeof(AcReal), 1, save_ptr); + } + fclose(save_ptr); + } +} +*/ + + // This funtion reads a run state from a set of C binaries. static inline void read_mesh(AcMesh& read_mesh, const int step, AcReal* t_step) @@ -181,7 +243,7 @@ read_mesh(AcMesh& read_mesh, const int step, AcReal* t_step) // appends an ascii file to contain all the result. static inline void print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file, - const AcReal sink_mass, const AcReal accreted_mass) + const AcReal sink_mass, const AcReal accreted_mass, int* found_nan) { AcReal buf_rms, buf_max, buf_min; @@ -200,6 +262,24 @@ print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* di fprintf(diag_file, "%d %e %e %e %e %e ", step, double(t_step), double(dt), double(buf_min), double(buf_rms), double(buf_max)); +#if LBFIELD + buf_max = acReduceVec(RTYPE_MAX, BFIELDX, BFIELDY, BFIELDZ); + buf_min = acReduceVec(RTYPE_MIN, BFIELDX, BFIELDY, BFIELDZ); + buf_rms = acReduceVec(RTYPE_RMS, BFIELDX, BFIELDY, BFIELDZ); + + printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, "bb total", double(buf_min), + double(buf_rms), double(buf_max)); + fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max)); + + buf_max = acReduceVecScal(RTYPE_ALFVEN_MAX, BFIELDX, BFIELDY, BFIELDZ, VTXBUF_LNRHO); + buf_min = acReduceVecScal(RTYPE_ALFVEN_MIN, BFIELDX, BFIELDY, BFIELDZ, VTXBUF_LNRHO); + buf_rms = acReduceVecScal(RTYPE_ALFVEN_RMS, BFIELDX, BFIELDY, BFIELDZ, VTXBUF_LNRHO); + + printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, "vA total", double(buf_min), + double(buf_rms), double(buf_max)); + fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max)); +#endif + // Calculate rms, min and max from the variables as scalars for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) { buf_max = acReduceScal(RTYPE_MAX, VertexBufferHandle(i)); @@ -209,6 +289,11 @@ print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* di printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, vtxbuf_names[i], double(buf_min), double(buf_rms), double(buf_max)); fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max)); + + if (isnan(buf_max) || isnan(buf_min) || isnan(buf_rms)) { + *found_nan = 1; + } + } if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) { @@ -252,11 +337,16 @@ run_simulation(const char* config_path) acLoad(*mesh); FILE* diag_file; + int found_nan = 0, found_stop = 0; // Nan or inf finder to give an error signal diag_file = fopen("timeseries.ts", "a"); // Generate the title row. if (start_step == 0) { fprintf(diag_file, "step t_step dt uu_total_min uu_total_rms uu_total_max "); +#if LBFIELD + fprintf(diag_file, "bb_total_min bb_total_rms bb_total_max "); + fprintf(diag_file, "vA_total_min vA_total_rms vA_total_max "); +#endif for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) { fprintf(diag_file, "%s_min %s_rms %s_max ", vtxbuf_names[i], vtxbuf_names[i], vtxbuf_names[i]); @@ -271,9 +361,9 @@ run_simulation(const char* config_path) if (start_step == 0) { #if LSINK - print_diagnostics(0, AcReal(.0), t_step, diag_file, mesh_info.real_params[AC_M_sink_init], 0.0); + print_diagnostics(0, AcReal(.0), t_step, diag_file, mesh_info.real_params[AC_M_sink_init], 0.0, &found_nan); #else - print_diagnostics(0, AcReal(.0), t_step, diag_file, -1.0, -1.0); + print_diagnostics(0, AcReal(.0), t_step, diag_file, -1.0, -1.0, &found_nan); #endif } @@ -297,9 +387,17 @@ run_simulation(const char* config_path) /* Step the simulation */ AcReal accreted_mass = 0.0; AcReal sink_mass = 0.0; + AcReal dt_typical = 0.0; + int dtcounter = 0; for (int i = start_step + 1; i < max_steps; ++i) { const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); - const AcReal dt = host_timestep(umax, mesh_info); +#if LBFIELD + const AcReal vAmax = acReduceVecScal(RTYPE_ALFVEN_MAX, BFIELDX, BFIELDY, BFIELDZ, VTXBUF_LNRHO); + const AcReal uref = max(umax, vAmax); + const AcReal dt = host_timestep(uref, vAmax, mesh_info); +#else + const AcReal dt = host_timestep(umax, 0.0l, mesh_info); +#endif #if LSINK @@ -332,7 +430,11 @@ run_simulation(const char* config_path) t_step += dt; - + /* Get the sense of a typical timestep */ + if (i < start_step + 100) { + dt_typical = dt; + } + /* Save the simulation state and print diagnostics */ if ((i % save_steps) == 0) { @@ -343,7 +445,7 @@ run_simulation(const char* config_path) timeseries.ts. */ - print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass); + print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass, &found_nan); #if LSINK printf("sink mass is: %.15e \n", double(sink_mass)); printf("accreted mass is: %.15e \n", double(accreted_mass)); @@ -387,6 +489,45 @@ run_simulation(const char* config_path) } } + // End loop if dt is too low + if (dt < dt_typical/AcReal(1e5)) { + if (dtcounter > 10) { + printf("dt = %e TOO LOW! Ending run at t = %#e \n", double(dt), double(t_step)); + acBoundcondStep(); + acStore(mesh); + save_mesh(*mesh, i, t_step); + break; + } else { + dtcounter += 1; + } + } else { + dtcounter = 0; + } + + // End loop if nan is found + if (found_nan > 0) { + printf("Found nan at t = %e \n", double(t_step)); + acBoundcondStep(); + acStore(mesh); + save_mesh(*mesh, i, t_step); + break; + } + + // End loop if STOP file is found + if( access("STOP", F_OK ) != -1 ) { + found_stop = 1; + } else { + found_stop = 0; + } + + if (found_stop == 1) { + printf("Found STOP file at t = %e \n", double(t_step)); + acBoundcondStep(); + acStore(mesh); + save_mesh(*mesh, i, t_step); + break; + } + } //////Save the final snapshot diff --git a/src/core/astaroth.cc b/src/core/astaroth.cc index 3a38240..1c8d0f7 100644 --- a/src/core/astaroth.cc +++ b/src/core/astaroth.cc @@ -126,6 +126,15 @@ acReduceVec(const ReductionType rtype, const VertexBufferHandle a, const VertexB return result; } +AcReal +acReduceVecScal(const ReductionType rtype, const VertexBufferHandle a, const VertexBufferHandle b, + const VertexBufferHandle c, const VertexBufferHandle d) +{ + AcReal result; + acNodeReduceVecScal(nodes[0], STREAM_DEFAULT, rtype, a, b, c, d, &result); + return result; +} + AcResult acStoreWithOffset(const int3 dst, const size_t num_vertices, AcMesh* host_mesh) { diff --git a/src/core/device.cc b/src/core/device.cc index 3e5515a..a6ec793 100644 --- a/src/core/device.cc +++ b/src/core/device.cc @@ -474,6 +474,28 @@ acDeviceReduceVec(const Device device, const Stream stream, const ReductionType return AC_SUCCESS; } +AcResult +acDeviceReduceVecScal(const Device device, const Stream stream, const ReductionType rtype, + const VertexBufferHandle vtxbuf0, const VertexBufferHandle vtxbuf1, + const VertexBufferHandle vtxbuf2, const VertexBufferHandle vtxbuf3, AcReal* result) +{ + cudaSetDevice(device->id); + + const int3 start = (int3){device->local_config.int_params[AC_nx_min], + device->local_config.int_params[AC_ny_min], + device->local_config.int_params[AC_nz_min]}; + + const int3 end = (int3){device->local_config.int_params[AC_nx_max], + device->local_config.int_params[AC_ny_max], + device->local_config.int_params[AC_nz_max]}; + + *result = acKernelReduceVecScal(device->streams[stream], rtype, start, end, device->vba.in[vtxbuf0], + device->vba.in[vtxbuf1], device->vba.in[vtxbuf2], device->vba.in[vtxbuf3], + device->reduce_scratchpad, device->reduce_result); + return AC_SUCCESS; +} + + #if AC_MPI_ENABLED /** Quick overview of the MPI implementation: @@ -2047,4 +2069,7 @@ acGridReduceVec(const Stream stream, const ReductionType rtype, const VertexBuff return acMPIReduceScal(local_result, rtype, result); } + +//MV: for MPI we will need acGridReduceVecScal() to get Alfven speeds etc. TODO + #endif // AC_MPI_ENABLED diff --git a/src/core/kernels/kernels.h b/src/core/kernels/kernels.h index 577567b..fc9c745 100644 --- a/src/core/kernels/kernels.h +++ b/src/core/kernels/kernels.h @@ -72,6 +72,12 @@ AcReal acKernelReduceVec(const cudaStream_t stream, const ReductionType rtype, c const int3 end, const AcReal* vtxbuf0, const AcReal* vtxbuf1, const AcReal* vtxbuf2, AcReal* scratchpad, AcReal* reduce_result); +/** */ +AcReal acKernelReduceVecScal(const cudaStream_t stream, const ReductionType rtype, const int3 start, + const int3 end, const AcReal* vtxbuf0, const AcReal* vtxbuf1, + const AcReal* vtxbuf2, const AcReal* vtxbuf3, AcReal* scratchpad, AcReal* reduce_result); + + #ifdef __cplusplus } // extern "C" #endif diff --git a/src/core/kernels/reductions.cuh b/src/core/kernels/reductions.cuh index 1f40df3..ba9e823 100644 --- a/src/core/kernels/reductions.cuh +++ b/src/core/kernels/reductions.cuh @@ -10,6 +10,7 @@ Reduction steps: // Function pointer definitions typedef AcReal (*FilterFunc)(const AcReal&); typedef AcReal (*FilterFuncVec)(const AcReal&, const AcReal&, const AcReal&); +typedef AcReal (*FilterFuncVecScal)(const AcReal&, const AcReal&, const AcReal&, const AcReal&); typedef AcReal (*ReduceFunc)(const AcReal&, const AcReal&); // clang-format off @@ -41,6 +42,13 @@ dsquared_vec(const AcReal& a, const AcReal& b, const AcReal& c) { return dsquare static __device__ inline AcReal dexp_squared_vec(const AcReal& a, const AcReal& b, const AcReal& c) { return dexp_squared(a) + dexp_squared(b) + dexp_squared(c); } + +static __device__ inline AcReal +dlength_alf(const AcReal& a, const AcReal& b, const AcReal& c, const AcReal& d) { return sqrt(a*a + b*b + c*c)/sqrt(AcReal(4.0)*M_PI*exp(d)); } + +static __device__ inline AcReal +dsquared_alf(const AcReal& a, const AcReal& b, const AcReal& c, const AcReal& d) { return (dsquared(a) + dsquared(b) + dsquared(c))/(AcReal(4.0)*M_PI*exp(d)); } + // clang-format on #include @@ -97,6 +105,35 @@ kernel_filter_vec(const __restrict__ AcReal* src0, const __restrict__ AcReal* sr src2[IDX(src_idx)]); } +template +static __global__ void +kernel_filter_vec_scal(const __restrict__ AcReal* src0, const __restrict__ AcReal* src1, + const __restrict__ AcReal* src2, const __restrict__ AcReal* src3, + const int3 start, const int3 end, AcReal* dst) +{ + const int3 src_idx = (int3){start.x + threadIdx.x + blockIdx.x * blockDim.x, + start.y + threadIdx.y + blockIdx.y * blockDim.y, + start.z + threadIdx.z + blockIdx.z * blockDim.z}; + + const int nx = end.x - start.x; + const int ny = end.y - start.y; + const int nz = end.z - start.z; + (void)nz; // Suppressed unused variable warning when not compiling with debug flags + + const int3 dst_idx = (int3){threadIdx.x + blockIdx.x * blockDim.x, + threadIdx.y + blockIdx.y * blockDim.y, + threadIdx.z + blockIdx.z * blockDim.z}; + + assert(src_idx.x < DCONST(AC_nx_max) && src_idx.y < DCONST(AC_ny_max) && + src_idx.z < DCONST(AC_nz_max)); + assert(dst_idx.x < nx && dst_idx.y < ny && dst_idx.z < nz); + assert(dst_idx.x + dst_idx.y * nx + dst_idx.z * nx * ny < nx * ny * nz); + + dst[dst_idx.x + dst_idx.y * nx + dst_idx.z * nx * ny] = filter( + src0[IDX(src_idx)], src1[IDX(src_idx)], src2[IDX(src_idx)], src3[IDX(src_idx)]); +} + + template static __global__ void kernel_reduce(AcReal* scratchpad, const int num_elems) @@ -253,3 +290,54 @@ acKernelReduceVec(const cudaStream_t stream, const ReductionType rtype, const in cudaMemcpy(&result, reduce_result, sizeof(AcReal), cudaMemcpyDeviceToHost); return result; } + +AcReal +acKernelReduceVecScal(const cudaStream_t stream, const ReductionType rtype, const int3 start, + const int3 end, const AcReal* vtxbuf0, const AcReal* vtxbuf1, + const AcReal* vtxbuf2, const AcReal* vtxbuf3, AcReal* scratchpad, AcReal* reduce_result) +{ + const unsigned nx = end.x - start.x; + const unsigned ny = end.y - start.y; + const unsigned nz = end.z - start.z; + const unsigned num_elems = nx * ny * nz; + + const dim3 tpb_filter(32, 4, 1); + const dim3 bpg_filter((unsigned int)ceil(nx / AcReal(tpb_filter.x)), + (unsigned int)ceil(ny / AcReal(tpb_filter.y)), + (unsigned int)ceil(nz / AcReal(tpb_filter.z))); + + const int tpb_reduce = 128; + const int bpg_reduce = num_elems / tpb_reduce; + + ERRCHK(nx >= tpb_filter.x); + ERRCHK(ny >= tpb_filter.y); + ERRCHK(nz >= tpb_filter.z); + ERRCHK(tpb_reduce <= num_elems); + ERRCHK(nx * ny * nz % 2 == 0); + + //NOTE: currently this has been made to only calculate afven speeds from the diagnostics. + + // clang-format off + if (rtype == RTYPE_ALFVEN_MAX) { + kernel_filter_vec_scal<<>>(vtxbuf0, vtxbuf1, vtxbuf2, vtxbuf3, start, end, scratchpad); + kernel_reduce<<>>(scratchpad, num_elems); + kernel_reduce_block<<<1, 1, 0, stream>>>(scratchpad, bpg_reduce, tpb_reduce, reduce_result); + } else if (rtype == RTYPE_ALFVEN_MIN) { + kernel_filter_vec_scal<<>>(vtxbuf0, vtxbuf1, vtxbuf2, vtxbuf3, start, end, scratchpad); + kernel_reduce<<>>(scratchpad, num_elems); + kernel_reduce_block<<<1, 1, 0, stream>>>(scratchpad, bpg_reduce, tpb_reduce, reduce_result); + } else if (rtype == RTYPE_ALFVEN_RMS) { + kernel_filter_vec_scal<<>>(vtxbuf0, vtxbuf1, vtxbuf2, vtxbuf3, start, end, scratchpad); + kernel_reduce<<>>(scratchpad, num_elems); + kernel_reduce_block<<<1, 1, 0, stream>>>(scratchpad, bpg_reduce, tpb_reduce, reduce_result); + } else { + ERROR("Unrecognized rtype"); + } + // clang-format on + + cudaStreamSynchronize(stream); + AcReal result; + cudaMemcpy(&result, reduce_result, sizeof(AcReal), cudaMemcpyDeviceToHost); + return result; +} + diff --git a/src/core/node.cc b/src/core/node.cc index 7a0ca03..70c3115 100644 --- a/src/core/node.cc +++ b/src/core/node.cc @@ -797,13 +797,13 @@ simple_final_reduce_scal(const Node node, const ReductionType& rtype, const AcRe { AcReal res = results[0]; for (int i = 1; i < n; ++i) { - if (rtype == RTYPE_MAX) { + if (rtype == RTYPE_MAX || rtype == RTYPE_ALFVEN_MAX) { res = max(res, results[i]); } - else if (rtype == RTYPE_MIN) { + else if (rtype == RTYPE_MIN || rtype == RTYPE_ALFVEN_MIN) { res = min(res, results[i]); } - else if (rtype == RTYPE_RMS || rtype == RTYPE_RMS_EXP || rtype == RTYPE_SUM) { + else if (rtype == RTYPE_RMS || rtype == RTYPE_RMS_EXP || rtype == RTYPE_SUM || rtype == RTYPE_ALFVEN_RMS) { res = sum(res, results[i]); } else { @@ -811,7 +811,7 @@ simple_final_reduce_scal(const Node node, const ReductionType& rtype, const AcRe } } - if (rtype == RTYPE_RMS || rtype == RTYPE_RMS_EXP) { + if (rtype == RTYPE_RMS || rtype == RTYPE_RMS_EXP || rtype == RTYPE_ALFVEN_RMS) { const AcReal inv_n = AcReal(1.) / (node->grid.n.x * node->grid.n.y * node->grid.n.z); res = sqrt(inv_n * res); } @@ -850,3 +850,21 @@ acNodeReduceVec(const Node node, const Stream stream, const ReductionType rtype, *result = simple_final_reduce_scal(node, rtype, results, node->num_devices); return AC_SUCCESS; } + +AcResult +acNodeReduceVecScal(const Node node, const Stream stream, const ReductionType rtype, + const VertexBufferHandle a, const VertexBufferHandle b, const VertexBufferHandle c, + const VertexBufferHandle d, AcReal* result) +{ + acNodeSynchronizeStream(node, STREAM_ALL); + + AcReal results[node->num_devices]; + // #pragma omp parallel for + for (int i = 0; i < node->num_devices; ++i) { + acDeviceReduceVecScal(node->devices[i], stream, rtype, a, b, c, d, &results[i]); + } + + *result = simple_final_reduce_scal(node, rtype, results, node->num_devices); + return AC_SUCCESS; +} + diff --git a/src/utils/modelsolver.c b/src/utils/modelsolver.c index 4117168..3347d75 100644 --- a/src/utils/modelsolver.c +++ b/src/utils/modelsolver.c @@ -702,17 +702,19 @@ momentum(const VectorData uu, const ScalarData lnrho static inline Vector induction(const VectorData uu, const VectorData aa) { - Vector ind; // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla // x A)) in order to avoid taking the first derivative twice (did the math, // yes this actually works. See pg.28 in arXiv:astro-ph/0109497) - // u cross B - ETA * mu0 * (mu0^-1 * [- laplace A + grad div A ]) + // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ]) const Vector B = curl(aa); - const Vector grad_div = gradient_of_divergence(aa); + //MV: Due to gauge freedom we can reduce the gradient of scalar (divergence) from the equation + //const Vector grad_div = gradient_of_divergence(aa); const Vector lap = laplace_vec(aa); - // Note, mu0 is cancelled out - ind = cross(vecvalue(uu), B) - getReal(AC_eta) * (grad_div - lap); + // Note, AC_mu0 is cancelled out + //MV: Due to gauge freedom we can reduce the gradient of scalar (divergence) from the equation + //const Vector ind = cross(value(uu), B) - getReal(AC_eta) * (grad_div - lap); + const Vector ind = cross(vecvalue(uu), B) + getReal(AC_eta) * lap; return ind; }