bubble sim

# TO DO #-Figure out plausible units for initial values and constants #-Determine correct coefficient for arbitrary order RK (i.e. RK4 -> 1/6, RK5 -> ?) #-Try to find a way to make RK1 = Euler method and RK2 = Midpoint method #-Dynamic number of steps (# of oscillations, etc.) #-International System units! @interact def solver(method_choice = selector(['Euler', 'Midpoint', 'Runge-Kutta'], default='Runge-Kutta', nrows=1, label='Method'), step_size = slider(vmin=0.0000001, vmax=0.1, step_size=0.0000001, default=1, label='Step size'), rk_order = input_box(default="4", label='Order of Runge-Kutta'), init_radius = input_box(default="1", label='Initial radius'), init_rad_vel = input_box(default="0", label='Initial radial velocity'), init_rad_acc = input_box(default="0", label='Initial radial acceleration'), f_bubble_pressure = input_box(default="1", label='Initial gas pressure inside bubble', width=10), f_v_pressure = input_box(default="1", label='Vapor pressure (p_v(T_inf))', width=10), f_i_pressure = input_box(default="1", label='Fluid pressure far away from cavitation', width=10), f_density = input_box(default="1", label='Liquid density', width=10), #f_var indicates fluid variable f_viscosity = input_box(default="1", label='Liquid viscosity', width=10), f_surface_tension = input_box(default="1", label='Bubble surface tension', width=10), f_constant_k = input_box(default="1", label='Constant K', width=10)): h = step_size rk_order = int(rk_order) R_0 = int(init_radius) RV_0 = int(init_rad_vel) RA_0 = int(init_rad_acc) PB = int(f_bubble_pressure) PV = int(f_v_pressure) PI = int(f_i_pressure) D = int(f_density) V = int(f_viscosity) S = int(f_surface_tension) K = int(f_constant_k) #Function for linearized Rayleigh-Plesset def equation(r, rv): return sin(r)#return rv*((3*RV_0/r) + (4*V/r^2)) + (R_0^K)/(r^(3*K+1)) + (2*S/(3*r^3)) + ((PB-PV-PI)/D) - (3*RV_0^2)/(2*r) #Create some arrays to store values R_array = [[0,R_0]] RV_array = [[0,RV_0]] RA_array = [[0,equation(R_0,RV_0)]] slope_array = [] extrema_array = [] trend = True #True = ascending, False = descending osc_goal = 10 osc_num = 0 #and variables to use in computations r = R_0 rv = RV_0 ra = RA_0 i = 1 #BEGIN MAIN LOOP while i < 1000 or osc_goal < osc_num: if method_choice == 'Euler': #calculate the radius based on previous radius and velocity r = r + rv*h R_array.append([i*h,r]) #calculate the velocity based on acceleration 'ra' and previous velocity rv = rv + ra*h RV_array.append([i*h,rv]) #find acceleration based on linearized equation function ra = equation(r,rv) RA_array.append([i*h,ra]) elif method_choice == 'Midpoint': r = r + rv*h rv = rv + ra*h ra = equation(r,rv) ra = equation(r+h/2, rv+ra/2) #append to arrays RA_array.append([i*h,ra]) RV_array.append([i*h,rv]) R_array.append([i*h,r]) elif method_choice == 'Runge-Kutta': r = r + rv*h R_array.append([i*h,r]) rv = rv + ra*h RV_array.append([i*h,rv]) #calculation of 'ra' for arbitrary RK order with delta values stored in 'k_array' k_array = [] k_array.append(h*equation(r,rv)) for a in range(rk_order-2): k_array.append(2*h*equation(r+h/2, rv+k_array[a]/2)) k_array.append(h*equation(r+h,rv+k_array[rk_order-2])) ra = ra + sum(k_array)/6 RA_array.append([i*h,ra]) #calculate slope and detect extrema slope = (R_array[i][1] - R_array[i-1][1])/h slope_array.append(slope.n()) if len(extrema_array) == 0: trend = R_array[1][1] > R_array[0][1] #on first iteration, set 'trend' elif R_array[i][1] > R_array[i-1][1] and trend == False: #detect the new slope; if it changes, declare an extrema and switch 'trend' trend = True extrema_array.append(R_array[i-1]) elif R_array[i][1] < R_array[i-1][1] and trend == True: #same as above but for the opposite direction trend = False extrema_array.append(R_array[i-1]) #check for low amounts of activity (overflow safeguard #1) avg_change = 0 tot_change = 0 #not really the 'total' change - takes sign changes into account, unlike avg_change l = len(slope_array) if l > 5: for a in range(5): avg_change += abs(slope_array[l-a-1] - slope_array[l-a-2]) tot_change += slope_array[l-a-1] - slope_array[l-a-2] avg_change /= 5 tot_change /= 5 if avg_change < 0.000005: print "Computation interrupted due to lack of activity\n\t(Average change = %s)" % avg_change break if avg_change > 100 and avg_change == abs(tot_change): print "Computation interrupted due to exponential behavior" break i += 1 #END MAIN LOOP #Create a line that connects all of the points in r_array L = line([(R_array[j][0], R_array[j][1]) for j in range(len(R_array))], color=(1,0,0), legend_label="Radial position") show(L) #Same thing, RV_array M = line([(RV_array[j][0], RV_array[j][1]) for j in range(len(RV_array))], color=(0,1,0), legend_label="Radial velocity") show(M) #ctrl-c ctrl-v RA_array N = line([(RA_array[j][0], RA_array[j][1]) for j in range(len(RA_array))], color=(0,0,1), legend_label="Radial accel") show(N) show(L+M+N) 

