Add eccentric shearing sheet example simulation

examples/shearing_sheet/problem.c

@@ -103,5 +103,4 @@ void heartbeat(struct reb_simulation* const r){
     if (reb_simulation_output_check(r, 2.*M_PI/r->ri_sei.OMEGA)){
         //reb_simulation_output_ascii("position.txt");
     }
-}
-
+}

examples/shearing_sheet_eccentric/Makefile

@@ -0,0 +1,35 @@
+export OPENGL=1# Set this to 1 to enable OpenGL
+export SERVER=0# Set this to 1 to enable the visualization web server
+include ../../src/Makefile.defs
+
+# CCPROBLEM is defined in Makefile.defs to allow for
+# a compact cross platform Makefile
+.PHONY: all librebound
+all: problem.c librebound
+	@echo "Compiling $< ..."
+	$(CCPROBLEM)
+	@echo ""
+	@echo "Compilation successful. To run REBOUND, execute the file '$(EXEREBOUND)'."
+	@echo ""
+
+librebound:
+	@echo "Compiling shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/
+	@-$(RM) $(LIBREBOUND)
+	@$(LINKORCOPYLIBREBOUND)
+	@echo ""
+
+clean:
+	@echo "Cleaning up shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/ clean
+	@echo "Cleaning up local directory ..."
+	@-$(RM) $(LIBREBOUND)
+	@-$(RM) $(EXEREBOUND)
+
+rebound_webgl.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL enabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -DOPENGL=1 -sSTACK_SIZE=655360 -s USE_GLFW=3 -s FULL_ES3=1 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_webgl.html -o rebound_webgl.html
+
+rebound_console.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL disabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -sSTACK_SIZE=655360 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_console.html -o rebound_console.html

examples/shearing_sheet_eccentric/problem.c

