Trying to put constraints on particle movement, percolation of water through soil.

Dev Shed Forums Sponsor:

December 10th, 2012, 02:01 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

Trying to put constraints on particle movement, percolation of water through soil.

Hello,
I want to make it so my particles do not follow the random walk unless two or more occupy the same position on the 3d array, I am not sure what to do.

Code:

#percolation.py
#percolation algorithm
from __future__ import division
from visual import*
from visual.graph import*
from random import*


npoint = 50
prob = .001
point = resize(array([0]), (npoint,npoint,npoint)) #making a 3d array.
perc = display(background=color.white, foreground=color.black, width=600, height=600, depth=600, xmin=0, xmax=npoint, ymin=0, ymax=npoint, zmin=0, zmax=npoint, title="Percolation", center=(npoint/2,npoint/2,npoint/2))
step_length = .5

#set up a random array
for i in range(1,(npoint)-1):
    for j in range(1,npoint-1):
        for k in range(1,npoint-1):
            if random()>prob:
                point[i,j,k] = 1
            if random()<prob:
                point[i,j,k] = 0
                box(pos=(i,j,k)) #puts boxes where there is debris like rocks and sand.
particle_list = []

def particle_release():
    for i in range(10):#10 is multiple of particle generation
        particle = sphere(color=color.blue, opacity=.1)
        particle.pos = (npoint/2,npoint-3,npoint/2)#change to origin
        particle_list.append(particle)
particle_release()
          
while True:
    for particle in particle_list:        
        rate(10)
        rand = random()
        if rand <= 1/6: 
            particle.pos.x = particle.pos.x + step_length
        elif rand <= 1/3:
            particle.pos.x = particle.pos.x - step_length
        elif rand <= 2/3:
            particle.pos.y = particle.pos.y - step_length
        elif rand <= 5/6:
            particle.pos.z = particle.pos.z + step_length
        elif rand <= 1:
            particle.pos.z = particle.pos.z - step_length

    particle_release()

December 10th, 2012, 02:31 PM

rrashkin

Contributing User

Join Date: May 2012

Location: 39N 104.28W

Posts: 90

Time spent in forums: 1 Day 13 h 41 m 29 sec
Reputation Power: 2

Quote:

I am not sure what to do.

First of all, use code tags.
Second of all, you may think your code reads like a narrative but at least for me it doesn't. Where is the code behaving incorrectly (use print statements to see what values are being assigned)? What would the correct behavior be?

Comments on this post

b49P23TIvg agrees!

December 10th, 2012, 07:17 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

Sorry, first time I have asked for python help on a forum. Well the program runs fine with respect to generating a random array, generating particles, and random walk. However, I want the random walk to only act on particles if two or more particles occupy an i,j,k coordinate. This way particles will always be in contact, currently they will run off because they do not have code that binds them.
Here is the random walk code that I want to use on particles, but I need to add some code so it only runs if 2 or more occupy the same location.

Code:

while True:
    for particle in particle_list:        
        rate(10)
        rand = random()
        if rand <= 1/6: 
            particle.pos.x = particle.pos.x + step_length
        elif rand <= 1/3:
            particle.pos.x = particle.pos.x - step_length
        elif rand <= 2/3:
            particle.pos.y = particle.pos.y - step_length
        elif rand <= 5/6:
            particle.pos.z = particle.pos.z + step_length
        elif rand <= 1:
            particle.pos.z = particle.pos.z - step_length

December 10th, 2012, 07:30 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

Code:

while True:
    for particle in particle_list:        
        if any other particle also has my coordinates:
	    one way
	else:
            rate(10)
            rand = random()
            if rand <= 1/6: 
                particle.pos.x += step_length
            elif rand <= 1/3:
                particle.pos.x -= step_length

	    # I'm sure you intented to skip 1/2.  Accounts for gravity maybe?

            elif rand <= 2/3:
                particle.pos.y = particle.pos.y - step_length
            elif rand <= 5/6:
                particle.pos.z = particle.pos.z + step_length
            else:
                particle.pos.z = particle.pos.z - step_length

__________________
[code]Code tags[/code] are essential for python code!

Last edited by b49P23TIvg : December 10th, 2012 at 07:32 PM.

December 10th, 2012, 07:45 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

Yeah I skipped it for gravity, I will probably tweak the probabilities to make them more realistic, but first I need to make particles only move if another is in the same position.

Quote:

Originally Posted by b49P23TIvg

Code:

while True:
    for particle in particle_list:        
        if any other particle also has my coordinates:
	    one way
	else:
            rate(10)
            rand = random()
            if rand <= 1/6: 
                particle.pos.x += step_length
            elif rand <= 1/3:
                particle.pos.x -= step_length

	    # I'm sure you intented to skip 1/2.  Accounts for gravity maybe?

            elif rand <= 2/3:
                particle.pos.y = particle.pos.y - step_length
            elif rand <= 5/6:
                particle.pos.z = particle.pos.z + step_length
            else:
                particle.pos.z = particle.pos.z - step_length

