Coverage for gamdpy/calculators/calculator_hydrodynamic_profile.py: 89%

1import numpy as np 

2import gamdpy as gp 

3from numba import jit 

4 

5@jit(nopython=True) 

6def __dosample__(cpart, temp, momc, pos, vel, velT, mass, ptype, thistype, npart, wbin, hbox): 

7 

8 for n in range(npart): 

9 if ptype[n]==thistype: 

10 

11 idx = int((pos[n]+hbox)/wbin) 

12 

13 cpart[idx] = cpart[idx] + 1 

14 momc[idx] = momc[idx] + mass*vel[n] 

15 temp[idx] = temp[idx] + mass*(velT[n,0]*velT[n,0] + velT[n,1]*velT[n,1]) 

16 

17 

18class CalculatorHydrodynamicProfile: 

19 """ 

20 Calculates the density, streaming velocity, and kinetic temperature profiles.  

21  

22 Example: See atomistic_wall example 

23  

24 Initialization variables: 

25 - configuration: Instance of the configuration class 

26 - ptype: The particle type for which the profile is calculated 

27 - bins: Number of bins to use in the profiles; more bins higher resolution. (default 100) 

28 - profdir: Profile spatial direction (default 2 - or z - direction) 

29 - veldir: The streaming velocity component used (default 0 - or x - direction) 

30 - verbose: If true print some information (default True) 

31 

32 Output:  

33 With the method read() you can retrieve the current data. By default a data file is printed with the same  

34 """ 

35 

36 def __init__(self, configuration, ptype, bins=100, profdir=2, veldir=0, verbose=True): 

37 self.conf = configuration 

38 self.ptype = ptype 

39 self.bins = bins 

40 self.pdir = profdir 

41 self.vdir = veldir 

42 

43 self.cpart = np.zeros(bins, dtype=np.int64) 

44 self.dens = np.zeros(bins, dtype=np.float64) 

45 self.temp = np.zeros(bins, dtype=np.float64) 

46 self.momc = np.zeros(bins, dtype=np.float64) 

47 

48 self.volbin = configuration.get_volume()/bins 

49 self.widthbin = configuration.simbox.get_lengths()[profdir]/bins 

50 self.hbox = configuration.simbox.get_lengths()[profdir]*0.5 

51 self.nsample = 0 

52 

53 self.Tdir = [0, 1, 2] 

54 self.Tdir.remove(veldir) 

55 

56 for n in range(configuration.N): 

57 if configuration.ptype[n] == self.ptype: 

58 self.mass = configuration['m'][n] 

59 break 

60 

61 if verbose: 

62 print(f"Types {self.ptype}, no. bin {self.bins}, profile dir {self.pdir}, vel dir {self.vdir}") 

63 

64 def update(self): 

65 

66 __dosample__(self.cpart, self.temp, self.momc, 

67 self.conf['r'][:,self.pdir], self.conf['v'][:,self.vdir], self.conf['v'][:,self.Tdir], 

68 self.mass, self.conf.ptype, self.ptype, self.conf.N, self.widthbin, self.hbox) 

69 

70 self.nsample = self.nsample + 1 

71 

72 def read(self, save=True): 

73 

74 density = np.zeros(self.bins) 

75 svel = np.zeros(self.bins) 

76 temp = np.zeros(self.bins) 

77 x = np.zeros(self.bins) 

78 

79 fac = 1.0/(self.nsample*self.volbin) 

80 for n in range(self.bins): 

81 x[n] = (n+0.5)*self.widthbin 

82 if self.cpart[n] > 0: 

83 density[n] = self.mass*self.cpart[n]*fac 

84 svel[n] = self.momc[n]*fac/density[n] 

85 temp[n] = self.temp[n]/(2.0*self.cpart[n]) 

86 

87 if save: 

88 fp = open("HydrodynamicProfile.dat","w") 

89 for n in range(self.bins): 

90 fp.write("%f %f %f %f\n" % (x[n], density[n], svel[n], temp[n])) 

91 

92 fp.close() 

93 

94 return (x, density, svel, temp)