Codecomparison

build/Makefile

@@ -253,13 +253,14 @@ endif
 
 ifeq ($(SETUP), disc)
 #   locally isothermal gas disc
-    DISC_VISCOSITY=yes
+    DISC_VISCOSITY=no
     SETUPFILE= setup_disc.f90
     ANALYSIS= analysis_disc.f90
     ISOTHERMAL=yes
     KNOWN_SETUP=yes
     MULTIRUNFILE= multirun.f90
     IND_TIMESTEPS=yes
+    MODDUMP=moddump_removeparticles_radius.f90
 endif
 
 ifeq ($(SETUP), grtde)
@@ -1884,7 +1885,8 @@ SRCSETUP= utils_omp.F90 utils_sort.f90 utils_timing.f90 utils_summary.F90 utils_
           mpi_balance.F90 set_dust_options.f90 \
           utils_indtimesteps.F90 partinject.F90 stack.F90 mpi_dens.F90 mpi_force.F90 mpi_derivs.F90 \
           ${SRCTURB} ${SRCNIMHD} ${SRCCHEM} \
-          ptmass.F90 energies.F90 density_profiles.f90 set_slab.f90 set_disc.F90 relax_star.f90 \
+          ptmass.F90 energies.F90 density_profiles.f90 set_slab.f90 load_from_file.f90 interpolate_grid.f90 \
+          grids_for_setup.f90 set_disc.F90 relax_star.f90 \
 	  set_cubic_core.f90 set_fixedentropycore.f90 set_softened_core.f90 \
           set_vfield.f90 sort_particles.F90 dust_formation.F90 ptmass_radiation.F90 ${SRCINJECT} \
           ${SETUPFILE} checksetup.F90 \
@@ -1912,7 +1914,8 @@ SRCTESTS=utils_testsuite.f90 ${TEST_FASTMATH} test_kernel.f90 \
          test_dust.F90 test_growth.F90 test_smol.F90 \
          test_nonidealmhd.F90 directsum.f90 test_gravity.F90 \
          test_derivs.F90 test_cooling.f90 test_eos.f90 test_externf.f90 test_rwdump.f90 \
-         test_step.F90 test_indtstep.F90 set_disc.F90 test_setdisc.F90 \
+         test_step.F90 test_indtstep.F90 load_from_file.f90 interpolate_grid.f90 \
+         grids_for_setup.f90 set_disc.F90 test_setdisc.F90 \
          test_link.F90 test_kdtree.F90 test_part.f90 test_ptmass.F90 test_luminosity.F90\
          test_gnewton.f90 test_corotate.f90 test_geometry.f90 \
          test_sedov.F90 test_radiation.F90 \
@@ -2020,7 +2023,8 @@ MODDUMPBIN=phantommoddump
 endif
 OBJMOD1 = utils_omp.F90 utils_summary.f90 utils_indtimesteps.F90 \
           utils_sort.f90 set_Bfield.f90 \
-          random.f90 set_disc.F90 set_dust.f90 \
+          random.f90 load_from_file.f90 interpolate_grid.f90 \
+          grids_for_setup.f90 set_disc.F90 set_dust.f90 \
           ${SRCTURB} ${SRCNIMHD} ${SRCCHEM} dust_formation.F90 ptmass_radiation.F90 \
           icosahedron.f90 radiation_utils.f90 cons2prim.f90 \
           partinject.F90 utils_inject.f90 ${SRCINJECT} \

scripts/sinkanalysis.py

@@ -16,6 +16,59 @@
 import sys
 import fnmatch
 
