Monaghan viscosity for Gadget2

src/amuse/community/gadget2/interface.cc

@@ -62,6 +62,7 @@ void set_default_parameters(){
     All.DesNumNgb = 50;
     All.MaxNumNgbDeviation = 5.;
     All.ArtBulkViscConst = 0.5;
+    All.ArtBulkViscBeta = 1.;
     All.MinGasTemp = 0;
     All.UnitLength_in_cm = 3.085678e21;
     All.UnitMass_in_g = 1.989e43;
@@ -1007,6 +1008,15 @@ int set_alpha(double artificial_viscosity_alpha){
     All.ArtBulkViscConst = artificial_viscosity_alpha;
     return 0;
 }
+int get_beta(double *artificial_viscosity_beta){
+    if (ThisTask) {return 0;}
+    *artificial_viscosity_beta = All.ArtBulkViscBeta;
+    return 0;
+}
+int set_beta(double artificial_viscosity_beta){
+    All.ArtBulkViscBeta = artificial_viscosity_beta;
+    return 0;
+}
 int get_courant(double *courant){
     if (ThisTask) {return 0;}
     *courant = All.CourantFac*2.0;

src/amuse/community/gadget2/interface.py

@@ -692,7 +692,23 @@ def get_alpha():
         function.addParameter('alpha', dtype='d', direction=function.OUT)
         function.result_type = 'i'
         return function
+
+    @legacy_function
+    def set_beta():
+        """ SPH artificial viscosity beta parameter (2*viscosity alpha parameter) """        
+        function = LegacyFunctionSpecification()
+        function.addParameter('beta', dtype='d', direction=function.IN)
+        function.result_type = 'i'
+        return function
 
+    @legacy_function
+    def get_beta():
+        """ SPH artificial viscosity beta parameter (2*viscosity alpha parameter) """        
+        function = LegacyFunctionSpecification()
+        function.addParameter('beta', dtype='d', direction=function.OUT)
+        function.result_type = 'i'
+        return function   
+
     @legacy_function
     def set_courant():
         """ SPH courant condition parameter (0.3) """            
@@ -1394,7 +1410,16 @@ def define_parameters(self, object):
             "SPH artificial viscosity alpha parameter (0.5)", 
             default_value = 0.5
         )
-
+
+
+        object.add_method_parameter(
+            "get_beta", 
+            "set_beta",
+            "artificial_viscosity_beta", 
+            "SPH artificial viscosity beta parameter (2*viscosity alpha parameter)", 
+            default_value = 1.0
+        )
+
         object.add_method_parameter(
             "get_courant", 
             "set_courant",
@@ -2112,7 +2137,21 @@ def define_methods(self, object):
             (object.NO_UNIT, ),
             (object.ERROR_CODE,)
         )
+
+
+        object.add_method(
+            "get_beta",
+            (),
+            (object.NO_UNIT, object.ERROR_CODE,)
+        )
 
+        object.add_method(
+            "set_beta",
+            (object.NO_UNIT, ),
+            (object.ERROR_CODE,)
+        )
+
+
         object.add_method(
             "get_courant",
             (),

src/amuse/community/gadget2/src/allvars.h

@@ -287,6 +287,7 @@ extern struct global_data_all_processes
   double MaxNumNgbDeviation;    /*!< Maximum allowed deviation neighbour number */
 
   double ArtBulkViscConst;      /*!< Sets the parameter \f$\alpha\f$ of the artificial viscosity */
+  double ArtBulkViscBeta;      /*!< Sets the parameter \f$\beta\f$ of the artificial viscosity */
   double InitGasTemp;		/*!< may be used to set the temperature in the IC's */
   double MinGasTemp;		/*!< may be used to set a floor for the gas temperature */
   double MinEgySpec;            /*!< the minimum allowed temperature expressed as energy per unit mass */

src/amuse/community/gadget2/src/begrun.c

@@ -103,6 +103,7 @@ void begrun(void)
       All.NumFilesPerSnapshot = all.NumFilesPerSnapshot;
       All.MaxNumNgbDeviation = all.MaxNumNgbDeviation;
       All.ArtBulkViscConst = all.ArtBulkViscConst;
+      All.ArtBulkViscBeta = 2*all.ArtBulkViscConst;
 
 
       All.OutputListOn = all.OutputListOn;
@@ -682,7 +683,7 @@ void read_parameter_file(char *fname)
   if(errorFlag)
     {
 #ifndef NOMPI
-      MPI_Finalize();
+      MPI_Finalize();
 #endif
       exit(0);
     }

src/amuse/community/gadget2/src/hydra.c

@@ -1,9 +1,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include <math.h>
+#include <math.h>
 #ifndef NOMPI
-#include <mpi.h>
+#include <mpi.h>
 #endif
 #include <gsl/gsl_math.h>
 #include "allvars.h"
@@ -250,9 +250,9 @@ void hydro_force(void)
 	  timecomp += timediff(tstart, tend);
 
 	  /* do a block to measure imbalance */
-	  tstart = second();
+	  tstart = second();
 #ifndef NOMPI
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(MPI_COMM_WORLD);
 #endif
 	  tend = second();
 	  timeimbalance += timediff(tstart, tend);
@@ -389,6 +389,9 @@ void hydro_evaluate(int target, int mode)
   double dx, dy, dz, dvx, dvy, dvz;
   double h_i2, hinv, hinv4;
   double p_over_rho2_i, p_over_rho2_j, soundspeed_i, soundspeed_j;
+#ifdef MONAGHANVISC
+  double soundspeed_ij, h_ij;
+#endif
   double hfc, dwk_i, vdotr, vdotr2, visc, mu_ij, rho_ij, vsig;
   double h_j, dwk_j, r, r2, u, hfc_visc;
 
@@ -523,19 +526,32 @@ void hydro_evaluate(int target, int mode)
 
 		  if(vdotr2 < 0)	/* ... artificial viscosity */
 		    {
+#ifndef MONAGHANVISC
 		      mu_ij = fac_mu * vdotr2 / r;	/* note: this is negative! */
-
+#else		      
+		      h_ij = 0.5 * (h_i + h_j);
+		      mu_ij = fac_mu * h_ij * vdotr2 / (r2 + 0.0001 * h_ij * h_ij);
+#endif
                       vsig = soundspeed_i + soundspeed_j - 3 * mu_ij;
 
+
+
 		      if(vsig > maxSignalVel)
 			maxSignalVel = vsig;
 
 		      rho_ij = 0.5 * (rho + SphP[j].Density);
+
 		      f2 =
 			fabs(SphP[j].DivVel) / (fabs(SphP[j].DivVel) + SphP[j].CurlVel +
 						0.0001 * soundspeed_j / fac_mu / SphP[j].Hsml);
 
+#ifndef MONAGHANVISC
 		      visc = 0.25 * All.ArtBulkViscConst * vsig * (-mu_ij) / rho_ij * (f1 + f2);
+#else			      
+		      soundspeed_ij = (soundspeed_i+soundspeed_j) * 0.5;
+
+		      visc = ((-All.ArtBulkViscConst) * soundspeed_ij * mu_ij + All.ArtBulkViscBeta * mu_ij * mu_ij) / rho_ij;
+#endif
 
 		      /* .... end artificial viscosity evaluation */
 #ifndef NOVISCOSITYLIMITER

src/amuse/community/gadget2/src/makefile_options

@@ -53,6 +53,7 @@ OPT   +=  -DNOSTOP_WHEN_BELOW_MINTIMESTEP
 #OPT   +=  -DISOTHERM_EQS
 #OPT   +=  -DADAPTIVE_GRAVSOFT_FORGAS
 #OPT   +=  -DSELECTIVE_NO_GRAVITY=2+4+8+16
+#OPT   +=  -DMONAGHANVISC
 
 #--------------------------------------- Testing and Debugging options
 #OPT   +=  -DFORCETEST=0.1

test/codes_tests/test_gadget2.py

@@ -485,7 +485,7 @@ def test7(self):
                 ('courant',0.3), ('type_of_timestep_criterion',0), ('hubble_parameter', 0.7),
                 ('omega_zero',0.0), ('omega_lambda',0.0), ('omega_baryon',0.0), 
                 ('min_gas_hsmooth_fractional',0.0), ('timestep_accuracy_parameter',0.025), 
-                ('tree_domain_update_frequency',0.05)]:
+                ('tree_domain_update_frequency',0.05),('artificial_viscosity_beta',1.0)]:
             self.assertEquals(value, getattr(instance.parameters, par))
             setattr(instance.parameters, par, 1)
             self.assertEquals(1, getattr(instance.parameters, par))