#Implement this code into the main bubble simulation #-Convert main simulations to 'while' loops #-Add in the overflow/nonactivity detector #-Create function to evaluate an extreme point given a slice of the array #-HAVE THE RANGES FOR OVERFLOW DETECTION SCALE WITH STEP SIZE #Finish this code #-Account for a triple peak with middle peak being lowest #-Implement run-time extrema evaluation (see below) array = [[0,0]] slope_array = [] extrema_array = [] n = 0.1 #step size trend = True #True = ascending, False = descending osc_goal = 10 osc_num = 0 i = 1 while osc_num < osc_goal: #calculate new x, y, and slope; add them to respective arrays x = i*n y = sin(x) + 0.1*sin(10*x)# + 0.2*sin(20*x)# + 0.3*sin(30*x) array.append([x,y]) slope = (array[i][1] - array[i-1][1])/(4*pi/n) slope_array.append(slope.n()) #detect extrema if i == 0: trend = y > array[0][1] #on first iteration, set 'trend' elif array[i][1] > array[i-1][1] and trend == False: #detect the new slope; if it changes, declare an extrema and switch 'trend' trend = True extrema_array.append(array[i-1]) elif array[i][1] < array[i-1][1] and trend == True: #same as above but for the opposite direction trend = False extrema_array.append(array[i-1]) #guess at a typical length of oscillation #possible things to look for: distance traveled in y direction sorted(extrema_array) #check for low amounts of activity or exponential behavior avg_change = 0 tot_change = 0 #not really the 'total' change - takes sign changes into account, unlike avg_change l = len(slope_array) if l > 5: for a in range(5): avg_change += abs(slope_array[l-a-1] - slope_array[l-a-2]) tot_change += slope_array[l-a-1] - slope_array[l-a-2] avg_change /= 5 tot_change /= 5 if avg_change < 0.000005: print "Computation interrupted due to lack of activity\n\t(Average change = %s)" % avg_change break if avg_change > 100 and avg_change == abs(tot_change): print "Computation interrupted due to exponential behavior" break i += 1 #Convert this into a run-time evaluation of extrema #REMEMBER: You can't have 2 consecutive maxima - must be a minimum in between! #Think about comparing y-distance of extrema to known min and max extrema_array_2 = [] #this will hold all the "significant" extrema prev_value = "" for a in range(len(extrema_array)-2): if a < 2: continue search_point = list(extrema_array[a]) prev_point = list(extrema_array[a-1]) prev_ext = list(extrema_array[a-2]) next_ext = list(extrema_array[a+2]) if (search_point[1] > prev_point[1] and search_point[1] >= prev_ext[1] and search_point[1] >= next_ext[1]): search_point.append("max") if (len(extrema_array_2) != 0): if (extrema_array_2[len(extrema_array_2)-1][2] == "min"): extrema_array_2.append(search_point) elif (search_point[1].n() > extrema_array_2[len(extrema_array_2)-1][1].n()): extrema_array_2.pop() extrema_array_2.append(search_point) elif (search_point[1].n() < extrema_array_2[len(extrema_array_2)-1][1].n()): pass else: extrema_array_2.append(search_point) elif (search_point[1] < prev_point[1] and search_point[1] <= prev_ext[1] and search_point[1] <= next_ext[1]): search_point.append("min") if (len(extrema_array_2) > 0): if (extrema_array_2[len(extrema_array_2)-1][2] == "max"): extrema_array_2.append(search_point) elif (search_point[1].n() < extrema_array_2[len(extrema_array_2)-1][1].n()): extrema_array_2.pop() extrema_array_2.append(search_point) elif (search_point[1].n() > extrema_array_2[len(extrema_array_2)-1][1].n()): pass else: extrema_array_2.append(search_point) osc_num = len(extrema_array_2)/2 print len(extrema_array_2) graph = line([(array[j][0].n(), array[j][1].n()) for j in range(len(array))]) points = point([(extrema_array_2[j][0].n(), extrema_array_2[j][1].n()) for j in range(len(extrema_array_2))], rgbcolor=(1,0,0), pointsize=18) points2 = point([(extrema_array[j][0].n(), extrema_array[j][1].n()) for j in range(len(extrema_array))], rgbcolor=(0,1,0), pointsize=18) show(points2 + points + graph) 

