Small sigma leads to dt=0 crash

[Jintram]

The following code leads to an error.

This is (probably) because s.t >> dt, which makes it impossible to calculate dt. (More precise: dt is calculated to be 0, and if dt = 0 zero to often, the simulator is halted by a check; see error message below.)

My hypothesis (confirmed by simulation) however is that this is caused by a small sigma. The limiting case is where you have two particles, one with a big domain, and one with a really small domain approaching the big domain. (This doesn't lead to a burst.) What then happens is that the time step of the small domain gets really small, and if some time has already passed in the simulator this leads to a crash [because s.t>>dt].

# Modules
# ===============================
import sys
#Also, set relative egfrd directory path
sys.path.append('/storage1/wehrens/myfork_egfrd')
import os
import shutil
import datetime
import random

import myrandom
import math
import egfrd
import model
import gfrdbase
import _gfrd

# Parameters
# ===============================
matrix_size=1
sigma=1e-10
D=5e-11
k1= (4*math.pi*sigma*(D*2))/100 # Diffusion limit
k2=1e2
world_size = 1e-3

allowedtime = 10

# Functions
# ===============================
def timeseconds():
    """ Returns current time in seconds """
    return (long(datetime.datetime.now().year*3600*24*365+
    datetime.datetime.now().month*30.5*24*3600+
    datetime.datetime.now().day*24*3600+datetime.datetime.now().hour*3600+
    datetime.datetime.now().minute*60+datetime.datetime.now().second))

def put_in_particles(N=1):
    """ Place 1 particle in center, rest in plane through center
    """
    global A, B, w
    gfrdbase.throw_in_particles(w, A, N)
    gfrdbase.throw_in_particles(w, B, N)

# Create "unique" seed
# ===============================
stop_clock = timeseconds()
myrandom.seed(stop_clock)
random.seed(stop_clock)

# Basic set up simulator
# ===============================
# Model
m = model.ParticleModel(world_size)
# Species
A = model.Species('A', D, sigma/2)                                   
m.add_species_type(A) 
B = model.Species('B', D, sigma/2)                                   
m.add_species_type(B) 
C = model.Species('C', D, sigma/2)                                   
m.add_species_type(C) 
# Reaction rules
r1 = model.create_binding_reaction_rule(A, B, C, k1)
m.network_rules.add_reaction_rule(r1)
r2 = model.create_unbinding_reaction_rule(C, A, B, k2)
m.network_rules.add_reaction_rule(r2)

# World
# Default matrix_size = 3
w = gfrdbase.create_world(m, matrix_size) 
# Simulator   
s = egfrd.EGFRDSimulator(w, myrandom.rng)

# Put in particles
put_in_particles()


# Running 
# ===============================

for bla in range(10000):
    s.step()

Dividing up the simulation in smaller runs resolves the issue, supporting the s.t>>dt hypothesis.

for cycle in range (100):
    s.reset()
    for bla in range(100):
        s.step() 

The error message:

Traceback (most recent call last):
  File "simpletest_errorregime.py", line 91, in <module>
    s.step()
  File "/storage1/wehrens/myfork_egfrd/egfrd.py", line 313, in step
    raise RuntimeError('too many dt=zero steps. '
RuntimeError: too many dt=zero steps. Simulator halted?

Note that in this regime the critical time before the simulation crashes is in the order of hours.

[Jintram]

I've adapted some lines in the code to avoid the problem longer:

    # assert if not too many successive dt=0 steps occur.
    if __debug__:
        if self.dt == 0:
            log.warning('dt=zero step, working in s.t >> dt~0 Python limit.')
            self.zero_steps += 1
            # TODO Changed from 10 to 10000, because only a problem 
            # when reaching certain magnitude.
            if self.zero_steps >= max(self.scheduler.size * 3, 10000): 
                raise RuntimeError('too many dt=zero steps. '
                                   'Simulator halted?'
                                'dt= %.300g-%.300g' % (self.scheduler.top[1].time, self.t))

(egfrd.py, ~line 314)

[Jintram]

Note that at the moment also initialisation of domains will give this "warning".

¡We should fix this!