December 10th, 2012, 07:59 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

I hope my post wasn't too subtle. Suggested changes here:

Code:

        if any other particle also has my coordinates:
	    one way
	else:
            rate(10)
            rand = random()
            if rand <= 1/6: 
                particle.pos.x += step_length

December 10th, 2012, 08:15 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

Quote:

Originally Posted by b49P23TIvg

I hope my post wasn't too subtle. Suggested changes here:

Code:

        if any other particle also has my coordinates:
	    one way
	else:
            rate(10)
            rand = random()
            if rand <= 1/6: 
                particle.pos.x += step_length

My problem is mainly how we put "if any other particle also has my coordinates" in code?

December 10th, 2012, 08:29 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

Code:

    for (i,particle,) in enumerate(particle_list):
        P = particle.pos
        crowded = any(P == p.pos for p in particle_list[i+1:])

#if that fails, you might have to


TOLERANCE = 1e-8

def same_cell(p,q):
    # look up relative and absolute error and fix this test
    return ((abs(p.x-q.x) < TOLERANCE) and
            (abs(p.y-q.y) < TOLERANCE) and
            (abs(p.z-q.z) < TOLERANCE)


         crowded = any(same_cell(P == p.pos) for p in particle_list[i+1:])

Comments on this post

WynnDeezl agrees!

December 11th, 2012, 08:18 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

So I see the second set of code allows for a little wiggle room. So now I want to make only the particles that satisfy crowded move. So, I was trying out the following but particles were still moving by themselves and spreading out. I must not understand what the code is doing well enough.

Code:

TOLERANCE = 1e-8

def same_cell(p,q):
    # look up relative and absolute error and fix this test
    return ((abs(p.x-q.x) < TOLERANCE) and
            (abs(p.y-q.y) < TOLERANCE) and
            (abs(p.z-q.z) < TOLERANCE))
            
while True:
    for particle in particle_list:
        for (i,particle,) in enumerate(particle_list):
            P = particle.pos
            crowded = any(same_cell(P == p.pos) for p in particle_list[i+1:])
            
                    rate(100)
                    rand = random()
                    if rand <= 1/6:
                        particle.pos.x = particle.pos.x + step_length
                    elif rand <= 1/3:
                        particle.pos.x = particle.pos.x - step_length
                    elif rand <= 2/3:
                        particle.pos.y = particle.pos.y - step_length
                    elif rand <= 5/6:
                        particle.pos.z = particle.pos.z + step_length
                    elif rand <= 1:
                        particle.pos.z = particle.pos.z - step_length
                    print particle_list

    particle_release()

#10

December 11th, 2012, 08:31 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

Is there a missing line of code here??

if crowded then:

Code:

            crowded = any(same_cell(P == p.pos) for p in particle_list[i+1:])
            if crowded then:

                    rate(100)

#11

December 11th, 2012, 08:41 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

yep

Yeah, I accidentaly deleted that line before I copied it, I had
if crowded:
So that makes it if crowded is true it will run, but I only want it to run on particles that satisfy crowded.

#12

December 11th, 2012, 09:18 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

Hey, looks interesting!

Actually, I think the logic is bad.
I'd first list all the particles in crowded cells.
Second step move them.
repeat the two steps until there are few (or none) crowded cells. Repeat for a fixed time. The time step could be inversely proportional to the number of crowded cells, since the crowded cell tally is akin to pressure. Also I'd require that a particle cannot move to a cell with more particles in it. Maybe I'll use this project to test j sparse arrays. www.jsoftware.com .

The requirement "No crowded cells" could take too long. Water seepage under a damn shouldn't interest a practical engineer except as a spare time hobby.

Last edited by b49P23TIvg : December 11th, 2012 at 09:24 PM.

#13

December 11th, 2012, 09:37 PM

red79phoenix

Registered User

Join Date: Dec 2012

Posts: 7

Time spent in forums: 1 h 42 m 57 sec
Reputation Power: 0

Alright. I am not entirely sure what you mean by listing the particles in crowded cells.

#14

December 12th, 2012, 10:38 AM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,380

Time spent in forums: 1 Month 2 Weeks 3 Days 13 h 7 m 19 sec
Reputation Power: 383

The algorithm looks much easier to me if you maintained a list of cells. Each cell has a particle count and pointers to the 6 neighboring cells. Or use a 3D array.

That's why I said the j sparse array is ideal for this project.

Your original post had step length 0.5 . Why? Use integers, then scale everything at the end if integers aren't correct. With integers you could directly compare p.pos == q.pos instead using tolerant comparisons per axis.

Last edited by b49P23TIvg : December 12th, 2012 at 10:41 AM.

Trying to put constraints on particle movement, percolation of water through soil.

Developer Shed Advertisers and Affiliates