20

20

#DETERMINING 'SIGNIFICANT' EXTREMA #-An extrema is significant if the two surrounding extrema of the same type (minimum/maximum) are both higher or both lower, respectively #-An extrema is not significant if its neighboring extrema is of the same type and of greater magnitude #-An extrema is *probably* not significant if it the change in Y between it and the neighboring extrema of the opposite type is significantly less than that of the furthest separated pair # -Thought: Margin of acceptability defined by: (distance between absolute min and max) - (distance between furthest apart extrema) 

def runge_kutta_two(a,b,N,alpha): #create functions we will need f = lambda x,y: 1.0 + (y/x) + (y/x)^2 #this is the first equation g = lambda x: x*tan(ln(x)) # this is the exact solution for the IVP #continuing with page 279 h = (b-a)/N t = a omega = alpha #make a list to store values w = [] w.append((0,omega, numerical_approx(abs(g(1)-omega)))) # this stores y(1) = 0 & actual error for i in range(1,N): omega = omega + h * f(t + (h/2), omega + (h/t)*f(t,omega)) # compute omega_i w.append( ( i, omega, numerical_approx(abs(g(1)-omega))) ) #store values t = a + i * h # compute t_i return w print runge_kutta_two(1,100,100,0) #a = start point #b = endpoint #N = number of steps #alpha = #h = step size 

if method_choice == 'Euler': for i in range(math.ceil(n/h)): #calculate the radius based on previous radius and velocity r = r + rv*h R_array.append([i*h,r]) #calculate the velocity based on acceleration 'ra' and previous velocity rv = rv + ra*h RV_array.append([i*h,rv]) #find acceleration based on linearized equation function ra = equation(r,rv) RA_array.append([i*h,ra]) elif method_choice == 'Midpoint': for i in range(math.ceil(n/h)): r = r + rv*h rv = rv + ra*h ra = equation(r,rv) ra = equation(r+h/2, rv+ra/2) #append to arrays RA_array.append([i*h,ra]) RV_array.append([i*h,rv]) R_array.append([i*h,r]) elif method_choice == 'Runge-Kutta': for i in range(math.ceil(n/h)): r = r + rv*h R_array.append([i*h,r]) rv = rv + ra*h RV_array.append([i*h,rv]) #calculation of 'ra' for arbitrary RK order with delta values stored in 'k_array' k_array = [] k_array.append(h*equation(r,rv)) for a in range(rk_order-2): k_array.append(2*h*equation(r+h/2, rv+k_array[a]/2)) k_array.append(h*equation(r+h,rv+k_array[rk_order-2])) ra = ra + sum(k_array)/6 RA_array.append([i*h,ra]) 

Traceback (click to the left of this block for traceback)
...
IndentationError: unindent does not match any outer indentation level

Traceback (most recent call last):                    rv = rv + ra*h
  File "", line 1, in <module>
    
  File "/tmp/tmpL_gBWK/___code___.py", line 14
    elif method_choice == 'Midpoint':
                                    ^
IndentationError: unindent does not match any outer indentation level