+def loadevfile():
+        #find sink files in directory
+        listFiles1=[]
+        files=os.listdir("./") #can alaso use os.listdir("./")
+        for evfiles in files:
+                if fnmatch.fnmatch(evfiles, '*[!N][0-9][0-9].ev'):
+                        listFiles1.append(evfiles)
+
+#        filename= sys.argv[j]+"Sink"+
+#       print "Prefix: ",filename
+
+        ######################################################################
+        #--------------------------------------------
+        #    THIS IS JUST A HACK TO COMPUTE NUMBER OF COLUMNS
+        #    Get column number and header lines 
+        #    Print the number of columns for each line in a prov file.
+        #    Number of columns given when ncols is equal for 4 lines.
+        #    Assuming that the header lenght is less than 20 lines
+        #----------------------------------------------
+
+        listFiles1.sort()
+        for filename in listFiles1:
+                print(filename)
+                ncols = 0
+                nrowh = 0
+                n = 21
+
+                command1="awk 'NR<"+ str(n)+" {print NF}' "+filename+">prov"
+                os.system(command1)
+
+                prov=np.loadtxt("prov")
+
+                for i in range(n-4):
+                    if (prov[i]==ncols and prov[i+1]==ncols and prov[i+2]==ncols and prov[i+3]==ncols): 
+                        nrowh=i+1
+                        break
+                    ncols=prov[i]
+
+                ncols=prov[i]
+                os.system("rm prov")  
+                skiph=nrowh-3
+        #########################################################
+
+                if "data1" in locals():
+                        dataProv=np.loadtxt(filename,skiprows=int(skiph))
+                        data1=np.vstack((data1,dataProv))
+                else:
+                        data1=np.loadtxt(filename,skiprows=int(skiph))
+
+        return data1
+
+
+
 def loadSink():
         #find sink files in directory
         listFiles1=[]
@@ -42,7 +95,7 @@ def loadSink():
         listFiles1.sort()
         listFiles2.sort()
         for filename in listFiles1:
-                print filename
+                print(filename)
                 ncols = 0
                 nrowh = 0
                 n = 21
@@ -60,16 +113,17 @@ def loadSink():
 
                 ncols=prov[i]
                 os.system("rm prov")  
+                skiph=nrowh-3
         #########################################################
 
                 if "data1" in locals():
-                        dataProv=np.loadtxt(filename,skiprows=int(nrowh))
+                        dataProv=np.loadtxt(filename,skiprows=int(skiph))
                         data1=np.vstack((data1,dataProv))
                 else:
-                        data1=np.loadtxt(filename,skiprows=int(nrowh))
+                        data1=np.loadtxt(filename,skiprows=int(skiph))
 
         for filename in listFiles2:
-                print filename
+                print(filename)
                 ncols = 0
                 nrowh = 0
                 n = 21
@@ -89,13 +143,14 @@ def loadSink():
 
                 ncols=prov[i]
                 os.system("rm prov")  
+                skiph=nrowh-3
         ###################################
 
                 if "data2" in locals():
-                        dataProv=np.loadtxt(filename,skiprows=int(nrowh))
+                        dataProv=np.loadtxt(filename,skiprows=int(skiph))
                         data2=np.vstack((data2,dataProv))
                 else:
-                        data2=np.loadtxt(filename,skiprows=int(nrowh))
+                        data2=np.loadtxt(filename,skiprows=int(skiph))
 
 
         return data1,data2
@@ -116,6 +171,12 @@ def FindOrbEvo(sink1,sink2):
         vy2=sink2[:,6]
         vz2=sink2[:,7]
         M2=sink2[:,4]
+        f1=sink1[:,12:15]
+        f2=sink2[:,12:15]
+        f1R1=sink1[:,20:23] # note that 23 is beyon the last element but syntax require 20:23 or 20: for loading 
+        f2R1=sink2[:,20:23] # the last 3 columns.
+
+
 
         xx1=np.array([x1,y1,z1])
         xx2=np.array([x2,y2,z2])
@@ -155,11 +216,15 @@ def FindOrbEvo(sink1,sink2):
 
 
         j1=x1*vy1-y1*vx1
+        Tor1=(x1*f1[:,1]-y1*f1[:,0])*M1
+        Tor1R1=(x1*f1R1[:,1]-y1*f1R1[:,0])*M1
         R1=np.sqrt(x1**2+y1**2+z1**2)
         v1Tot=np.sqrt(vx1**2+vy1**2+vz1**2)
 
 
         j2=x2*vy2-y2*vx2
+        Tor2=(x2*f2[:,1]-y2*f2[:,0])*M2
+        Tor2R1=(x2*f2R1[:,1]-y2*f2R1[:,0])*M2
         R2=np.sqrt(x2**2+y2**2+z2**2)
         v2Tot=np.sqrt(vx2**2+vy2**2+vz2**2)
 
@@ -175,55 +240,64 @@ def FindOrbEvo(sink1,sink2):
         eccConf=np.sqrt(1-j1**2/((M1+M2)*a1**4)*a**3)
         #needs to be computed in the CM centred in the focus of the orbit (i.e. M_*)
         Phase=np.arctan2(-(j2S*vx2S/(M1+M2)+y2S/R2S),j2S*vy2S/(M1+M2)-x2S/R2S)+np.pi
