diff --git a/examples/mhs_atmosphere/Bz_Heliocentric.fits b/examples/mhs_atmosphere/Bz_Heliocentric.fits new file mode 100755 index 0000000..7539bd3 Binary files /dev/null and b/examples/mhs_atmosphere/Bz_Heliocentric.fits differ diff --git a/examples/mhs_atmosphere/hmi_atmosphere.py b/examples/mhs_atmosphere/hmi_atmosphere.py index 20dbfb6..47bdf15 100644 --- a/examples/mhs_atmosphere/hmi_atmosphere.py +++ b/examples/mhs_atmosphere/hmi_atmosphere.py @@ -31,6 +31,10 @@ import pysac.mhs_atmosphere as atm import astropy.units as u from pysac.mhs_atmosphere.parameters.model_pars import hmi_model as model_pars +import astropy +l_astropy_0=True +if astropy.__version__[0]=='1': + l_astropy_0=False #============================================================================== #check whether mpi is required and the number of procs = size #============================================================================== @@ -41,9 +45,11 @@ rank = comm.Get_rank() size = comm.Get_size() + collective=True l_mpi = True l_mpi = l_mpi and (size != 1) except ImportError: + collective=False l_mpi = False rank = 0 size = 1 @@ -57,11 +63,29 @@ scales, physical_constants = \ atm.units_const.get_parameters() +if not l_astropy_0: + dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram_fits' +else: + dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram.fits' +l_newdata = True # change this to False if a local copy already exists at ~/sunpy/data/ #obtain code coordinates and model parameters in astropy units -coords = atm.model_pars.get_hmi_map(model_pars, option_pars, indx = [1787,1798,1818,1822], - dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram_fits', - l_newdata = False +model_pars['Nxyz'] = [64,64,64] # 3D grid +model_pars['xyz'] = [-0.63*u.Mm,0.63*u.Mm,-0.63*u.Mm,0.63*u.Mm,0.0*u.Mm,3.8*u.Mm] +indx = [1785,1800,1814,1829] +interpfactor=5 +frac = [0.25,0.25] +coords = atm.model_pars.get_hmi_coords( + model_pars['Nxyz'], + model_pars['xyz'], + indx = indx, + dataset = dataset, + l_newdata = False, + interpfactor=interpfactor, + frac=frac, + rank=rank, + lmpi=l_mpi ) +model_pars['xyz'][0:4] = coords['xmin'],coords['xmax'],coords['ymin'],coords['ymax'] #interpolate the hs 1D profiles from empirical data source[s] empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu']) @@ -92,12 +116,19 @@ # load flux tube footpoint parameters #============================================================================== # axial location and value of Bz at each footpoint -model_pars['B_corona']/=model_pars['nftubes'] -xi, yi, Si = atm.flux_tubes.get_flux_tubes( - model_pars, - coords, - option_pars - ) +Stmp,xtmp,ytmp,FWHM,sdummy,xdummy,ydummy,cmax,cmin \ + =atm.parameters.model_pars.get_hmi_map( + indx, + dataset = dataset, + l_newdata = False, + interpfactor=interpfactor, + frac=frac, + rank=rank, + lmpi=l_mpi + ) +xi, yi, Si = xtmp.to(u.Mm).reshape(xtmp.size),\ + ytmp.to(u.Mm).reshape(ytmp.size), Stmp.reshape(Stmp.size) +model_pars['radial_scale'] = 0.5*FWHM.to(u.Mm)/np.sqrt(2*np.log(2)) #============================================================================== # split domain into processes if mpi #============================================================================== @@ -105,17 +136,21 @@ coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1], coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]] +axindex=np.arange(model_pars['Nxyz'][0]) # split the grid between processes for mpi if l_mpi: x_chunks = np.array_split(ax, size, axis=0) y_chunks = np.array_split(ay, size, axis=0) z_chunks = np.array_split(az, size, axis=0) + i_chunks = np.array_split(axindex, size, axis=0) x = comm.scatter(x_chunks, root=0) y = comm.scatter(y_chunks, root=0) z = comm.scatter(z_chunks, root=0) + xindex=i_chunks[rank] else: x, y, z = ax, ay, az + xindex=axindex x = u.Quantity(x, unit=coords['xmin'].unit) y = u.Quantity(y, unit=coords['ymin'].unit) @@ -137,8 +172,8 @@ #============================================================================== #calculate the magnetic field and pressure/density balancing expressions #============================================================================== -for i in range(0,model_pars['nftubes']): - for j in range(i,model_pars['nftubes']): +for i in range(0,Si.size): + for j in range(i,Si.size): if rank == 0: print'calculating ij-pair:',i,j if i == j: @@ -175,6 +210,19 @@ # clear some memory del pressure_mi, rho_mi, Bxi, Byi ,Bzi, B2x, B2y + +import matplotlib.pyplot as plt +plt.figure() +plt.pcolormesh(x[:,:,0].T.value,y[:,:,0].T.value,Bz[:,:,0].T.value,vmin=cmin,vmax=cmax) +plt.xlabel('lon [Mm]') +plt.ylabel('lat [Mm]') +plt.axis([x.min().value,x.max().value,y.min().value,y.max().value]) +cbar = plt.colorbar() +cbar.ax.set_ylabel(r'$B_z$ [T]') +cbar.solids.set_edgecolor("face") +plt.savefig('model_derived_hmi.png') +plt.close() +print 'Bz[:,:,0].min(), Bz[:,:,0].max()=', Bz[:,:,0].min(), Bz[:,:,0].max() #============================================================================== # Construct 3D hs arrays and then add the mhs adjustments to obtain atmosphere #============================================================================== @@ -192,11 +240,27 @@ ) magp = (Bx**2 + By**2 + Bz**2)/(2.*physical_constants['mu0']) if rank ==0: - print'max B corona = ',magp[:,:,-1].max().decompose() - print'min B corona = ',magp[:,:,-1].min().decompose() + print'max pB corona = ',magp[:,:,-1].max().decompose() + print'min pB corona = ',magp[:,:,-1].min().decompose() energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants) +Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) +Rgas[:] = Rgas_z +temperature = pressure/rho/Rgas +if not option_pars['l_hdonly']: + inan = np.where(magp <=1e-7*pressure.min()) + magpbeta = magp + magpbeta[inan] = 1e-7*pressure.min() # low pressure floor to avoid NaN + pbeta = pressure/magpbeta +else: + pbeta = magp+1.0 #dummy to avoid NaN +alfven = np.sqrt(2.*magp/rho) +#if rank == 0: +# print'Alfven speed Z.min to Z.max =',\ +# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2, 0].decompose(),\ +# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2,-1].decompose() +cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) #============================================================================ # Save data for SAC and plotting #============================================================================ @@ -212,66 +276,67 @@ # save the variables for the initialisation of a SAC simulation atm.mhs_snapshot.save_SACvariables( - filename, - rho, - Bx, - By, - Bz, - energy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + filename, + rho, + Bx, + By, + Bz, + energy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) # save the balancing forces as the background source terms for SAC simulation atm.mhs_snapshot.save_SACsources( - sourcefile, - Fx, - Fy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + sourcefile, + Fx, + Fy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) # save auxilliary variable and 1D profiles for plotting and analysis -Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) -Rgas[:] = Rgas_z -temperature = pressure/rho/Rgas -if not option_pars['l_hdonly']: - inan = np.where(magp <=1e-7*pressure.min()) - magpbeta = magp - magpbeta[inan] = 1e-7*pressure.min() # low pressure floor to avoid NaN - pbeta = pressure/magpbeta -else: - pbeta = magp+1.0 #dummy to avoid NaN -alfven = np.sqrt(2.*magp/rho) -#if rank == 0: -# print'Alfven speed Z.min to Z.max =',\ -# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2, 0].decompose(),\ -# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2,-1].decompose() -cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) atm.mhs_snapshot.save_auxilliary3D( - aux3D, - pressure_m, - rho_m, - temperature, - pbeta, - alfven, - cspeed, - Btensx, - Btensy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) -atm.mhs_snapshot.save_auxilliary1D( - aux1D, - pressure_Z, - rho_Z, - Rgas_Z, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + aux3D, + pressure_m, + rho_m, + temperature, + pbeta, + alfven, + cspeed, + Btensx, + Btensy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) +if rank==0: + atm.mhs_snapshot.save_auxilliary1D( + aux1D, + pressure_Z, + rho_Z, + Rgas_Z, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + rank=rank, + collective=False + ) +if rho.min()<0 or pressure.min()<0: + print"FAIL rank {}: negative rho.min() {} and/or pressure.min() {}.".format( + rank,rho.min(),pressure.min()) +FWHM = 2*np.sqrt(np.log(2))*model_pars['radial_scale'] +print'FWHM(0) =',FWHM diff --git a/examples/mhs_atmosphere/observed_test.py b/examples/mhs_atmosphere/observed_test.py new file mode 100755 index 0000000..51402e6 --- /dev/null +++ b/examples/mhs_atmosphere/observed_test.py @@ -0,0 +1,319 @@ +# -*- coding: utf-8 -*- +""" +Created on Monday Feb 18 2019 + +@author: BS + +Load a fits file, find magnetic sources and create 3D atmosphere. + +TO DO +1) Fix length scale and magnetic field for observation +2) Check boundary +3) Tweak parameters to create stable atmosphere +4) Check SAC stability + +""" + +import os +import numpy as np +import pysac.mhs_atmosphere as atm +import astropy.units as u +from pysac.mhs_atmosphere.parameters.model_pars import obs_test as model_pars #mfe_setup + +import pysac.mhs_atmosphere.obs_test + + +#============================================================================== +#check whether mpi is required and the number of procs = size +#============================================================================== +try: + from mpi4py import MPI + + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + size = comm.Get_size() + + collective=True + l_mpi = True + l_mpi = l_mpi and (size != 1) + print 'rank = ',rank +except ImportError: + collective=False + l_mpi = False + rank = 0 + size = 1 +#============================================================================== +#set up model parameters +#============================================================================== +local_procs=1 +#standard set of logical switches +option_pars = atm.options.set_options(model_pars, l_mpi, l_gdf=True) +#standard conversion to dimensionless units and physical constants +scales, physical_constants = \ + atm.units_const.get_parameters() + +physical_constants['gravity'] = -274.0*u.m/u.s**2 + +#obtain code coordinates and model parameters in astropy units +coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz'])) +#interpolate the hs 1D profiles from empirical data source[s] +empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu']) + +table = \ + atm.hs_atmosphere.interpolate_atmosphere(empirical_data, + coords['Zext'] + ) +if model_pars['model'] == 'mfe_setup': + table['rho'] = table['rho'] + table['rho'].min()*3.6 + + +#============================================================================== +# load magnetogram and find flux tubes +#============================================================================== + +#load the fits file +obv_bz=atm.obs_test.read_fits_bz('Bz_Heliocentric.fits') + +#find flux tubes +xi,yi,Si=atm.obs_test.flat_fits(obv_bz,2) + +xi *= 0.0769230769*u.Mm #0.1 +xi -= 5.0*u.Mm +yi *= 0.0769230769*u.Mm #0.1 +yi -= 5.0*u.Mm +Si *= 0.00001*u.T #0.0001 from G to T + +print xi +print yi +print Si + + +#============================================================================== +# split domain into processes if mpi +#============================================================================== +ax, ay, az = np.mgrid[coords['xmin']:coords['xmax']:1j*model_pars['Nxyz'][0], + coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1], + coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]] + +axindex=np.arange(model_pars['Nxyz'][0]) +# split the grid between processes for mpi +if l_mpi: + x_chunks = np.array_split(ax, size, axis=0) + y_chunks = np.array_split(ay, size, axis=0) + z_chunks = np.array_split(az, size, axis=0) + i_chunks = np.array_split(axindex, size, axis=0) + + x = comm.scatter(x_chunks, root=0) + y = comm.scatter(y_chunks, root=0) + z = comm.scatter(z_chunks, root=0) + xindex=i_chunks[rank] +else: + x, y, z = ax, ay, az + xindex=axindex + +x = u.Quantity(x, unit=coords['xmin'].unit) +y = u.Quantity(y, unit=coords['ymin'].unit) +z = u.Quantity(z, unit=coords['zmin'].unit) + + +#============================================================================== +# initialize zero arrays in which to add magnetic field and mhs adjustments +#============================================================================== +Bx = u.Quantity(np.zeros(x.shape), unit=u.T) # magnetic x-component +By = u.Quantity(np.zeros(x.shape), unit=u.T) # magnetic y-component +Bz = u.Quantity(np.zeros(x.shape), unit=u.T) # magnetic z-component +pressure_m = u.Quantity(np.zeros(x.shape), unit=u.Pa) # magneto-hydrostatic adjustment to pressure +rho_m = u.Quantity(np.zeros(x.shape), unit=u.kg/u.m**3) # magneto-hydrostatic adjustment to density +# initialize zero arrays in which to add balancing forces and magnetic tension +Fx = u.Quantity(np.zeros(x.shape), unit=u.N/u.m**3) # balancing force x-component +Fy = u.Quantity(np.zeros(x.shape), unit=u.N/u.m**3) # balancing force y-component +# total tension force for comparison with residual balancing force +Btensx = u.Quantity(np.zeros(x.shape), unit=u.N/u.m**3) +Btensy = u.Quantity(np.zeros(x.shape), unit=u.N/u.m**3) +#============================================================================== +#calculate the magnetic field and pressure/density balancing expressions +#============================================================================== +for i in range(0,len(Si)): + for j in range(i,len(Si)): + if rank == 0: + print'calculating ij-pair:',i,j + if i == j: + pressure_mi, rho_mi, Bxi, Byi ,Bzi, B2x, B2y =\ + atm.flux_tubes.construct_magnetic_field( + x, y, z, + xi[i], yi[i], Si[i], + model_pars, option_pars, + physical_constants, + scales + ) + Bx, By, Bz = Bxi+Bx, Byi+By ,Bzi+Bz + Btensx += B2x + Btensy += B2y + pressure_m += pressure_mi + rho_m += rho_mi + else: + pressure_mi, rho_mi, Fxi, Fyi, B2x, B2y =\ + atm.flux_tubes.construct_pairwise_field( + x, y, z, + xi[i], yi[i], + xi[j], yi[j], Si[i], Si[j], + model_pars, + option_pars, + physical_constants, + scales + ) + pressure_m += pressure_mi + rho_m += rho_mi + Fx += Fxi + Fy += Fyi + Btensx += B2x + Btensy += B2y + +# clear some memory +del pressure_mi, rho_mi, Bxi, Byi ,Bzi, B2x, B2y + +print 'magnetic pressure has nans: ', np.isnan(pressure_m).any() +print 'magnetic density has nans: ', np.isnan(rho_m).any() + +#============================================================================== +#calculate 1d hydrostatic balance from empirical density profile +#============================================================================== +# the hs pressure balance is enhanced by pressure equivalent to the +# residual mean coronal magnetic pressure contribution once the magnetic +# field has been applied +magp_meanz = np.ones(len(coords['Z'])) * u.one +#magp_meanz *= model_pars['pBplus']**2/(2*physical_constants['mu0']) +magp_meanz *= np.abs(np.mean(pressure_m[:,:,-1]))#**2/(2*physical_constants['mu0']) + +magrho_meanz = np.ones(len(coords['Z'])) * u.one +magrho_meanz *= np.abs(np.mean(rho_m[:,:,0])) + +#print magp_meanz + +pressure_Z, rho_Z, Rgas_Z = atm.hs_atmosphere.vertical_profile_test( + coords['Z'], + table, + magp_meanz, + physical_constants, + coords['dz'], + magrho_meanz + ) + +#============================================================================== +# Construct 3D hs arrays and then add the mhs adjustments to obtain atmosphere +#============================================================================== +# select the 1D array spanning the local mpi process; the add/sub of dz to +# ensure all indices are used, but only once +indz = np.where(coords['Z'] >= z.min()-0.1*coords['dz']) and \ + np.where(coords['Z'] <= z.max()+0.1*coords['dz']) +pressure_z, rho_z, Rgas_z = pressure_Z[indz], rho_Z[indz], Rgas_Z[indz] +# local proc 3D mhs arrays +pressure, rho = atm.mhs_3D.mhs_3D_profile(z, + pressure_z, + rho_z, + pressure_m, + rho_m + ) +magp = (Bx**2 + By**2 + Bz**2)/(2.*physical_constants['mu0']) +if rank ==0: + print'max B corona = ',magp[:,:,-1].max().decompose() + print'min B corona = ',magp[:,:,-1].min().decompose() +energy = atm.mhs_3D.get_internal_energy(pressure, + magp, + physical_constants) +Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) +Rgas[:] = Rgas_z +temperature = pressure/rho/Rgas +if not option_pars['l_hdonly']: + inan = np.where(magp <=1e-7*pressure.min()) + magpbeta = magp + magpbeta[inan] = 1e-7*pressure.min() # low pressure floor to avoid NaN + pbeta = pressure/magpbeta +else: + pbeta = magp+1.0 #dummy to avoid NaN +alfven = np.sqrt(2.*magp/rho) +cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) + +#print np.isnan(pbeta).any() + +#============================================================================ +# Save data for SAC and plotting +#============================================================================ +# set up data directory and file names +# may be worthwhile locating on /data if files are large +#datadir = os.path.expanduser('~/Documents/mhs_atmosphere/'+model_pars['model']+'/') +#datadir = os.path.expanduser('/fastdata-sharc/sm1bjs2/mhs_atmosphere/'+model_pars['model']+'/') +datadir = os.path.expanduser('./mhs_atmosphere/'+model_pars['model']+'/') +filename = datadir + model_pars['model'] + option_pars['suffix'] +if not os.path.exists(datadir): + os.makedirs(datadir) +sourcefile = datadir + model_pars['model'] + '_sources' + option_pars['suffix'] +aux3D = datadir + model_pars['model'] + '_3Daux' + option_pars['suffix'] +aux1D = datadir + model_pars['model'] + '_1Daux' + option_pars['suffix'] +# save the variables for the initialisation of a SAC simulation +atm.mhs_snapshot.save_SACvariables( + filename, + rho, + Bx, + By, + Bz, + energy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) +# save the balancing forces as the background source terms for SAC simulation +atm.mhs_snapshot.save_SACsources( + sourcefile, + Fx, + Fy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) +# save auxilliary variable and 1D profiles for plotting and analysis +atm.mhs_snapshot.save_auxilliary3D( + aux3D, + pressure_m, + rho_m, + temperature, + pbeta, + alfven, + cspeed, + Btensx, + Btensy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) +if rank==0: + atm.mhs_snapshot.save_auxilliary1D( + aux1D, + pressure_Z, + rho_Z, + Rgas_Z, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + rank=rank, + collective=False + ) +if rho.min()<0 or pressure.min()<0: + print"FAIL rank {}: negative rho.min() {} and/or pressure.min() {}.".format( + rank,rho.min(),pressure.min()) +FWHM = 2*np.sqrt(np.log(2))*model_pars['radial_scale'] +print'FWHM(0) =',FWHM + diff --git a/examples/mhs_atmosphere/plots/val_mtw_plot_3d.py b/examples/mhs_atmosphere/plots/val_mtw_plot_3d.py index 074d691..d574318 100644 --- a/examples/mhs_atmosphere/plots/val_mtw_plot_3d.py +++ b/examples/mhs_atmosphere/plots/val_mtw_plot_3d.py @@ -86,4 +86,4 @@ if nix[0].size > 0: seeds.T[1][nix] = - maxr * .99 atm.mhs_plot.make_3d_plot(ds, figname, fields=fkeys, figxy=figxy, view=view, - seeds=seeds) + seeds=seeds, offscreen=False) diff --git a/examples/mhs_atmosphere/spruit_atmosphere.py b/examples/mhs_atmosphere/spruit_atmosphere.py index c0da04d..bbf8133 100644 --- a/examples/mhs_atmosphere/spruit_atmosphere.py +++ b/examples/mhs_atmosphere/spruit_atmosphere.py @@ -40,9 +40,11 @@ rank = comm.Get_rank() size = comm.Get_size() + collective=True l_mpi = True l_mpi = l_mpi and (size != 1) except ImportError: + collective=False l_mpi = False rank = 0 size = 1 @@ -130,17 +132,21 @@ coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1], coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]] + axindex=np.arange(model_pars['Nxyz'][0]) # split the grid between processes for mpi if l_mpi: x_chunks = np.array_split(ax, size, axis=0) y_chunks = np.array_split(ay, size, axis=0) z_chunks = np.array_split(az, size, axis=0) + i_chunks = np.array_split(axindex, size, axis=0) x = comm.scatter(x_chunks, root=0) y = comm.scatter(y_chunks, root=0) z = comm.scatter(z_chunks, root=0) + xindex=i_chunks[rank] else: x, y, z = ax, ay, az + xindex=axindex x = u.Quantity(x, unit=coords['xmin'].unit) y = u.Quantity(y, unit=coords['ymin'].unit) @@ -221,43 +227,6 @@ energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants) - #============================================================================ - # Save data for SAC and plotting - #============================================================================ - # set up data directory and file names - # may be worthwhile locating on /data if files are large - datadir = os.path.expanduser('~/Documents/mhs_atmosphere/'+model_pars['model']+'/') - filename = datadir + model_pars['model'] + option_pars['suffix'] - if not os.path.exists(datadir): - os.makedirs(datadir) - sourcefile = datadir + model_pars['model'] + '_sources' + option_pars['suffix'] - aux3D = datadir + model_pars['model'] + '_3Daux' + option_pars['suffix'] - aux1D = datadir + model_pars['model'] + '_1Daux' + option_pars['suffix'] - print 'units of Bx =',Bx.unit - # save the variables for the initialisation of a SAC simulation - atm.mhs_snapshot.save_SACvariables( - filename, - rho, - Bx, - By, - Bz, - energy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) - # save the balancing forces as the background source terms for SAC simulation - atm.mhs_snapshot.save_SACsources( - sourcefile, - Fx, - Fy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) - # save auxilliary variable and 1D profiles for plotting and analysis Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) Rgas[:] = Rgas_z temperature = pressure/rho/Rgas @@ -274,28 +243,81 @@ alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2, 0].decompose(),\ alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2,-1].decompose() cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) +#============================================================================ +# Save data for SAC and plotting +#============================================================================ +# set up data directory and file names +# may be worthwhile locating on /data if files are large + datadir = os.path.expanduser('~/Documents/mhs_atmosphere/'+model_pars['model']+'/') + filename = datadir + model_pars['model'] + option_pars['suffix'] + if not os.path.exists(datadir): + os.makedirs(datadir) + sourcefile = datadir + model_pars['model'] + '_sources' + option_pars['suffix'] + aux3D = datadir + model_pars['model'] + '_3Daux' + option_pars['suffix'] + aux1D = datadir + model_pars['model'] + '_1Daux' + option_pars['suffix'] + # save the variables for the initialisation of a SAC simulation + atm.mhs_snapshot.save_SACvariables( + filename, + rho, + Bx, + By, + Bz, + energy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) + # save the balancing forces as the background source terms for SAC simulation + atm.mhs_snapshot.save_SACsources( + sourcefile, + Fx, + Fy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) + # save auxilliary variable and 1D profiles for plotting and analysis atm.mhs_snapshot.save_auxilliary3D( - aux3D, - pressure_m, - rho_m, - temperature, - pbeta, - alfven, - cspeed, - Btensx, - Btensy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) - atm.mhs_snapshot.save_auxilliary1D( - aux1D, - pressure_Z, - rho_Z, - Rgas_Z, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + aux3D, + pressure_m, + rho_m, + temperature, + pbeta, + alfven, + cspeed, + Btensx, + Btensy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) + if rank==0: + atm.mhs_snapshot.save_auxilliary1D( + aux1D, + pressure_Z, + rho_Z, + Rgas_Z, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + rank=rank, + collective=False + ) + if rho.min()<0 or pressure.min()<0: + print"FAIL rank {}: negative rho.min() {} and/or pressure.min() {}."\ + .format(rank,rho.min(),pressure.min()) + FWHM = 2*np.sqrt(np.log(2))*model_pars['radial_scale'] + print'FWHM(0) =',FWHM diff --git a/examples/mhs_atmosphere/val_mtw_atmosphere.py b/examples/mhs_atmosphere/val_mtw_atmosphere.py index 26c18ca..0c5ac85 100644 --- a/examples/mhs_atmosphere/val_mtw_atmosphere.py +++ b/examples/mhs_atmosphere/val_mtw_atmosphere.py @@ -29,7 +29,7 @@ import numpy as np import pysac.mhs_atmosphere as atm import astropy.units as u -from pysac.mhs_atmosphere.parameters.model_pars import paper2d as model_pars +from pysac.mhs_atmosphere.parameters.model_pars import mfe_setup as model_pars #============================================================================== #check whether mpi is required and the number of procs = size #============================================================================== @@ -40,9 +40,11 @@ rank = comm.Get_rank() size = comm.Get_size() + collective=True l_mpi = True l_mpi = l_mpi and (size != 1) except ImportError: + collective=False l_mpi = False rank = 0 size = 1 @@ -58,7 +60,6 @@ #obtain code coordinates and model parameters in astropy units coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz'])) - #interpolate the hs 1D profiles from empirical data source[s] empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu']) @@ -103,17 +104,21 @@ coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1], coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]] +axindex=np.arange(model_pars['Nxyz'][0]) # split the grid between processes for mpi if l_mpi: x_chunks = np.array_split(ax, size, axis=0) y_chunks = np.array_split(ay, size, axis=0) z_chunks = np.array_split(az, size, axis=0) + i_chunks = np.array_split(axindex, size, axis=0) x = comm.scatter(x_chunks, root=0) y = comm.scatter(y_chunks, root=0) z = comm.scatter(z_chunks, root=0) + xindex=i_chunks[rank] else: x, y, z = ax, ay, az + xindex=axindex x = u.Quantity(x, unit=coords['xmin'].unit) y = u.Quantity(y, unit=coords['ymin'].unit) @@ -195,6 +200,18 @@ energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants) +Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) +Rgas[:] = Rgas_z +temperature = pressure/rho/Rgas +if not option_pars['l_hdonly']: + inan = np.where(magp <=1e-7*pressure.min()) + magpbeta = magp + magpbeta[inan] = 1e-7*pressure.min() # low pressure floor to avoid NaN + pbeta = pressure/magpbeta +else: + pbeta = magp+1.0 #dummy to avoid NaN +alfven = np.sqrt(2.*magp/rho) +cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) #============================================================================ # Save data for SAC and plotting #============================================================================ @@ -209,71 +226,67 @@ aux1D = datadir + model_pars['model'] + '_1Daux' + option_pars['suffix'] # save the variables for the initialisation of a SAC simulation atm.mhs_snapshot.save_SACvariables( - filename, - rho, - Bx, - By, - Bz, - energy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + filename, + rho, + Bx, + By, + Bz, + energy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) # save the balancing forces as the background source terms for SAC simulation atm.mhs_snapshot.save_SACsources( - sourcefile, - Fx, - Fy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + sourcefile, + Fx, + Fy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) # save auxilliary variable and 1D profiles for plotting and analysis -Rgas = u.Quantity(np.zeros(x.shape), unit=Rgas_z.unit) -Rgas[:] = Rgas_z -temperature = pressure/rho/Rgas -if not option_pars['l_hdonly']: - inan = np.where(magp <=1e-7*pressure.min()) - magpbeta = magp - magpbeta[inan] = 1e-7*pressure.min() # low pressure floor to avoid NaN - pbeta = pressure/magpbeta -else: - pbeta = magp+1.0 #dummy to avoid NaN -alfven = np.sqrt(2.*magp/rho) -#if rank == 0: -# print'Alfven speed Z.min to Z.max =',\ -# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2, 0].decompose(),\ -# alfven[model_pars['Nxyz'][0]/2,model_pars['Nxyz'][1]/2,-1].decompose() -cspeed = np.sqrt(physical_constants['gamma']*pressure/rho) atm.mhs_snapshot.save_auxilliary3D( - aux3D, - pressure_m, - rho_m, - temperature, - pbeta, - alfven, - cspeed, - Btensx, - Btensy, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) -atm.mhs_snapshot.save_auxilliary1D( - aux1D, - pressure_Z, - rho_Z, - Rgas_Z, - option_pars, - physical_constants, - coords, - model_pars['Nxyz'] - ) + aux3D, + pressure_m, + rho_m, + temperature, + pbeta, + alfven, + cspeed, + Btensx, + Btensy, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + xindex, + rank=rank, + collective=collective + ) +if rank==0: + atm.mhs_snapshot.save_auxilliary1D( + aux1D, + pressure_Z, + rho_Z, + Rgas_Z, + option_pars, + physical_constants, + coords, + model_pars['Nxyz'], + rank=rank, + collective=False + ) if rho.min()<0 or pressure.min()<0: - print"FAIL: negative rho.min() {} and/or pressure.min() {}.".format( - rho.min(),pressure.min()) + print"FAIL rank {}: negative rho.min() {} and/or pressure.min() {}.".format( + rank,rho.min(),pressure.min()) FWHM = 2*np.sqrt(np.log(2))*model_pars['radial_scale'] print'FWHM(0) =',FWHM diff --git a/pysac/io/gdf_writer.py b/pysac/io/gdf_writer.py index c8be807..b021211 100644 --- a/pysac/io/gdf_writer.py +++ b/pysac/io/gdf_writer.py @@ -11,6 +11,20 @@ import h5py from h5py import h5s, h5p, h5fd +try: + from mpi4py import MPI + + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + size = comm.Get_size() + + l_mpi = True + l_mpi = l_mpi and (size != 1) +except ImportError: + l_mpi = False + rank = 0 + size = 1 + __all__ = ['write_field', 'write_field_raw', 'create_file', 'SimulationParameters'] @@ -68,7 +82,7 @@ def __repr__(self): def create_file(f, simulation_parameters, grid_dimensions, - data_author=None, data_comment=None): + data_author=None, data_comment=None, driver=None, comm=None): """ Do all the structral creation of a gdf file. @@ -100,8 +114,12 @@ def create_file(f, simulation_parameters, grid_dimensions, ----- GDF is defined here: https://bitbucket.org/yt_analysis/grid_data_format/ """ - if isinstance(f, basestring): - f = h5py.File(f, 'a') + if l_mpi: + if isinstance(f, basestring): + f = h5py.File(f, 'a', driver='mpio', comm=comm) + else: + if isinstance(f, basestring): + f = h5py.File(f, 'a') # "gridded_data_format" group g = f.create_group("gridded_data_format") @@ -266,10 +284,9 @@ def write_field_raw(gdf_file, field, field_shape=None, arr_slice=np.s_[:], def _write_dset_high(dset, data, arr_slice,collective=False): if collective: - with dset.collective: - dset[arr_slice] = np.ascontiguousarray(data) + dset[arr_slice,...] = data else: - dset[arr_slice] = np.ascontiguousarray(data) + dset[arr_slice,...] = np.ascontiguousarray(data) def _write_dset_low(dset, data, arr_slice, collective=False): memory_space = h5s.create_simple(data.shape) diff --git a/pysac/mhs_atmosphere/hs_model/VALIIIC.dat b/pysac/mhs_atmosphere/hs_model/VALIIIC.dat old mode 100755 new mode 100644 diff --git a/pysac/mhs_atmosphere/hs_model/hs_atmosphere.py b/pysac/mhs_atmosphere/hs_model/hs_atmosphere.py index 9658f53..4d57d61 100644 --- a/pysac/mhs_atmosphere/hs_model/hs_atmosphere.py +++ b/pysac/mhs_atmosphere/hs_model/hs_atmosphere.py @@ -246,6 +246,7 @@ def vertical_profile(Z, # - 84.*(g0*rdata[-i-5])*dz # + 6.*(g0*rdata[-i-6])*dz # )/162. + magp[-i-1] - magp[-i-0] + for i in range(1,Z.size): linp[-i-1] = (1152.*linp[-i] + 35.*(g0*rdata[-i-7])*dz @@ -254,5 +255,121 @@ def vertical_profile(Z, -784.*(g0*rdata[-i-4])*dz + 77.*(g0*rdata[-i-3])*dz )/1152. + magp[-i-1] - magp[-i-0] + thermalp_Z = linp + + +# try something with energy + eng_t = u.Quantity(np.ones(table['Z'].size), u.one) + eng_t *= Rgas*table_T/(physical_constants['gamma'] -1.0) + eng_Z = eng_t[4:-4].copy() + rho_ref = 200.0*rdata[1] #20.0*rdata[1] + rho_t = u.Quantity(np.ones(table['Z'].size), u.one) + rho_t *= rho_ref + rho_Z = rho_t[4:-4].copy() + table_TZ= table_T[4:-4].copy() +# magp0_Z=magp0[4:-4].copy + +# print magp0 + + for i in range(1,Z.size-1): + dg = 2.0*(physical_constants['gamma']-1.0)*eng_Z[i] + g0*dz.to('m') + dg = dg/(2.0*(physical_constants['gamma']-1.0)*eng_Z[i+1] -g0*dz ) + rho_Z[i+1] = rho_Z[i]*dg +# print rho_Z[i+1], Z[i] + + thermalp_Z = Rgas_Z*rho_Z*table_TZ + magp0 + rdata_Z = rho_Z + +# print thermalp_Z +# print (rdata_Z*g0-np.gradient(thermalp_Z)/dz.to('m'))/rdata_Z + return thermalp_Z, rdata_Z, Rgas_Z + +#============================================================================ +# Construct 3D hydrostatic test profile SNOW +#============================================================================ +def vertical_profile_test(Z, + table, + magp0, + physical_constants, dz, rhom_mean + ): + """Return the vertical profiles for thermal pressure and density in 1D. + Integrate in reverse from the corona to the photosphere to remove + sensitivity to larger chromospheric gradients.""" + g0 = physical_constants['gravity'].to('m s-2') + Rgas = u.Quantity(np.ones(table['Z'].size), u.one) + Rgas *= (physical_constants['boltzmann']/\ + physical_constants['proton_mass']/table['mu']).to('m2 K-1 s-2') + Rgas_Z = Rgas[4:-4].copy() + rdata = u.Quantity(table['rho'], copy=True).to('kg m-3') + rdata_Z = rdata[4:-4].copy() + magp = magp0.to('kg m-1 s-2') + # inverted SAC 4th order derivative scheme to minimise numerical error + """evaluate upper boundary pressure from equation of state + enhancement, + magp, which will be replaced by the mean magnetic pressure in the + corona, then integrate from inner next pressure + """ + table_T = u.Quantity(table['T']) + linp_1 = table_T[-1]*rdata[-1]*Rgas[-1] + magp[-1] + linp = u.Quantity(np.ones(len(Z)), unit=linp_1.unit) + linp[-1] = table_T[-5]*rdata[-5]*Rgas[-5] + magp[-1] + +# for i in range(1,Z.size): +# linp[-i-1] = (144.*linp[-i]+18.*linp[-i+1] +# -102.*(g0*rdata[-i-4] )*dz +# - 84.*(g0*rdata[-i-5])*dz +# + 6.*(g0*rdata[-i-6])*dz +# )/162. + magp[-i-1] - magp[-i-0] + + for i in range(1,Z.size): + linp[-i-1] = (1152.*linp[-i] + + 35.*(g0*rdata[-i-7])*dz + -112.*(g0*rdata[-i-6])*dz + -384.*(g0*rdata[-i-5])*dz + -784.*(g0*rdata[-i-4])*dz + + 77.*(g0*rdata[-i-3])*dz + )/1152. + magp[-i-1] - magp[-i-0] + + thermalp_Z = linp + + +# try something with energy + eng_t = u.Quantity(np.ones(table['Z'].size), u.one) + eng_t *= Rgas*table_T/(physical_constants['gamma'] -1.0) + eng_Z = eng_t[4:-4].copy() + rho_ref = np.max(rdata) #200.0*rdata[1] + rho_t = u.Quantity(np.ones(table['Z'].size), u.one) + rho_t *= rho_ref #20* +# rho_t -= rhom_min +# print rho_ref - rhom_min + rho_Z = rho_t[4:-4].copy() + rho_Z += rhom_mean + table_TZ= table_T[4:-4].copy() +# magp0_Z=magp0[4:-4].copy + +# print magp0 + + for i in range(1,Z.size-1): + dg = 2.0*(physical_constants['gamma']-1.0)*eng_Z[i] + g0*dz.to('m') + dg = dg/(2.0*(physical_constants['gamma']-1.0)*eng_Z[i+1] -g0*dz ) + rho_Z[i+1] = rho_Z[i]*dg +# print rho_Z[i+1], Z[i] + +# calculate thermal pressure + thermalp_Z = Rgas_Z*rho_Z*table_TZ + +# Correct to value at base + tcorrtemp = u.Quantity(np.ones(table['Z'].size), u.one) + tcorr = tcorrtemp[4:-4].copy() + tcorr *= table_TZ[1]*Rgas_Z[1]*rho_Z[1] - thermalp_Z[1] + thermalp_Z += tcorr + + print 'temperature at z=1 ', table_TZ[1], thermalp_Z[1]/(Rgas_Z[1]*rho_Z[1]) + +# Add the magnetic reference pressure + thermalp_Z += magp0 + rdata_Z = rho_Z + +# print thermalp_Z +# print (rdata_Z*g0-np.gradient(thermalp_Z)/dz.to('m'))/rdata_Z return thermalp_Z, rdata_Z, Rgas_Z diff --git a/pysac/mhs_atmosphere/hs_model/mcwhirter.dat b/pysac/mhs_atmosphere/hs_model/mcwhirter.dat old mode 100755 new mode 100644 diff --git a/pysac/mhs_atmosphere/mhs_model/__init__.py b/pysac/mhs_atmosphere/mhs_model/__init__.py index c030f3d..0963f4f 100644 --- a/pysac/mhs_atmosphere/mhs_model/__init__.py +++ b/pysac/mhs_atmosphere/mhs_model/__init__.py @@ -2,3 +2,4 @@ #from mhs_3D import * import flux_tubes import mhs_3D +import obs_test diff --git a/pysac/mhs_atmosphere/mhs_model/flux_tubes.py b/pysac/mhs_atmosphere/mhs_model/flux_tubes.py index cc6884f..0748b1a 100644 --- a/pysac/mhs_atmosphere/mhs_model/flux_tubes.py +++ b/pysac/mhs_atmosphere/mhs_model/flux_tubes.py @@ -39,11 +39,14 @@ def get_flux_tubes( ) # parameters for matching Mumford,Fedun,Erdelyi 2014 + if option_pars['l_sunspot']: + Si = [[0.5]]*u.T # 128.5mT SI units + # parameters for matching Mumford,Fedun,Erdelyi 2014 if option_pars['l_mfe']: - Si = [[0.1436]]*u.T # 128.5mT SI units + Si = [[0.15]]*u.T # [[0.1436]]*u.T 128.5mT SI units # parameters for matching Gent,Fedun,Mumford,Erdelyi 2014 elif option_pars['l_single']: - Si = [[0.1]]*u.T # 100mT SI units + Si = [[-0.1]]*u.T # 100mT SI units # parameters for matching Gent,Fedun,Erdelyi 2014 flux tube pair elif option_pars['l_tube_pair']: xi, yi, Si = ( @@ -66,6 +69,54 @@ def get_flux_tubes( [ 50e-3] ], unit=u.T) )# 50mT SI + elif option_pars['l_tube_pair_test']: + xi, yi, Si = ( + u.Quantity([ + [ 1.0], + [ 1.0], + [- 1.0], + [- 1.0] + ], unit=u.Mm), + u.Quantity([ + [ 0.00], + [ 0.00], + [ 0.00], + [ 0.00] + ], unit=u.Mm), + u.Quantity([ + [ 50e-3], + [ 50e-3], + [ 50e-3], + [ 50e-3] + ], unit=u.T) + )# 50mT SI + elif option_pars['l_tube_pair_test_2']: + xi, yi, Si = ( + u.Quantity([ + [ 1.0], + [- 1.0] + ], unit=u.Mm), + u.Quantity([ + [ 0.00], + [ 0.00] + ], unit=u.Mm), + u.Quantity([ + [ 100e-3], + [-100e-3] + ], unit=u.T) + )# 50mT SI + elif option_pars['l_singletube']: + xi, yi, Si = ( + u.Quantity([ + [ 0.0] + ], unit=u.Mm), + u.Quantity([ + [ 0.00] + ], unit=u.Mm), + u.Quantity([ + [ 50e-3] + ], unit=u.T) + )# 50mT SI # parameters for matching Gent,Fedun,Erdelyi 2014 twisted flux tubes elif option_pars['l_multi_twist']: xi, yi, Si = ( @@ -160,7 +211,6 @@ def get_hmi_flux_tubes( "2014/07/05 23:59:55"), vso.attrs.Instrument('HMI'), vso.attrs.Physobs('LOS_magnetic_field')) - #print results.show() if l_newdata: if not os.path.exits(sunpydir): @@ -266,6 +316,11 @@ def construct_magnetic_field( B10dz=- 2 * z *B1z**2/z1**2 - B2z**2/z2 - B3z/z3 B20dz= 8*z**2*B1z**3/z1**4 - 2* B1z**2/z1**2 +2*B2z**3/z2**2 +2*B3z/z3**2 B30dz=-48*z**3*B1z**4/z1**6 +24*z*B1z**3/z1**4 -6*B2z**4/z2**3 -6*B3z/z3**3 + elif option_pars['l_B0_test']: + B0z = Bf1 * np.exp(-z/z1) + B10dz= - B0z/z1 + B20dz= + B0z/z1**2 + B30dz= - B0z/z1**3 else: raise ValueError("in mhs_model.flux_tubes.construct_magnetic_field \ option_pars all False for axial strength Z dependence") @@ -309,7 +364,7 @@ def construct_magnetic_field( B2x = (Bx * dxBx + By * dyBx + Bz * dzBx)/mu0 B2y = (Bx * dxBy + By * dyBy + Bz * dzBy)/mu0 - warnings.warn("pbbal.max() = {}".format(pbbal.max().decompose()), Warning) +# warnings.warn("pbbal.max() = {}".format(pbbal.max().decompose()), Warning) return pbbal, rho_1, Bx, By, Bz, B2x, B2y #============================================================================ @@ -346,6 +401,11 @@ def construct_pairwise_field(x, y, z, B10dz= -2*z*B1z/z1**2 - B2z/z2 - B3z/z3 B20dz= -2* B1z/z1**2 + 4*z**2*B1z/z1**4 + B2z/z2**2 + B3z/z3**2 B30dz= 12*z*B1z/z1**4 - 8*z**3*B1z/z1**6 - B2z/z2**3 - B3z/z3**3 + elif option_pars['l_B0_test']: + B0z = Bf1 * np.exp(-z/z1) + B10dz= - B0z/z1 + B20dz= + B0z/z1**2 + B30dz= - B0z/z1**3 else: #if option_pars['l_BO_quadz']: B1z = Bf1 * z1**2 / (z**2 + z1**2) @@ -410,6 +470,8 @@ def construct_pairwise_field(x, y, z, (x-xi) * fxyzi + (x-xj) * fxyzj ) Fy = - 2*Si*Sj/mu0 * G0ij*BB10dz2/f02 * ( (y-yi) * fxyzi + (y-yj) * fxyzj ) +# Fx = (x-xi)/mu0 * Si*Sj*G0ij*B0z2*(-2.*fxyzi*fxyzi/f02)*BB10dz2 +(x-xi)/mu0 *Si*Sj*G0ij*B0z3*BB20dz +# Fx = (y-yi)/mu0 * Si*Sj*G0ij*B0z2*(-2.*fxyzi*fxyzi/f02)*BB10dz2 +(y-yi)/mu0 *Si*Sj*G0ij*B0z2*B0z*BB10dz2 #Define derivatives of Bx dxiBx = - Si * (BB10dz * G0i) \ + 2 * Si * (x-xi)**2 * B10dz * B0z3 * G0i/f02 @@ -433,6 +495,6 @@ def construct_pairwise_field(x, y, z, B2y = (Bxi * dxjBy + Byi * dyjBy + Bzi * dzjBy + Bxj * dxiBy + Byj * dyiBy + Bzj * dziBy)/mu0 - print"pbbal.max() = ",pbbal.max() +# print"pbbal.max() = ",pbbal.max() return pbbal, rho_1, Fx, Fy, B2x, B2y diff --git a/pysac/mhs_atmosphere/mhs_model/obs_test.py b/pysac/mhs_atmosphere/mhs_model/obs_test.py new file mode 100644 index 0000000..830d0ee --- /dev/null +++ b/pysac/mhs_atmosphere/mhs_model/obs_test.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- + +#Try to read the fits file + +import numpy as np +from numpy import unravel_index + +#import astropy.units as u +from astropy.io import fits + +import matplotlib.cm as cm +import matplotlib.pyplot as plt + +import operator + +import astropy.units as u + +def read_fits_bz(fits_name): +#Read fits file + + opened_fits=fits.open(fits_name) + + data= opened_fits[0].data + +# print np.size(data) +# fig, ax = plt.subplots() +# CS = ax.contour(data) +# plt.show() + + return data[345:475,370:500] #data[370:450,370:500] + +def flat_fits(data,nt): +#flatten fits data so only strong flux tubes are present + +# Set a threshold for B_z +# flat_dat=np.where((np.abs(data) < bmax), 0.0,data) +# flat_dat=np.where((data > -1250.0), 0.0,data) + flat_dat=np.where((np.abs(data) < 100.0), 0.0,data) + + nrows=len(flat_dat) + ncols=len(flat_dat[0]) + + print nrows,ncols + + flat_dat_grad=np.gradient(flat_dat) + flat_dat_grad_i=flat_dat_grad[0] + flat_dat_grad_j=flat_dat_grad[1] + + tube_loc_i=[] + tube_loc_j=[] + tube_loc_B=[] + + for ii in range(2,nrows-1): + for jj in range(2,ncols-1): + if (flat_dat_grad_i[ii-1,jj]*flat_dat_grad_i[ii+1,jj] < 0) and (flat_dat_grad_j[ii,jj-1]*flat_dat_grad_j[ii,jj+1] < 0) and (np.abs(flat_dat[ii,jj]) > 100): + + tube_loc_i.append([ii]) + tube_loc_j.append([jj]) + tube_loc_B.append([flat_dat[ii,jj]]) +# print ii, jj, flat_dat[ii,jj] + +# print tube_loc_i +# print tube_loc_j + + fti=[] + ftj=[] + ftB=[] + + for ntube in range(0,nt): + # Find the maximum value +# loc=unravel_index(np.abs(tube_loc_B).argmax(), tube_loc_B.shape) + +# print tube_loc_B + loc, value = max(enumerate(np.abs(tube_loc_B)), key=operator.itemgetter(1)) +# print loc + #print tube_loc_i[loc] + #print tube_loc_j[loc] + + fti.append(tube_loc_i[loc]) + ftj.append(tube_loc_j[loc]) + ftB.append(tube_loc_B[loc]) + + tube_loc_B[loc]=[0.0] + """# Remove said flux tube and find next maximum + # flat_dat[loc[0]-5:loc[0]+5,loc[1]-5:loc[1]+5]=0.0 + i=loc[0] + i2=loc[0] + j=loc[1] + j2=loc[1] + while flat_dat[i-1,loc[1]] > flat_dat[i,loc[1]]: + i=i-1 + while flat_dat[i2+1,loc[1]] > flat_dat[i2,loc[1]]: + i2=i2+1 + while flat_dat[loc[0],j-1] > flat_dat[loc[0],j]: + j=j-1 + while flat_dat[loc[0],j2+1] > flat_dat[loc[0],j2]: + j2=j2+1 + flat_dat[i:i2,j:j2]=0.0 + """ + +# loc=unravel_index(np.abs(flat_dat).argmax(), flat_dat.shape) + +# print loc +# print flat_dat[loc] + +# flat_dat[loc[0]-5:loc[0]+5,loc[1]-5:loc[1]+5]=0.0 + +# print np.size(flat_dat) +# fig, ax = plt.subplots() +# CS = ax.contour(flat_dat) +# cbar = fig.colorbar(CS) +# plt.show() + + +############################################################### +# This doesnt work yet.... +############################################################### +# sort out x,y locations scaled between 0 and 1 +# fti *= 0.00769230769 #(nrows-1.0) +# ftj *= 0.00769230769 #(ncols-1.0) + +# Scale x,y locations to Mm +# fti *= 10.0*u.Mm +# fti -= 5.0*u.Mm +# ftj *= 10.0*u.Mm +# ftj -= 5.0*u.Mm + +# Flux tubes units to Tesla +# ftB *= 0.0001*u.T + + return fti,ftj,ftB + diff --git a/pysac/mhs_atmosphere/parameters/model_pars.py b/pysac/mhs_atmosphere/parameters/model_pars.py index 0558d31..b7741a2 100644 --- a/pysac/mhs_atmosphere/parameters/model_pars.py +++ b/pysac/mhs_atmosphere/parameters/model_pars.py @@ -8,13 +8,17 @@ import numpy as np import os import astropy.units as u +import astropy +l_astropy_0=True +if astropy.__version__[0]=='1': + l_astropy_0=False hmi_model = {'photo_scale': 0.6*u.Mm, #scale height for photosphere 'chrom_scale': 0.1*u.Mm, #scale height for chromosphere 'corona_scale': 2.5e3*u.Mm, #scale height for the corona - 'coratio': 0.03*u.one, #footpoint portion scaling as corona + 'coratio': 0*u.one, #footpoint portion scaling as corona 'model': 'hmi_model', - 'phratio': 0.15*u.one, #footpoint portion scaling as photosphere + 'phratio': 0*u.one, #footpoint portion scaling as photosphere 'pixel': 0.36562475*u.Mm, #(HMI pixel) 'radial_scale': 0.10979002*u.Mm, #=> FWHM = half pixel 'nftubes': 1, @@ -22,24 +26,43 @@ 'pBplus': 1e-3*u.T} hmi_model['chratio'] = 1*u.one - hmi_model['coratio'] - hmi_model['phratio'] +sunspot = {'photo_scale': 0.60*u.Mm, + 'chrom_scale': 0.1*u.Mm, + 'corona_scale': 2.5e3*u.Mm, #scale height for the corona + 'coratio': 0.03*u.one, + 'model': 'sunspot', + 'phratio': 0.0*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 3.6096*u.Mm, + #'radial_scale': 0.14*u.Mm, + 'nftubes': 1, + #'B_corona': 4.85e-4*u.T, + 'B_corona': 5.5e-4*u.T, + 'pBplus': 12.0e-4*u.T} +sunspot['chratio'] = 1*u.one - sunspot['coratio'] - sunspot['phratio'] +#if 1D or 2D set unused dimensions to 0, and unrequired xyz limits to 1. +sunspot['Nxyz'] = [512,512,256] # 3D grid +#sunspot['Nxyz'] = [128,128,64] # 3D grid +sunspot['xyz'] = [-25*u.Mm,25*u.Mm,-25*u.Mm,25*u.Mm,2.5*u.Mm,22.5*u.Mm] #grid size + mfe_setup = {'photo_scale': 0.60*u.Mm, 'chrom_scale': 0.4*u.Mm, - 'corona_scale': 0.25*u.Mm, #scale height for the corona + 'corona_scale': 0.22*u.Mm, #scale height for the corona 'coratio': 0.0*u.one, 'model': 'mfe_setup', 'phratio': 0.0*u.one, 'pixel': 0.36562475*u.Mm, #(HMI pixel) - 'radial_scale': 0.03938*u.Mm, + 'radial_scale': 0.1*u.Mm, #0.03938*u.Mm, #'radial_scale': 0.14*u.Mm, 'nftubes': 1, #'B_corona': 4.85e-4*u.T, - 'B_corona': 5.5e-4*u.T, - 'pBplus': 12.0e-4*u.T} + 'B_corona': 5.5e-5*u.T, #5.5e-4*u.T, + 'pBplus': 12.0e-5*u.T} mfe_setup['chratio'] = 1*u.one - mfe_setup['coratio'] - mfe_setup['phratio'] #if 1D or 2D set unused dimensions to 0, and unrequired xyz limits to 1. -mfe_setup['Nxyz'] = [128,128,128] # 3D grid -mfe_setup['Nxyz'] = [129,129,128] # 3D grid -mfe_setup['xyz'] = [-1*u.Mm,1*u.Mm,-1*u.Mm,1*u.Mm,3.6641221e-2*u.Mm,1.5877863*u.Mm] #grid size +#mfe_setup['Nxyz'] = [128,128,128] # 3D grid +mfe_setup['Nxyz'] = [129,129,259] # 3D grid +mfe_setup['xyz'] = [-2.5*u.Mm,2.5*u.Mm,-2.5*u.Mm,2.5*u.Mm,3.6641221e-2*u.Mm,10.0*u.Mm] #grid size spruit = {'photo_scale': 1.5*u.Mm, 'chrom_scale': 0.5*u.Mm, @@ -76,7 +99,7 @@ #paper1['Nxyz'] = [127,128,128] # 3D grid #paper1['xyz'] = [-1*u.Mm,1*u.Mm,-1*u.Mm,1*u.Mm,3.5e-2*u.Mm,1.6*u.Mm] #grid size -paper2a = {'photo_scale': 0.6*u.Mm, +paper2ar = {'photo_scale': 0.6*u.Mm, 'chrom_scale': 0.1*u.Mm, 'corona_scale': 2.5e3*u.Mm, #scale height for the corona 'coratio': 0.03*u.one, @@ -87,14 +110,74 @@ 'nftubes': 4, 'B_corona': 8.4e-4*u.T, 'pBplus': 1e-3*u.T} +paper2ar['chratio'] = 1*u.one - paper2ar['coratio'] - paper2ar['phratio'] +paper2ar['Nxyz'] = [102,102,400]#[160,80,432] # 3D grid +paper2ar['xyz'] = [-1.59*u.Mm,1.59*u.Mm,-0.79*u.Mm,0.79*u.Mm,0.*u.Mm,8.62*u.Mm] #grid size + +tutu = {'photo_scale': 0.6*u.Mm, + 'chrom_scale': 0.26*u.Mm, #0.1*u.Mm, + 'corona_scale': 6.5e3*u.Mm, #2.5e3*u.Mm, #scale height for the corona + 'coratio': 0.08*u.one, #0.03*u.one, + 'model': 'tutu', + 'phratio': 0.2*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 0.10979002*u.Mm, + 'nftubes': 4, + 'B_corona': 4.4e-4*u.T, + 'pBplus': 1.5e-3*u.T} +tutu['chratio'] = 1*u.one - tutu['coratio'] - tutu['phratio'] +tutu['Nxyz'] = [102,102,200]#[160,80,432] # 3D grid +tutu['xyz'] = [-2.0*u.Mm,2.0*u.Mm,-1.4*u.Mm,1.4*u.Mm,0.*u.Mm,2.0*u.Mm] #grid size + +loop = {'photo_scale': 0.6*u.Mm, + 'chrom_scale': 0.26*u.Mm, #0.1*u.Mm, + 'corona_scale': 2.5e3*u.Mm, #2.5e3*u.Mm, #scale height for the corona + 'coratio': 0.08*u.one, #0.03*u.one, + 'model': 'loop', + 'phratio': 0.2*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 0.10979002*u.Mm, + 'nftubes': 2, + 'B_corona': 0.0*u.T, + 'pBplus': 1.5e-3*u.T} +loop['chratio'] = 1*u.one - loop['coratio'] - loop['phratio'] +loop['Nxyz'] = [101,101,200]#[160,80,432] # 3D grid +loop['xyz'] = [-2.0*u.Mm,2.0*u.Mm,-1.4*u.Mm,1.4*u.Mm,0.*u.Mm,2.0*u.Mm] #grid size + +singletube = {'photo_scale': 0.6*u.Mm, + 'chrom_scale': 0.1*u.Mm, #0.1*u.Mm, + 'corona_scale': 2.0e0*u.Mm, #2.5e3*u.Mm, #scale height for the corona + 'coratio': 0.06*u.one, #0.03*u.one, + 'model': 'singletube', + 'phratio': 0.3*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 0.3*u.Mm, + 'nftubes': 1, + 'B_corona': 1.0e-4*u.T, + 'pBplus': 1.5e-3*u.T} +singletube['chratio'] = 1*u.one - singletube['coratio'] - singletube['phratio'] +singletube['Nxyz'] = [101,101,200]#[160,80,432] # 3D grid +singletube['xyz'] = [-2.0*u.Mm,2.0*u.Mm,-2.0*u.Mm,2.0*u.Mm,0.*u.Mm,2.0*u.Mm] #grid size + +paper2a = {'photo_scale': 0.6*u.Mm, + 'chrom_scale': 0.26*u.Mm, #0.1*u.Mm, + 'corona_scale': 6.5e3*u.Mm, #2.5e3*u.Mm, #scale height for the corona + 'coratio': 0.08*u.one, #0.03*u.one, + 'model': 'paper2a', + 'phratio': 0.2*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 0.10979002*u.Mm, + 'nftubes': 4, + 'B_corona': 4.4e-4*u.T, + 'pBplus': 1.5e-3*u.T} paper2a['chratio'] = 1*u.one - paper2a['coratio'] - paper2a['phratio'] -paper2a['Nxyz'] = [160,80,432] # 3D grid -paper2a['xyz'] = [-1.59*u.Mm,1.59*u.Mm,-0.79*u.Mm,0.79*u.Mm,0.*u.Mm,8.62*u.Mm] #grid size +paper2a['Nxyz'] = [102,102,200]#[160,80,432] # 3D grid +paper2a['xyz'] = [-1.59*u.Mm,1.59*u.Mm,-0.79*u.Mm,0.79*u.Mm,0.*u.Mm,4.0*u.Mm] #[-1.59*u.Mm,1.59*u.Mm,-0.79*u.Mm,0.79*u.Mm,0.*u.Mm,8.62*u.Mm] #grid size paper2b = {'photo_scale': 0.6*u.Mm, 'chrom_scale': 0.1*u.Mm, 'corona_scale': 2.5e3*u.Mm, #scale height for the corona - 'coratio': 0.03*u.one, + 'coratio': 0.032*u.one, 'model': 'paper2b', 'phratio': 0.0*u.one, 'pixel': 0.36562475*u.Mm, #(HMI pixel) @@ -103,7 +186,7 @@ 'B_corona': 8.2e-4*u.T, 'pBplus': 1.0e-3*u.T} paper2b['chratio'] = 1*u.one - paper2b['coratio'] - paper2b['phratio'] -paper2b['Nxyz'] = [50,50,140] # 3D grid +paper2b['Nxyz'] = [50,50,200] # 3D grid paper2b['xyz'] = [-0.49*u.Mm,0.49*u.Mm,-0.49*u.Mm,0.49*u.Mm,0*u.Mm,2.78*u.Mm] #grid size paper2c = {'photo_scale': 0.6*u.Mm, @@ -133,9 +216,60 @@ 'B_corona': 5.95e-4*u.T, 'pBplus': 1.0e-3*u.T} paper2d['chratio'] = 1*u.one - paper2d['coratio'] - paper2d['phratio'] -paper2d['Nxyz'] = [224,224,140] # 3D grid +paper2d['Nxyz'] = [50,50,140] #[224,224,140] # 3D grid paper2d['xyz'] = [-2.23*u.Mm,2.23*u.Mm,-0.79*u.Mm,0.79*u.Mm,0*u.Mm,2.78*u.Mm] #grid size +hs_test = {'photo_scale': 0.60*u.Mm, + 'chrom_scale': 0.4*u.Mm, + 'corona_scale': 0.22*u.Mm, #scale height for the corona + 'coratio': 0.0*u.one, + 'model': 'hs_test', + 'phratio': 0.0*u.one, + 'pixel': 0.36562475*u.Mm, #(HMI pixel) + 'radial_scale': 0.1*u.Mm, #0.03938*u.Mm, + #'radial_scale': 0.14*u.Mm, + 'nftubes': 0, + #'B_corona': 4.85e-4*u.T, + 'B_corona':0.0*u.T, #5.5e-4*u.T, + 'pBplus': 0.0*u.T} +hs_test['chratio'] = 1*u.one - hs_test['coratio'] - hs_test['phratio'] +#if 1D or 2D set unused dimensions to 0, and unrequired xyz limits to 1. +#mfe_setup['Nxyz'] = [128,128,128] # 3D grid +hs_test['Nxyz'] = [32,32,259] # 3D grid +hs_test['xyz'] = [-2.5*u.Mm,2.5*u.Mm,-2.5*u.Mm,2.5*u.Mm,3.6641221e-2*u.Mm,2.0*u.Mm] #grid size + +mhs_test = {'photo_scale': 0.1*u.Mm, + 'chrom_scale': 0.8*u.Mm, + 'corona_scale': 2.0*u.Mm, #scale height for the corona + 'coratio': 0.0*u.one, + 'model': 'mhs_test', + 'phratio': 0.0*u.one, + 'pixel': 0.36562475*u.Mm, + 'radial_scale': 0.1*u.Mm, + 'nftubes': 1, + 'B_corona': 0.0*u.T, #5.5e-6 + 'pBplus': 12.0e-5*u.T} +mhs_test['chratio'] = 1*u.one - mhs_test['coratio'] - mhs_test['phratio'] +#if 1D or 2D set unused dimensions to 0, and unrequired xyz limits to 1. +mhs_test['Nxyz'] = [129,129,259] # 3D grid +mhs_test['xyz'] = [-2.5*u.Mm,2.5*u.Mm,-2.5*u.Mm,2.5*u.Mm,3.6641221e-2*u.Mm,2.0*u.Mm] #grid size + +obs_test = {'photo_scale': 0.25*u.Mm, # 0.3 + 'chrom_scale': 0.4*u.Mm, + 'corona_scale': 2.0*u.Mm, #scale height for the corona + 'coratio': 0.0*u.one, + 'model': 'obs_test', + 'phratio': 1.0*u.one, + 'pixel': 0.36562475*u.Mm, + 'radial_scale': 0.025*u.Mm, #0.2 + 'nftubes': 1, + 'B_corona': 0.0*u.T, #5.5e-6 + 'pBplus': 1.0e-4*u.T} +obs_test['chratio'] = 1*u.one - obs_test['coratio'] - obs_test['phratio'] +#if 1D or 2D set unused dimensions to 0, and unrequired xyz limits to 1. +obs_test['Nxyz'] = [32,32,32] # 3D grid +obs_test['xyz'] = [-6.0*u.Mm,6.0*u.Mm,-6.0*u.Mm,6.0*u.Mm,0.0*u.Mm,2.0*u.Mm] #grid size + def get_coords(Nxyz, xyz): """ @@ -161,14 +295,71 @@ def get_coords(Nxyz, xyz): return coords #----------------------------------------------------------------------------- + +def get_hmi_coords( + Nxyz,xyz, + indx, + dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram_fits', + sunpydir = os.path.expanduser(os.path.expanduser('~')+'/sunpy/data/'), + figsdir = os.path.expanduser(os.path.expanduser('~')+'/figs/hmi/'), + l_newdata = True, + interpfactor=5, + frac=[0,0], + rank=0, + lmpi=False + ): + """ + get_coords returns a non-dimensional dictionary describing the domain + coordinates. + """ + dz=(xyz[5]-xyz[4])/(Nxyz[2]-1) + Z = u.Quantity(np.linspace(xyz[4].value,xyz[5].value,Nxyz[2]), + unit=u.Mm) + Zext = u.Quantity(np.linspace(Z.min().value-4.*dz.value, + Z.max().value+4.*dz.value, Nxyz[2]+8), + unit=u.Mm) + s,x,y,FWHM,sdummy,xdummy,ydummy,c1,c2=get_hmi_map( + indx, + dataset = dataset, + sunpydir = sunpydir, + figsdir = figsdir, + l_newdata = l_newdata, + interpfactor=interpfactor, + frac=frac, + rank=rank, + lmpi=lmpi + ) + xmin=x.min()+(x.max()-x.min())*frac[0] + xmax=x.min()+(x.max()-x.min())*(1-frac[0]) + ymin=y.min()+(y.max()-y.min())*frac[1] + ymax=y.min()+(y.max()-y.min())*(1-frac[1]) + coords = { + 'dx':(xmax-xmin)/(Nxyz[0]-1), + 'dy':(ymax-ymin)/(Nxyz[1]-1), + 'dz':(xyz[5]-xyz[4])/(Nxyz[2]-1), + 'xmin':xmin.to(u.Mm), + 'xmax':xmax.to(u.Mm), + 'ymin':ymin.to(u.Mm), + 'ymax':ymax.to(u.Mm), + 'zmin':xyz[4], + 'zmax':xyz[5], + 'Z':Z, + 'Zext':Zext + } + + return coords +#----------------------------------------------------------------------------- # def get_hmi_map( - model_pars, option_pars, - indx, - dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram_fits', - sunpydir = os.path.expanduser('~/sunpy/data/'), - figsdir = os.path.expanduser('~/figs/hmi/'), - l_newdata = False + indx, + dataset = 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram.fits', + sunpydir = os.path.expanduser(os.path.expanduser('~')+'/sunpy/data/'), + figsdir = os.path.expanduser(os.path.expanduser('~')+'/figs/hmi/'), + l_newdata = False, + interpfactor=5, + frac=[0,0], + rank=0, + lmpi=False ): """ indx is 4 integers: lower and upper indices each of x,y coordinates # dataset of the form 'hmi_m_45s_2014_07_06_00_00_45_tai_magnetogram_fits' @@ -178,83 +369,145 @@ def get_hmi_map( import sunpy.map client = vso.VSOClient() results = client.query(vso.attrs.Time("2014/07/05 23:59:50", - "2014/07/05 23:59:55"), + "2014/07/05 23:59:55"), vso.attrs.Instrument('HMI'), vso.attrs.Physobs('LOS_magnetic_field')) #print results.show() + if lmpi: + from mpi4py import MPI + comm = MPI.COMM_WORLD if l_newdata: - if not os.path.exists(sunpydir): - raise ValueError("in get_hmi_map set 'sunpy' dir for vso data\n"+ - "for large files you may want link to local drive rather than network") - client.get(results).wait(progress=True) + if rank==0: + if not os.path.exists(sunpydir): + raise ValueError("in get_hmi_map set 'sunpy' dir for vso data\n"+ + "for large files you may want link to local drive rather than network" + ) + client.get(results).wait(progress=True) + lwait = False + else: + lwait = None + if lmpi: + lwait = comm.bcast(lwait, root=0) + if not os.path.exists(figsdir): os.makedirs(figsdir) + print sunpydir, dataset hmi_map = sunpy.map.Map(sunpydir+dataset) #hmi_map = hmi_map.rotate() #hmi_map.peek() s = hmi_map.data[indx[0]:indx[1],indx[2]:indx[3]] #units of Gauss Bz - s *= u.G nx = s.shape[0] ny = s.shape[1] - nx2, ny2 = 2*nx, 2*ny # size of interpolant + nx_int, ny_int = interpfactor*(nx-1)+1, interpfactor*(ny-1)+1 # size of interpolant #pixel size in arc seconds - dx, dy = hmi_map.scale.items()[0][1],hmi_map.scale.items()[1][1] - x, y = np.mgrid[ - hmi_map.xrange[0]+indx[0]*dx:hmi_map.xrange[0]+indx[1]*dx:1j*nx2, - hmi_map.xrange[0]+indx[2]*dy:hmi_map.xrange[0]+indx[3]*dy:1j*ny2 - ] - #arrays to interpolate s from/to - fx = u.Quantity(np.linspace(x.min().value,x.max().value,nx), unit=x.unit) - fy = u.Quantity(np.linspace(y.min().value,y.max().value,ny), unit=y.unit) - xnew = u.Quantity(np.linspace(x.min().value,x.max().value,nx2), unit=x.unit) - ynew = u.Quantity(np.linspace(y.min().value,y.max().value,ny2), unit=y.unit) - f = RectBivariateSpline(fx,fy,s.to(u.T)) - #The initial model assumes a relatively small region, so a linear + if not l_astropy_0: + dx, dy = hmi_map.scale.items()[0][1],hmi_map.scale.items()[1][1] + x_int, y_int = np.mgrid[ + hmi_map.xrange[0]+indx[0]*dx: + hmi_map.xrange[0]+(indx[1]-1)*dx:1j*nx_int, + hmi_map.yrange[0]+indx[2]*dy: + hmi_map.yrange[0]+(indx[3]-1)*dy:1j*ny_int + ] + x, y = np.mgrid[ + hmi_map.xrange[0]+indx[0]*dx: + hmi_map.xrange[0]+(indx[1]-1)*dx:1j*nx, + hmi_map.yrange[0]+indx[2]*dy: + hmi_map.yrange[0]+(indx[3]-1)*dy:1j*ny + ] + else: + dx, dy = hmi_map.scale.x.value,hmi_map.scale.y.value + x_int, y_int = np.mgrid[ + hmi_map.xrange[0].value+indx[0]*dx: + hmi_map.xrange[0].value+indx[1]*dx:1j*nx_int, + hmi_map.yrange[0].value+indx[2]*dy: + hmi_map.yrange[0].value+indx[3]*dy:1j*ny_int + ] + x, y = np.mgrid[ + hmi_map.xrange[0].value+indx[0]*dx: + hmi_map.xrange[0].value+indx[1]*dx:1j*nx, + hmi_map.yrange[0].value+indx[2]*dy: + hmi_map.yrange[0].value+indx[3]*dy:1j*ny + ] + #arrays to interpolate s from/to + fx = np.linspace(x_int.min(),x_int.max(),nx) + fy = np.linspace(y_int.min(),y_int.max(),ny) + xnew = np.linspace(x_int.min(),x_int.max(),nx_int) + ynew = np.linspace(y_int.min(),y_int.max(),ny_int) + f = RectBivariateSpline(fx,fy,s) + #The initial model assumes a relatively small region, so a linear #Cartesian map is applied here. Consideration may be required if larger #regions are of interest, where curvature or orientation near the lim #of the surface is significant. s_int = f(xnew,ynew) #interpolate s and convert units to Tesla - s_int /= 4. # rescale s as extra pixels will sum over FWHM - x_int = x * 7.25e5 * u.m #convert units to metres - y_int = y * 7.25e5 * u.m - dx_int = dx * 7.25e5 * u.m - dy_int = dy * 7.25e5 * u.m - FWHM = 0.5*(dx_SI+dy_SI) - smax = max(abs(s.min()),abs(s.max())) # set symmetric plot scale - cmin = -smax*1e-4 - cmax = smax*1e-4 -# -# filename = 'hmi_map' -# import loop_plots as mhs -# mhs.plot_hmi( -# s*1e-4,x_SI.min(),x_SI.max(),y_SI.min(),y_SI.max(), -# cmin,cmax,filename,savedir,annotate = '(a)' -# ) -# filename = 'hmi_2x2_map' -# mhs.plot_hmi( -# s_SI*4,x_SI.min(),x_SI.max(),y_SI.min(),y_SI.max(), -# cmin,cmax,filename,savedir,annotate = '(a)' -# ) -# -# return s_SI, x_SI, y_SI, nx2, ny2, dx_SI, dy_SI, cmin, cmax, FWHM - dz=(xyz[5]-xyz[4])/(Nxyz[2]-1) - Z = u.Quantity(np.linspace(xyz[4].value, xyz[5].value, Nxyz[2]), unit=xyz.unit) - Zext = u.Quantity(np.linspace(Z.min().value-4.*dz.value, Z.max().value+4.*dz.value, Nxyz[2]+8), unit=Z.unit) - coords = { - 'dx':(xyz[1]-xyz[0])/(Nxyz[0]-1), - 'dy':(xyz[3]-xyz[2])/(Nxyz[1]-1), - 'dz':(xyz[5]-xyz[4])/(Nxyz[2]-1), - 'xmin':xyz[0], - 'xmax':xyz[1], - 'ymin':xyz[2], - 'ymax':xyz[3], - 'zmin':xyz[4], - 'zmax':xyz[5], - 'Z':Z, - 'Zext':Zext - } + interp_scale = 1. #contribution to field strength from location + interp_scale += 2. #contribution to field strength from 4 lateral neighbours + interp_scale += 0.2706705664732254 #contribution from nearest 4 diagonal neighbours + interp_scale += 0.036631277777468357 #from 2nd nearest 4 lateral + interp_scale += 0.026951787996341868 #from 1st 8 offset diagonals + interp_scale += 0.00067092525580502371 #from 2nd 4 diagonals + xq = u.Quantity(x_int * 7.25e5, unit= u.m) + yq = u.Quantity(y_int * 7.25e5, unit= u.m) + xq_int = u.Quantity(x_int * 7.25e5, unit= u.m) + yq_int = u.Quantity(y_int * 7.25e5, unit= u.m) + sq = u.Quantity(s * 1e-4, unit= u.T) + sq_int = u.Quantity(s_int * 1e-4 / interp_scale, unit= u.T) - return coords + dx *= 7.25e5 * u.m + dy *= 7.25e5 * u.m + pixel = (dx+dy)*0.5 + FWHM = 2*pixel/interpfactor + + import matplotlib.pyplot as plt + xmin=xq.to(u.Mm).min().value+(xq.to(u.Mm).max().value-xq.to(u.Mm).min().value)*frac[0] + xmax=xq.to(u.Mm).min().value+(xq.to(u.Mm).max().value-xq.to(u.Mm).min().value)*(1-frac[0]) + ymin=yq.to(u.Mm).min().value+(yq.to(u.Mm).max().value-yq.to(u.Mm).min().value)*frac[0] + ymax=yq.to(u.Mm).min().value+(yq.to(u.Mm).max().value-yq.to(u.Mm).min().value)*(1-frac[0]) + + cmax = max(-s_int.min()*1e-4, s_int.max()*1e-4) + cmin=-cmax + + plt.figure() + plt.pcolormesh(xq_int.T.to(u.Mm).value, yq_int.T.to(u.Mm).value, s_int.T*1e-4, + vmin=cmin, vmax=cmax) + plt.axis([xq.to(u.Mm).value.min(),xq.to(u.Mm).value.max(),yq.to(u.Mm).value.min(),yq.to(u.Mm).value.max()]) + plt.xlabel('lon [Mm]') + plt.ylabel('lat [Mm]') + cbar = plt.colorbar() + cbar.ax.set_ylabel(r'$B_z$ [T]') + cbar.solids.set_edgecolor("face") + plt.savefig('interpolated_hmi.png') + plt.close() + + plt.figure() + plt.pcolormesh(xq.T.to(u.Mm).value, yq.T.to(u.Mm).value, sq.T.value, + vmin=cmin, vmax=cmax) + plt.axis([xq.to(u.Mm).value.min(),xq.to(u.Mm).value.max(),yq.to(u.Mm).value.min(),yq.to(u.Mm).value.max()]) + plt.xlabel('lon [Mm]') + plt.ylabel('lat [Mm]') + cbar = plt.colorbar() + cbar.ax.set_ylabel(r'$B_z$ [T]') + cbar.solids.set_edgecolor("face") + plt.savefig('hmi.png') + plt.close() + + cmax = max(-sq[s.shape[0]*frac[0]:s.shape[0]*(1-frac[0]), + s.shape[1]*frac[1]:s.shape[1]*(1-frac[1])].value.min(), + sq[s.shape[0]*frac[0]:s.shape[0]*(1-frac[0]), + s.shape[1]*frac[1]:s.shape[1]*(1-frac[1])].value.max())*1.1 + cmin=-cmax + + plt.figure() + plt.pcolormesh(xq_int.T.to(u.Mm).value, yq_int.T.to(u.Mm).value, s_int.T*1e-4, + vmin=cmin, vmax=cmax) + plt.xlabel('lon [Mm]') + plt.ylabel('lat [Mm]') + plt.axis([xmin,xmax,ymin,ymax]) + cbar = plt.colorbar() + cbar.ax.set_ylabel(r'$B_z$ [T]') + cbar.solids.set_edgecolor("face") + plt.savefig('model_domain_hmi.png') + plt.close() + return sq_int, xq_int, yq_int, FWHM, sq, xq, yq, cmax, cmin diff --git a/pysac/mhs_atmosphere/parameters/options.py b/pysac/mhs_atmosphere/parameters/options.py index 6c0c0cc..d9f7a12 100644 --- a/pysac/mhs_atmosphere/parameters/options.py +++ b/pysac/mhs_atmosphere/parameters/options.py @@ -28,65 +28,96 @@ def set_options(model, l_mpi, l_gdf=True): 'l_B0_expz': False,# Z-depend of Bz(r=0) exponentials 'l_B0_quadz': False,# Z-depend of Bz(r=0) polynomials + exponential 'l_B0_rootz': False,# Z-depend of Bz(r=0) sqrt polynomials + 'l_B0_test': False,# Z-depend of Bz(r=0) single exponential 'l_single': False,# only one flux tube 'l_hmi': False,# construct photopheric map of Bz from HMI/SDI 'l_tube_pair': False,# pair of flux tubes + 'l_tube_pair_test': False,# pair of identical flux tubes + 'l_tube_pair_test_2': False,# pair of opposite flux tubes + 'l_singletube': False,# single flux tube 'l_multi_netwk': False,# multiple flux tubes as described in GFE (2014) 'l_multi_lanes': False,# multiple flux tubes as described in GFE (2014) 'l_multi_twist': False,# multiple flux tubes as described in GFE (2014) 'l_2D_loop': False,# make a 2D loop with sinusoidal Bz(x,0,0) 'l_mfe': False,# model Viktor's model from MFE (2014) + 'l_sunspot': False,# model sunspot (2016) 'l_atmos_val3c_mtw':False,# interpolate composite VAL3c+MTW atmosphere 'suffix': '.gdf' } #revise optional parameters depending on configuration required if model['model'] == 'hmi_model': - option_pars['l_hmi'] = True - option_pars['l_B0_expz'] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_hmi'] = True + option_pars['l_B0_expz'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'sunspot': + option_pars['l_single'] = True + option_pars['l_sunspot'] = True + option_pars['l_B0_expz'] = True + option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'mfe_setup': - option_pars['l_single'] = True - option_pars['l_mfe'] = True - option_pars['l_B0_expz'] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_single'] = True + option_pars['l_mfe'] = True + option_pars['l_B0_expz'] = True + option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'spruit': - option_pars['l_single'] = True - option_pars['l_spruit'] = True + option_pars['l_single'] = True + option_pars['l_spruit'] = True if model['model'] == 'paper1': - option_pars['l_ambB'] = True + option_pars['l_ambB'] = True option_pars['l_B0_expz'] = True option_pars['l_single'] = True option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'paper2a': - option_pars['l_ambB'] = True - option_pars['l_B0_expz'] = True - option_pars['l_tube_pair'] = True + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_tube_pair'] = True option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'paper2b': - option_pars['l_ambB'] = True - option_pars['l_B0_expz'] = True - option_pars['l_multi_twist' ] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_multi_twist' ] = True + option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'paper2c': - option_pars['l_ambB'] = True - option_pars['l_B0_expz'] = True - option_pars['l_multi_netwk'] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_multi_netwk'] = True + option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'paper2d': - option_pars['l_ambB'] = True - option_pars['l_B0_expz'] = True - option_pars['l_multi_lanes' ] = True - option_pars['l_atmos_val3c_mtw'] = True - if model['model'] == 'hmi_model': - option_pars['l_B0_quadz'] = True - option_pars['l_single'] = True - option_pars['l_hmi'] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_multi_lanes' ] = True + option_pars['l_atmos_val3c_mtw'] = True if model['model'] == 'loop_model': - option_pars['l_B0_quadz'] = True - option_pars['l_single'] = True - option_pars['l_2D_loop'] = True - option_pars['l_atmos_val3c_mtw'] = True + option_pars['l_B0_quadz'] = True + option_pars['l_single'] = True + option_pars['l_2D_loop'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'tutu': + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_tube_pair_test'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'loop': + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_tube_pair_test_2'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'singletube': + option_pars['l_ambB'] = True + option_pars['l_B0_expz'] = True + option_pars['l_single'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'mhs_test': + option_pars['l_single'] = True +# option_pars['l_mfe'] = True + option_pars['l_B0_expz'] = True + option_pars['l_atmos_val3c_mtw'] = True + if model['model'] == 'obs_test': +# option_pars['l_single'] = True +# option_pars['l_mfe'] = True +# option_pars['l_B0_expz'] = True + option_pars['l_B0_test'] = True + option_pars['l_atmos_val3c_mtw'] = True if l_mpi: option_pars['l_mpi'] = True else: diff --git a/pysac/mhs_atmosphere/plot/mhs_plot.py b/pysac/mhs_atmosphere/plot/mhs_plot.py index d6f166d..b00ab38 100755 --- a/pysac/mhs_atmosphere/plot/mhs_plot.py +++ b/pysac/mhs_atmosphere/plot/mhs_plot.py @@ -415,12 +415,13 @@ def make_2d_plot(ds, var_field, figname, normal = ['y',64], ## Fieldline Generation ##============================================================================ def make_3d_plot(ds, figname, - fields= ['thermal_pressure','plasma_beta', + fields = ['thermal_pressure','plasma_beta', 'mag_field_x','mag_field_y','mag_field_z'], - figxy=[500,550], - view=(-45., 90., 20., np.array([0,0,3.75])), - seeds=np.array([[0,0,1]]) - ): + figxy = [500,550], + view = (-45., 90., 20., np.array([0,0,3.75])), + seeds = np.array([[0,0,1]]), + offscreen = False + ): """Make a 3D rendition of the atmosphere including volume filling, iso surfaces and field lines. ds: gdf data set @@ -441,6 +442,9 @@ def make_3d_plot(ds, figname, # vector_field = 'vorticity' # mlab.options.offscreen = True + if offscreen: + #mlab.engine.current_scene.scene.off_screen_rendering = True + mlab.options.offscreen = True scene = mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0.5, 0.5, 0.5),size=figxy) diff --git a/pysac/mhs_atmosphere/save_atmos/mhs_snapshot.py b/pysac/mhs_atmosphere/save_atmos/mhs_snapshot.py index 315d72a..1604a02 100644 --- a/pysac/mhs_atmosphere/save_atmos/mhs_snapshot.py +++ b/pysac/mhs_atmosphere/save_atmos/mhs_snapshot.py @@ -8,7 +8,6 @@ import pysac.io.gdf_writer as gdf import h5py import astropy.units as u -#import h5py as h5 ##============================================================================ ## Save a file!!! @@ -21,135 +20,169 @@ and add symlink to ${HOME} to avoid exceeding quota. """ def save_SACvariables( - filename, - rho, - Bx, - By, - Bz, - energy, - option_pars, - physical_constants, - coords, - Nxyz - ): + filename, + rho, + Bx, + By, + Bz, + energy, + option_pars, + physical_constants, + coords, + Nxyz, + xindex, + rank=0, + collective=False, + ): """ Save the background variables for a SAC model in hdf5 (gdf default) format after collating the data from mpi sub processes if necessary. """ - rank = 0 - if option_pars['l_mpi']: - from mpi4py import MPI - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - gather_vars = [ - rho.to(u.Unit('kg/m**3')), - Bx.to(u.Unit('T')), - By.to(u.Unit('T')), - Bz.to(u.Unit('T')), - energy.to(u.Unit('kg/(m s**2)')) - ] - concat_vars = [] - for var in gather_vars: - concat_vars.append(comm.gather(var, root=0)) - - if rank == 0: - out_vars = [] - for cvar in concat_vars: - out_vars.append(np.concatenate(cvar, axis=0)) - - rho,Bx,By,Bz,energy = out_vars if rank == 0: print'writing',filename print'SAC background atmosphere' - - grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] - left_edge = u.Quantity([coords['xmin'], - coords['ymin'], - coords['zmin']]).to(u.m) - right_edge = u.Quantity([coords['xmax'], - coords['ymax'], - coords['zmax']]).to(u.m) - g0 = physical_constants['gravity'] - gamma = 5./3. - dummy = np.zeros(rho.shape) - simulation_parameters = gdf.SimulationParameters([ - ['boundary_conditions', np.zeros(6) + 2], - ['cosmological_simulation', 0 ], - ['current_iteration', 0 ], - ['current_time', 0.0 ], - ['dimensionality', 3 ], - ['domain_dimensions', grid_dimensions ], - ['domain_left_edge', left_edge ], - ['domain_right_edge', right_edge ], - ['eta', 0.0 ], - ['field_ordering', 0 ], - ['gamma', gamma ], - ['gravity0', 0.0 ], - ['gravity1', 0.0 ], - ['gravity2', g0 ], - ['nu', 0.0 ], - ['num_ghost_zones', 0 ], - ['refine_by', 0 ], - ['unique_identifier', 'sacgdf2014' ] - ]) - + grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] + left_edge = u.Quantity([coords['xmin'], + coords['ymin'], + coords['zmin']]).to(u.m) + right_edge = u.Quantity([coords['xmax'], + coords['ymax'], + coords['zmax']]).to(u.m) + + gamma = 5./3. + dummy = np.zeros(rho.shape) + simulation_parameters = gdf.SimulationParameters([ + ['boundary_conditions', np.zeros(6) + 2], + ['cosmological_simulation', 0 ], + ['current_iteration', 0 ], + ['current_time', 0.0 ], + ['dimensionality', 3 ], + ['domain_dimensions', grid_dimensions ], + ['domain_left_edge', left_edge ], + ['domain_right_edge', right_edge ], + ['eta', 0.0 ], + ['field_ordering', 0 ], + ['gamma', gamma ], + ['gravity0', 0.0 ], + ['gravity1', 0.0 ], + ['gravity2', physical_constants['gravity'] ], + ['nu', 0.0 ], + ['num_ghost_zones', 0 ], + ['refine_by', 0 ], + ['unique_identifier', 'sacgdf2014' ] + ]) + if collective: + from mpi4py import MPI + gdf_file = gdf.create_file(h5py.File(filename,'w'), + simulation_parameters, grid_dimensions, + driver='mpio', comm=MPI.COMM_WORLD) + else: gdf_file = gdf.create_file(h5py.File(filename,'w'), simulation_parameters, grid_dimensions) - - gdf.write_field(gdf_file, rho, - 'density_bg', - 'Background Density' - ) - gdf.write_field(gdf_file, dummy*u.Unit('kg/m**3'), - 'density_pert', - 'Perturbation Density' - ) - gdf.write_field(gdf_file, - energy, - 'internal_energy_bg', - 'Background Internal Energy' - ) - gdf.write_field(gdf_file, dummy*u.Unit('kg/(m s**2)'), - 'internal_energy_pert', - 'Perturbation Internal Energy' - ) - gdf.write_field(gdf_file, Bx, - 'mag_field_x_bg', - 'x Component of Background Magnetic Field' - ) - gdf.write_field(gdf_file, dummy*u.Unit('T'), - 'mag_field_x_pert', - 'x Component of Pertubation Magnetic Field' - ) - gdf.write_field(gdf_file, By, - 'mag_field_y_bg', - 'y Component of Background Magnetic Field' - ) - gdf.write_field(gdf_file, dummy*u.Unit('T'), - 'mag_field_y_pert', - 'y Component of Pertubation Magnetic Field' - ) - gdf.write_field(gdf_file, Bz, - 'mag_field_z_bg', - 'z Component of Background Magnetic Field' - ) - gdf.write_field(gdf_file, dummy*u.Unit('T'), - 'mag_field_z_pert', - 'z Component of Pertubation Magnetic Field' - ) - gdf.write_field(gdf_file, dummy*u.Unit('m/s'), - 'velocity_x', - 'x Component of Velocity' - ) - gdf.write_field(gdf_file, dummy*u.Unit('m/s'), - 'velocity_y', - 'y Component of Velocity' - ) - gdf.write_field(gdf_file, dummy*u.Unit('m/s'), - 'velocity_z', - 'z Component of Velocity' - ) - - gdf_file.close() + + gdf.write_field(gdf_file, rho, + 'density_bg', + 'Background Density' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('kg/m**3'), + 'density_pert', + 'Perturbation Density' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + energy, + 'internal_energy_bg', + 'Background Internal Energy' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('kg/(m s**2)'), + 'internal_energy_pert', + 'Perturbation Internal Energy' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, Bx, + 'mag_field_x_bg', + 'x Component of Background Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('T'), + 'mag_field_x_pert', + 'x Component of Pertubation Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, By, + 'mag_field_y_bg', + 'y Component of Background Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('T'), + 'mag_field_y_pert', + 'y Component of Pertubation Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, Bz, + 'mag_field_z_bg', + 'z Component of Background Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('T'), + 'mag_field_z_pert', + 'z Component of Pertubation Magnetic Field' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('m/s'), + 'velocity_x', + 'x Component of Velocity' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('m/s'), + 'velocity_y', + 'y Component of Velocity' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, dummy*u.Unit('m/s'), + 'velocity_z', + 'z Component of Velocity' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + + gdf_file.close() + if collective: + mpi_save_SACvariables( + filename, + ((rho, 'density_bg' ), + (energy,'internal_energy_bg'), + (Bx, 'mag_field_x_bg' ), + (By, 'mag_field_y_bg' ), + (Bz, 'mag_field_z_bg' )), + xindex, + ) #============================================================================ @@ -161,77 +194,80 @@ def save_SACsources( option_pars, physical_constants, coords, - Nxyz + Nxyz, + xindex, + rank=0, + collective=False, ): """ Save the balancing forces for a SAC model with multiple flux tubes in hdf5 (gdf default) format after collating the data from mpi sub processes if necessary. """ - rank = 0 - if option_pars['l_mpi']: - from mpi4py import MPI - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - gather_vars = [ - Fx,Fy - ] - concat_vars = [] - for var in gather_vars: - concat_vars.append(comm.gather(var, root=0)) - - if rank == 0: - out_vars = [] - for cvar in concat_vars: - out_vars.append(np.concatenate(cvar, axis=0)) - - Fx,Fy = out_vars if rank == 0: print'writing',sourcesfile print'SAC background source terms' - grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] - left_edge = u.Quantity([coords['xmin'], - coords['ymin'], - coords['zmin']]).to(u.m) - right_edge = u.Quantity([coords['xmax'], - coords['ymax'], - coords['zmax']]).to(u.m) - g0 = physical_constants['gravity'] - - simulation_parameters = gdf.SimulationParameters([ - ['boundary_conditions', np.zeros(6) + 2], - ['cosmological_simulation', 0 ], - ['current_iteration', 0 ], - ['current_time', 0.0 ], - ['dimensionality', 3 ], - ['domain_dimensions', grid_dimensions ], - ['domain_left_edge', left_edge ], - ['domain_right_edge', right_edge ], - ['eta', 0.0 ], - ['field_ordering', 0 ], - ['gamma', physical_constants['gamma'] ], - ['gravity0', 0.0 ], - ['gravity1', 0.0 ], - ['gravity2', physical_constants['gravity']], - ['nu', 0.0 ], - ['num_ghost_zones', 0 ], - ['refine_by', 0 ], - ['unique_identifier', 'sacgdf2014' ] - ]) - - gdf_file = gdf.create_file(h5py.File(sourcesfile,'w'), simulation_parameters, grid_dimensions) - - gdf.write_field(gdf_file, - Fx, - 'balancing_force_x_bg', - 'x Component of Background Balancing Force' - ) - gdf.write_field(gdf_file, - Fy, - 'balancing_force_y_bg', - 'y Component of Background Balancing Force' + grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] + left_edge = u.Quantity([coords['xmin'], + coords['ymin'], + coords['zmin']]).to(u.m) + right_edge = u.Quantity([coords['xmax'], + coords['ymax'], + coords['zmax']]).to(u.m) + + simulation_parameters = gdf.SimulationParameters([ + ['boundary_conditions', np.zeros(6) + 2], + ['cosmological_simulation', 0 ], + ['current_iteration', 0 ], + ['current_time', 0.0 ], + ['dimensionality', 3 ], + ['domain_dimensions', grid_dimensions ], + ['domain_left_edge', left_edge ], + ['domain_right_edge', right_edge ], + ['eta', 0.0 ], + ['field_ordering', 0 ], + ['gamma', physical_constants['gamma'] ], + ['gravity0', 0.0 ], + ['gravity1', 0.0 ], + ['gravity2', physical_constants['gravity']], + ['nu', 0.0 ], + ['num_ghost_zones', 0 ], + ['refine_by', 0 ], + ['unique_identifier', 'sacgdf2014' ] + ]) + if collective: + from mpi4py import MPI + gdf_file = gdf.create_file(h5py.File(sourcesfile,'w'), + simulation_parameters, grid_dimensions, + driver='mpio', comm=MPI.COMM_WORLD) + else: + gdf_file = gdf.create_file(h5py.File(sourcesfile,'w'), + simulation_parameters, grid_dimensions) + + gdf.write_field(gdf_file, + Fx, + 'balancing_force_x_bg', + 'x Component of Background Balancing Force' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + Fy, + 'balancing_force_y_bg', + 'y Component of Background Balancing Force' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf_file.close() + if collective: + mpi_save_SACvariables( + sourcesfile, + ((Fx,'balancing_force_x_bg'), + (Fy,'balancing_force_y_bg')), + xindex, ) - gdf_file.close() #============================================================================ @@ -248,115 +284,143 @@ def save_auxilliary3D( option_pars, physical_constants, coords, - Nxyz + Nxyz, + xindex, + rank=0, + collective=False, ): """ Save auxilliary variables for use in plotting background setup in hdf5 (gdf default) format after collating the data from mpi sub processes if necessary. """ - rank = 0 - if option_pars['l_mpi']: - from mpi4py import MPI - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - gather_vars = [ - pressure_m, - rho_m, - temperature, - pbeta, - alfven, - cspeed, - dxB2, - dyB2 - ] - concat_vars = [] - for var in gather_vars: - concat_vars.append(comm.gather(var, root=0)) - - if rank == 0: - out_vars = [] - for cvar in concat_vars: - out_vars.append(np.concatenate(cvar, axis=0)) - - pressure_m, rho_m, temperature, pbeta,\ - alfven, cspeed, dxB2, dyB2 = out_vars if rank == 0: print'writing',auxfile print'non-SAC 3D auxilliary data for plotting' - grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] - left_edge = u.Quantity([coords['xmin'], - coords['ymin'], - coords['zmin']]).to(u.m) - right_edge = u.Quantity([coords['xmax'], - coords['ymax'], - coords['zmax']]).to(u.m) - - - simulation_parameters = gdf.SimulationParameters([ - ['boundary_conditions', np.zeros(6) + 2], - ['cosmological_simulation', 0 ], - ['current_iteration', 0 ], - ['current_time', 0.0 ], - ['dimensionality', 3 ], - ['domain_dimensions', grid_dimensions ], - ['domain_left_edge', left_edge ], - ['domain_right_edge', right_edge ], - ['eta', 0.0 ], - ['field_ordering', 0 ], - ['gamma', physical_constants['gamma'] ], - ['gravity0', 0.0 ], - ['gravity1', 0.0 ], - ['gravity2', physical_constants['gravity']], - ['nu', 0.0 ], - ['num_ghost_zones', 0 ], - ['refine_by', 0 ], - ['unique_identifier', 'sacgdf2014' ] - ]) - - gdf_file = gdf.create_file(h5py.File(auxfile,'w'), simulation_parameters, grid_dimensions) + grid_dimensions = [Nxyz[0], Nxyz[1], Nxyz[2]] + left_edge = u.Quantity([coords['xmin'], + coords['ymin'], + coords['zmin']]).to(u.m) + right_edge = u.Quantity([coords['xmax'], + coords['ymax'], + coords['zmax']]).to(u.m) + - gdf.write_field(gdf_file, - pressure_m, - 'pressure_mhs', - 'Background magneto-pressure balance' - ) - gdf.write_field(gdf_file, - rho_m, - 'density_mhs', - 'Background magneto-density balance' - ) - gdf.write_field(gdf_file, - temperature.to(u.K), - 'temperature', - 'Background temperature' - ) - gdf.write_field(gdf_file, - pbeta, - 'plasma_beta', - 'Background plasma beta' - ) - gdf.write_field(gdf_file, - alfven, - 'alfven_speed', - 'Background Alfven speed' - ) - gdf.write_field(gdf_file, - cspeed, - 'sound_speed', - 'Background sound speed' - ) - gdf.write_field(gdf_file, - dxB2, - 'mag_tension_x', - 'x-component background magnetic tension' - ) - gdf.write_field(gdf_file, - dyB2, - 'mag_tension_y', - 'y-component background magnetic tension' + simulation_parameters = gdf.SimulationParameters([ + ['boundary_conditions', np.zeros(6) + 2], + ['cosmological_simulation', 0 ], + ['current_iteration', 0 ], + ['current_time', 0.0 ], + ['dimensionality', 3 ], + ['domain_dimensions', grid_dimensions ], + ['domain_left_edge', left_edge ], + ['domain_right_edge', right_edge ], + ['eta', 0.0 ], + ['field_ordering', 0 ], + ['gamma', physical_constants['gamma'] ], + ['gravity0', 0.0 ], + ['gravity1', 0.0 ], + ['gravity2', physical_constants['gravity']], + ['nu', 0.0 ], + ['num_ghost_zones', 0 ], + ['refine_by', 0 ], + ['unique_identifier', 'sacgdf2014' ] + ]) + if collective: + from mpi4py import MPI + gdf_file = gdf.create_file(h5py.File(auxfile,'w'), + simulation_parameters, grid_dimensions, + driver='mpio', comm=MPI.COMM_WORLD) + else: + gdf_file = gdf.create_file(h5py.File(auxfile,'w'), + simulation_parameters, grid_dimensions) + + gdf.write_field(gdf_file, + pressure_m, + 'pressure_mhs', + 'Background magneto-pressure balance' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + rho_m, + 'density_mhs', + 'Background magneto-density balance' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + temperature, + 'temperature', + 'Background temperature' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + pbeta, + 'plasma_beta', + 'Background plasma beta' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + alfven, + 'alfven_speed', + 'Background Alfven speed' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + cspeed, + 'sound_speed', + 'Background sound speed' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + dxB2, + 'mag_tension_x', + 'x-component background magnetic tension' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf.write_field(gdf_file, + dyB2, + 'mag_tension_y', + 'y-component background magnetic tension' + ,field_shape=grid_dimensions + ,arr_slice=xindex + ,collective=collective + ) + gdf_file.close() + if collective: + mpi_save_SACvariables( + auxfile, + ((pressure_m, + 'pressure_mhs'), + ( rho_m, + 'density_mhs'), + (temperature, + 'temperature'), + (pbeta, + 'plasma_beta'), + (alfven, + 'alfven_speed'), + (cspeed, + 'sound_speed'), + (dxB2, + 'mag_tension_x'), + (dyB2, + 'mag_tension_y')), + xindex, ) - gdf_file.close() #============================================================================ @@ -368,40 +432,36 @@ def save_auxilliary1D( option_pars, physical_constants, coords, - Nxyz + Nxyz, + rank=0, + collective=False, ): """ Save auxilliary variables for use in plotting background setup in hdf5 (gdf default) format after collating the data from mpi sub processes if necessary. """ - rank = 0 - if option_pars['l_mpi']: - from mpi4py import MPI - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - if rank == 0: print'writing',auxfile print'non-SAC 1D auxilliary data for plotting' - grid_dimensions = [2, 2, Nxyz[2]] #dims > 1 to be read by yt - left_edge = u.Quantity([coords['xmin'], - coords['ymin'], - coords['zmin']]).to(u.m) - right_edge = u.Quantity([coords['xmax'], - coords['ymax'], - coords['zmax']]).to(u.m) - - pressureHS = u.Quantity(np.zeros(grid_dimensions), - unit=pressure_Z.to("Pa").unit) - rhoHS = u.Quantity(np.zeros(grid_dimensions), - unit=rho_Z.to('kg/m**3').unit) - RgasHS = u.Quantity(np.zeros(grid_dimensions), - unit=Rgas_Z.to('m**2/(K s**2)').unit) - pressureHS[:] = pressure_Z - rhoHS[:] = rho_Z - RgasHS[:] = Rgas_Z - simulation_parameters = gdf.SimulationParameters([ + grid_dimensions = [2, 2, Nxyz[2]] #dims > 1 to be read by yt + left_edge = u.Quantity([coords['xmin'], + coords['ymin'], + coords['zmin']]).to(u.m) + right_edge = u.Quantity([coords['xmax'], + coords['ymax'], + coords['zmax']]).to(u.m) + + pressureHS = u.Quantity(np.zeros(grid_dimensions), + unit=pressure_Z.to("Pa").unit) + rhoHS = u.Quantity(np.zeros(grid_dimensions), + unit=rho_Z.to('kg/m**3').unit) + RgasHS = u.Quantity(np.zeros(grid_dimensions), + unit=Rgas_Z.to('m**2/(K s**2)').unit) + pressureHS[:] = pressure_Z + rhoHS[:] = rho_Z + RgasHS[:] = Rgas_Z + simulation_parameters = gdf.SimulationParameters([ ['boundary_conditions', np.zeros(6) + 2], ['cosmological_simulation', 0 ], ['current_iteration', 0 ], @@ -422,25 +482,54 @@ def save_auxilliary1D( ['unique_identifier', 'sacgdf2014' ] ]) -# import pdb; pdb.set_trace() - gdf_file = gdf.create_file(h5py.File(auxfile,'w'), simulation_parameters, grid_dimensions) - - gdf.write_field(gdf_file, - pressureHS, - 'pressure_HS', - 'Background 1D hydrostatic-pressure' - ) - gdf.write_field(gdf_file, - rhoHS, - 'density_HS', - 'Background 1D hydrostatic-density' - ) - gdf.write_field(gdf_file, - RgasHS, - 'ideal_gas_constant_HS', - 'Background 1D hydrostatic-R_gas' - ) - gdf_file.close() +# import pdb; pdb.set_trace() + gdf_file = gdf.create_file(h5py.File(auxfile,'w'), simulation_parameters, grid_dimensions) + + gdf.write_field(gdf_file, + pressureHS, + 'pressure_HS', + 'Background 1D hydrostatic-pressure' + ,collective=collective + ) + gdf.write_field(gdf_file, + rhoHS, + 'density_HS', + 'Background 1D hydrostatic-density' + ,collective=collective + ) + gdf.write_field(gdf_file, + RgasHS, + 'ideal_gas_constant_HS', + 'Background 1D hydrostatic-R_gas' + ,collective=collective + ) + gdf_file.close() #============================================================================= +## Save a file!!! +##============================================================================ +""" For large data arrays this has been producing overfull memory so moved to +dedicated serial function, which should be parallel and moved ahead of gather +for plotting -- maybe handle plotting from hdf5 also + +This file is potentially large - recommended to mkdir hdf5 in /data/${USER} +and add symlink to ${HOME} to avoid exceeding quota. +""" +def mpi_save_SACvariables( + filename, + fields, + xindex, + ): + + """ Save the background variables for a SAC model in hdf5 (gdf default) + format after collating the data from mpi sub processes if necessary. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + + ds = h5py.File(filename, 'a', driver='mpio', comm=comm) + for var, field in fields: + ds['data']['grid_0000000000'][field][xindex,...] = var + ds.close()