@@ -0,0 +1,115 @@
+/**
+ * Shearing sheet (Hill's approximation)
+ *
+ * This example simulates a small patch of Saturn's
+ * Rings in shearing sheet coordinates. If you have OpenGL enabled, 
+ * you'll see one copy of the computational domain. Press `g` to see
+ * the ghost boxes which are used to calculate gravity and collisions.
+ * Particle properties resemble those found in Saturn's rings. 
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "rebound.h"
+
+double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v);
+void heartbeat(struct reb_simulation* const r);
+
+int main(int argc, char* argv[]) {
+    struct reb_simulation* r = reb_simulation_create();
+
+    // This allows you to connect to the simulation using
+    // a web browser. Simply go to http://localhost:1234
+    reb_simulation_start_server(r, 1234);
+
+    // Setup constants
+    r->opening_angle2    = .5;                  // This determines the precission of the tree code gravity calculation.
+    r->integrator        = REB_INTEGRATOR_SEI;
+    r->boundary          = REB_BOUNDARY_SHEAR_E;
+    r->gravity           = REB_GRAVITY_BASIC;
+    r->collision         = REB_COLLISION_DIRECT;
+    r->collision_resolve = reb_collision_resolve_hardsphere;
+    double OMEGA         = 0.00013143527;       // 1/s
+    double Q_NL          = 0.1;
+    r->ri_sei.OMEGA      = OMEGA;
+    r->ri_sei.Q_NL       = Q_NL;
+    r->G                 = 6.67428e-11;         // N / (1e-5 kg)^2 m^2
+    r->softening         = 0.1;                 // m
+    r->dt                = 1e-3*2.*M_PI/OMEGA;  // s
+    r->heartbeat         = heartbeat;           // function pointer for heartbeat
+    // This example uses two root boxes in the x and y direction. 
+    // Although not necessary in this case, it allows for the parallelization using MPI. 
+    // See Rein & Liu for a description of what a root box is in this context.
+    double surfacedensity          = 400;     // kg/m^2
+    double particle_density        = 400;     // kg/m^3
+    double particle_radius_min     = 1;       // m
+    double particle_radius_max     = 4;       // m
+    double particle_radius_slope   = -3;    
+    double boxsize                 = 50;         // m
+    if (argc>1){                              // Try to read boxsize from command line
+        boxsize = atof(argv[1]);
+    }
+    reb_simulation_configure_box(r, boxsize, 2, 2, 1);
+    r->N_ghost_x = 2;
+    r->N_ghost_y = 2;
+    r->N_ghost_z = 0;
+
+    // Initial conditions
+    printf("Toomre wavelength: %f\n",4.*M_PI*M_PI*surfacedensity/OMEGA/OMEGA*r->G);
+    // Use Bridges et al coefficient of restitution.
+    r->coefficient_of_restitution = coefficient_of_restitution_bridges;
+    // When two particles collide and the relative velocity is zero, the might sink into each other in the next time step.
+    // By adding a small repulsive velocity to each collision, we prevent this from happening.
+    r->minimum_collision_velocity = particle_radius_min*OMEGA*0.001;  // small fraction of the shear accross a particle
+
+
+    // Add all ring paricles
+    double total_mass = surfacedensity*r->boxsize.x*r->boxsize.y;
+    double mass = 0;
+    while(mass<total_mass){
+        struct reb_particle pt;
+        double x_0 = reb_random_uniform(r, -r->boxsize.x / 2., r->boxsize.x / 2.);
+        double y_0 = reb_random_uniform(r, -r->boxsize.y / 2., r->boxsize.y / 2.);
+
+        pt.x         = x_0-Q_NL*x_0*cos(OMEGA*r->t);
+        pt.y         = y_0-1.5*x_0*OMEGA*r->t+2.0*Q_NL*x_0*sin(OMEGA*r->t);
+        pt.z         = reb_random_normal(r, 1.);                    // m
+        pt.vx         = Q_NL*pt.x*OMEGA*sin(OMEGA*r->t);
+        pt.vy         = -1.5*pt.x*OMEGA+2.0*Q_NL*pt.x*OMEGA*cos(OMEGA*r->t);
+        pt.vz         = 0;
+        pt.ax         = 0;
+        pt.ay         = 0;
+        pt.az         = 0;
+        double radius     = reb_random_powerlaw(r, particle_radius_min,particle_radius_max,particle_radius_slope);
+        pt.r         = radius;                        // m
+        double        particle_mass = particle_density*4./3.*M_PI*radius*radius*radius;
+        pt.m         = particle_mass;     // kg
+        reb_simulation_add(r, pt);
+        mass += particle_mass;
+    }
+    reb_simulation_integrate(r, INFINITY);
+    reb_simulation_free(r);
+}
+
+// This example is using a custom velocity dependend coefficient of restitution
+double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v){
+    // assumes v in units of [m/s]
+    double eps = 0.32*pow(fabs(v)*100.,-0.234);
+    if (eps>1) eps=1;
+    if (eps<0) eps=0;
+    return eps;
+}
+
+void heartbeat(struct reb_simulation* const r){
+    if (reb_simulation_output_check(r, 1e-1*2.*M_PI/r->ri_sei.OMEGA)){
+        reb_simulation_output_timing(r, 0);
+        //reb_output_append_velocity_dispersion("veldisp.txt");
+    }
+    if (reb_simulation_output_check(r, 2.*M_PI/r->ri_sei.OMEGA)){
+        //reb_simulation_output_ascii("position.txt");
+    }
+}
+
+

src/boundary.c

@@ -138,6 +138,50 @@ void reb_boundary_check(struct reb_simulation* const r){
 				}
 			}
 		break;
+		case REB_BOUNDARY_SHEAR_E:
+		{
+			// The offset of ghostcell is time dependent.
+			const double OMEGA = r->ri_sei.OMEGA;
+			const double q = r->ri_sei.Q_NL; // Nonlinearity parameter, 0 < q < 1
+			const double Lx_t = boxsize.x*(1-q*cos(OMEGA*r->t)); // Time dependent box size
+			const double offsetp1 = -fmod(-1.5*OMEGA*boxsize.x*r->t+boxsize.y/2.+2.0*q*boxsize.x*sin(OMEGA*r->t),boxsize.y)-boxsize.y/2.; 
+    		const double offsetm1 = -fmod(1.5*OMEGA*boxsize.x*r->t-boxsize.y/2.-2.0*q*boxsize.x*sin(OMEGA*r->t),boxsize.y)+boxsize.y/2.; 
+			struct reb_particle* const particles = r->particles;
+
+#pragma omp parallel for schedule(guided)
+			for (int i=0;i<N;i++){
+				// Radial
+				while (particles[i].x > Lx_t / 2.) {
+					particles[i].x -= Lx_t;
+					particles[i].y += offsetp1;
+					// Apply epicyclic adjustments to velocity
+					particles[i].vx -= q*boxsize.x*OMEGA*sin(OMEGA*r->t);
+					particles[i].vy += (1.5*boxsize.x*OMEGA-2.0*q*boxsize.x*OMEGA*cos(OMEGA*r->t));
+				}
+				while (particles[i].x < -Lx_t / 2.) {
+					particles[i].x += Lx_t;
+					particles[i].y += offsetm1;
+					// Apply epicyclic adjustments to velocity
+					particles[i].vx += q*boxsize.x*OMEGA*sin(OMEGA*r->t);
+					particles[i].vy += -(1.5*boxsize.x*OMEGA-2.0*q*boxsize.x*OMEGA*cos(OMEGA*r->t));
+				}
+				// Azimuthal
+				while(particles[i].y>boxsize.y/2.){
+					particles[i].y -= boxsize.y;
+				}
+				while(particles[i].y<-boxsize.y/2.){
+					particles[i].y += boxsize.y;
+				}
+				// Vertical (there should be no boundary, but periodic makes life easier)
+				while(particles[i].z>boxsize.z/2.){
+					particles[i].z -= boxsize.z;
+				}
+				while(particles[i].z<-boxsize.z/2.){
+					particles[i].z += boxsize.z;
+				}
+			}
+		}
+		break;
 		default:
 		break;
 	}
@@ -193,6 +237,35 @@ struct reb_vec6d reb_boundary_get_ghostbox(struct reb_simulation* const r, int i
 			gb.vz = 0;
 			return gb;
 		}
+		case REB_BOUNDARY_SHEAR_E:
+		{
+			const double OMEGA = r->ri_sei.OMEGA;
+			const double q = r->ri_sei.Q_NL; // Nonlinearity parameter, 0 < q < 1
+			const double Lx_t = r->boxsize.x*(1-q*cos(OMEGA*r->t)); // Time dependent box size
+
+			struct reb_vec6d gb;
+
+			// Ghostboxes habe a finite velocity.
+			gb.vx = 0.;
+			gb.vy = -1.5*(double)i*OMEGA*r->boxsize.x + 2.0*q*i*r->boxsize.x*OMEGA*cos(OMEGA * r->t); 
+			gb.vz = 0.;
+
+			// The shift in the y direction is time dependent. 
+			double shift;
+			if (i==0){
+				shift = -fmod(gb.vy*r->t,r->boxsize.y); 
+			}else{
+				if (i>0){
+					shift = -fmod(-1.5*(double)i*OMEGA*r->boxsize.x*r->t + 2*q*i*r->boxsize.x*sin(OMEGA*r->t),r->boxsize.y); 
+				}else{
+					shift = -fmod(-1.5*(double)i*OMEGA*r->boxsize.x*r->t + 2*q*i*r->boxsize.x*sin(OMEGA*r->t),r->boxsize.y); 
+				}	
+			}
+			gb.x = Lx_t*(double)i;
+			gb.y = r->boxsize.y*(double)j-shift;
+			gb.z = r->boxsize.z*(double)k;
+			return gb;
+		}
 		default:
 			return nan_ghostbox;
 	}
@@ -222,6 +295,20 @@ int reb_boundary_particle_is_in_box(const struct reb_simulation* const r, struct
 				return 0;
 			}
 			return 1;
+		case REB_BOUNDARY_SHEAR_E:
+		{
+			const double OMEGA = r->ri_sei.OMEGA;
+			const double q = r->ri_sei.Q_NL;  // Nonlinearity parameter
+			double Lx_t = r->boxsize.x*(1-q*cos(OMEGA*r->t)); // Time dependent box size
+
+			// Check boundaries with time dependent changing Lx_t
+			if (p.x > Lx_t / 2. || p.x < -Lx_t / 2.) return 0;
+			if (p.y > r->boxsize.y / 2.) return 0;
+			if (p.y < -r->boxsize.y / 2.) return 0;
+			if (p.z > r->boxsize.z / 2.) return 0;
+			if (p.z < -r->boxsize.z / 2.) return 0;
+			return 1;
+		}
         default:
 		case REB_BOUNDARY_NONE:
 			return 1;

src/output.c

@@ -96,6 +96,7 @@ const struct reb_binary_field_descriptor reb_binary_field_descriptor_list[]= {
     { 53, REB_INT,          "gravity",                      offsetof(struct reb_simulation, gravity), 0, 0},
     { 54, REB_DOUBLE,       "ri_sei.OMEGA",                 offsetof(struct reb_simulation, ri_sei.OMEGA), 0, 0},
     { 55, REB_DOUBLE,       "ri_sei.OMEGAZ",                offsetof(struct reb_simulation, ri_sei.OMEGAZ), 0, 0},
+    { 170, REB_DOUBLE,       "ri_sei.Q_NL",                offsetof(struct reb_simulation, ri_sei.Q_NL), 0, 0},
     { 56, REB_DOUBLE,       "ri_sei.lastdt",                offsetof(struct reb_simulation, ri_sei.lastdt), 0, 0},
     { 57, REB_DOUBLE,       "ri_sei.sindt",                 offsetof(struct reb_simulation, ri_sei.sindt), 0, 0},
     { 58, REB_DOUBLE,       "ri_sei.tandt",                 offsetof(struct reb_simulation, ri_sei.tandt), 0, 0},

src/rebound.h

@@ -248,6 +248,7 @@ struct reb_integrator_mercurius {
 struct reb_integrator_sei {
     double OMEGA;       // Epicyclic frequency
     double OMEGAZ;      // Epicyclic frequency in z direction (if not set, use OMEGA)
+    double Q_NL;        // Nonlinearity parameter
 
     // Internal use
     double lastdt;      // Cached sin(), tan() for this value of dt.
@@ -628,6 +629,7 @@ struct reb_simulation {
         REB_BOUNDARY_OPEN = 1,          // Open boundary conditions. Removes particles if they leave the box 
         REB_BOUNDARY_PERIODIC = 2,      // Periodic boundary conditions
         REB_BOUNDARY_SHEAR = 3,         // Shear periodic boundary conditions, needs OMEGA variable
+        REB_BOUNDARY_SHEAR_E = 4,       // Shear periodic eccentric boundary conditions, needs OMEGA variable
         } boundary;
     enum {
         REB_GRAVITY_NONE = 0,           // Do not calculate graviational forces