+        Torque=Tor1+Tor2
+        TorqueR1=Tor1R1+Tor2R1
 
-        return Time,ecc,a,Phase
+        return Time,ecc,a,Phase,Torque,TorqueR1
 
 if __name__=="__main__":
 
         drawdirect=True
         sink1,sink2=loadSink()
-        Time,ecc,a,Phase=FindOrbEvo(sink1,sink2)
+        Time,ecc,a,Phase,Torque,TorqueR1=FindOrbEvo(sink1,sink2)
 
+        if(drawdirect):
 
-        plt.figure(1)
-        plt.plot(Time,ecc)
-        plt.title("ecc vs time")
-        plt.xlabel("$t$")
-        plt.ylabel("$e$")
-
-        plt.figure(2)
-        plt.plot(Time,a)
-        plt.title("a vs time")
-        plt.xlabel("$t$")
-        plt.ylabel("$a$")
+            plt.figure(1)
+            plt.plot(Time,ecc)
+            plt.title("ecc vs time")
+            plt.xlabel("$t$")
+            plt.ylabel("$e$")
+
+            plt.figure(2)
+            plt.plot(Time,a)
+            plt.title("a vs time")
+            plt.xlabel("$t$")
+            plt.ylabel("$a$")
+
+            plt.figure(3)
+            plt.plot(Time,Phase,marker=".",linestyle="")
+            plt.title("Phase vs time")
+            plt.xlabel("$t$")
+            plt.ylabel("$\\Phi$")
+
+            plt.figure(4)
+            plt.plot(Time,Torque,marker=".",linestyle="")
+            plt.title("Torque vs time")
+            plt.xlabel("$t$")
+            plt.ylabel("$(\\bar R x\\bar F)_z$")
 
-        plt.figure(3)
-        plt.plot(Time,Phase,marker=".",linestyle="")
-        plt.title("Phase vs time")
-        plt.xlabel("$t$")
-        plt.ylabel("$\\Phi$")
 
 
-        if(drawdirect):
             plt.draw()
             plt.pause(1)
 
-            eccen_name = raw_input("Enter file name to save eccen vs time plot (blank = don't save): ")
+            eccen_name = input("Enter file name to save eccen vs time plot (blank = don't save): ")
             if (eccen_name.strip() != ''):
                 plt.figure(1)
                 plt.savefig(eccen_name.strip())
 
-            semi_name = raw_input("Enter file name to save a vs time plot (blank = don't save): ")
+            semi_name = input("Enter file name to save a vs time plot (blank = don't save): ")
             if (semi_name.strip() != ''):
                 plt.figure(2)
                 plt.savefig(semi_name.strip())
 
-            phi_name = raw_input("Enter file name to save phi vs time plot (blank = don't save): ")
+            phi_name = input("Enter file name to save phi vs time plot (blank = don't save): ")
             if (phi_name.strip() != ''):
                 plt.figure(3)
                 plt.savefig(phi_name.strip())
 
-            raw_input("<Hit enter to close the plots>")
+            input("<Hit enter to close the plots>")
             plt.close('all')
 
 

src/main/cons2prim.f90

@@ -210,7 +210,6 @@ subroutine cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,&
        yi      = xyzh(2,i)
        zi      = xyzh(3,i)
        hi      = xyzh(4,i)
-
        if (maxphase==maxp) call get_partinfo(iphase(i),iactivei,iamgasi,iamdusti,iamtypei)
 
        pmassi  = massoftype(iamtypei)
@@ -239,6 +238,13 @@ subroutine cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,&
           else
              call equationofstate(ieos,p_on_rhogas,spsound,rhogas,xi,yi,zi,eni=vxyzu(4,i),tempi=temperaturei)
           endif
+
+       elseif (ieos == 20) then
+          !eccentric isothermal
+          call equationofstate(ieos,p_on_rhogas,spsound,rhogas,&
+                                xi,yi,zi,&
+                                vxi=vxyzu(1,i),vyi=vxyzu(2,i),vzi=vxyzu(3,i),&
+                                tempi=temperaturei)
        else
           !isothermal
           call equationofstate(ieos,p_on_rhogas,spsound,rhogas,xi,yi,zi,tempi=temperaturei)