#NEWLINE# def evolve( self ) : #NEWLINE# #INDENT# '\n        Evolve the incompressible equations through one timestep.\n        ' #NEWLINE# u = self . cc_data . get_var ( 'x-velocity' ) #NEWLINE# v = self . cc_data . get_var ( 'y-velocity' ) #NEWLINE# gradp_x = self . cc_data . get_var ( 'gradp_x' ) #NEWLINE# gradp_y = self . cc_data . get_var ( 'gradp_y' ) #NEWLINE# phi = self . cc_data . get_var ( 'phi' ) #NEWLINE# myg = self . cc_data . grid #NEWLINE# limiter = self . rp . get_param ( 'incompressible.limiter' ) #NEWLINE# if ( limiter == 0 ) : #NEWLINE# #INDENT# limitFunc = reconstruction_f . nolimit #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( limiter == 1 ) : #NEWLINE# #INDENT# limitFunc = reconstruction_f . limit2 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# limitFunc = reconstruction_f . limit4 #NEWLINE# #UNINDENT# #UNINDENT# ldelta_ux = limitFunc ( 1 , u . d , myg . qx , myg . qy , myg . ng ) #NEWLINE# ldelta_vx = limitFunc ( 1 , v . d , myg . qx , myg . qy , myg . ng ) #NEWLINE# ldelta_uy = limitFunc ( 2 , u . d , myg . qx , myg . qy , myg . ng ) #NEWLINE# ldelta_vy = limitFunc ( 2 , v . d , myg . qx , myg . qy , myg . ng ) #NEWLINE# '\n        the advective velocities are the normal velocity through each cell\n        interface, and are defined on the cell edges, in a MAC type\n        staggered form\n\n                         n+1/2\n                        v\n                         i,j+1/2\n                    +------+------+\n                    |             |\n            n+1/2   |             |   n+1/2\n           u        +     U       +  u\n            i-1/2,j |      i,j    |   i+1/2,j\n                    |             |\n                    +------+------+\n                         n+1/2\n                        v\n                         i,j-1/2\n\n        ' #NEWLINE# if ( self . verbose > 0 ) : #NEWLINE# #INDENT# print ( '  making MAC velocities' ) #NEWLINE# #UNINDENT# ( _um , _vm ) = incomp_interface_f . mac_vels ( myg . qx , myg . qy , myg . ng , myg . dx , myg . dy , self . dt , u . d , v . d , ldelta_ux , ldelta_vx , ldelta_uy , ldelta_vy , gradp_x . d , gradp_y . d ) #NEWLINE# u_MAC = patch . ArrayIndexer ( d = _um , grid = myg ) #NEWLINE# v_MAC = patch . ArrayIndexer ( d = _vm , grid = myg ) #NEWLINE# if ( self . verbose > 0 ) : #NEWLINE# #INDENT# print ( '  MAC projection' ) #NEWLINE# #UNINDENT# mg = MG . CellCenterMG2d ( myg . nx , myg . ny , xl_BC_type = 'periodic' , xr_BC_type = 'periodic' , yl_BC_type = 'periodic' , yr_BC_type = 'periodic' , xmin = myg . xmin , xmax = myg . xmax , ymin = myg . ymin , ymax = myg . ymax , verbose = 0 ) #NEWLINE# divU = mg . soln_grid . scratch_array ( ) #NEWLINE# divU . v ( ) [ : , : ] = ( ( ( u_MAC . ip ( 1 ) - u_MAC . v ( ) ) / myg . dx ) + ( ( v_MAC . jp ( 1 ) - v_MAC . v ( ) ) / myg . dy ) ) #NEWLINE# mg . init_zeros ( ) #NEWLINE# mg . init_RHS ( divU . d ) #NEWLINE# mg . solve ( rtol = 1e-12 ) #NEWLINE# phi_MAC = self . cc_data . get_var ( 'phi-MAC' ) #NEWLINE# solution = mg . get_solution ( ) #NEWLINE# phi_MAC . v ( buf = 1 ) [ : , : ] = solution . v ( buf = 1 ) #NEWLINE# b = ( 0 , 1 , 0 , 0 ) #NEWLINE# u_MAC . v ( buf = b ) [ : , : ] -= ( ( phi_MAC . v ( buf = b ) - phi_MAC . ip ( ( - 1 ) , buf = b ) ) / myg . dx ) #NEWLINE# b = ( 0 , 0 , 0 , 1 ) #NEWLINE# v_MAC . v ( buf = b ) [ : , : ] -= ( ( phi_MAC . v ( buf = b ) - phi_MAC . jp ( ( - 1 ) , buf = b ) ) / myg . dy ) #NEWLINE# if ( self . verbose > 0 ) : #NEWLINE# #INDENT# print ( '  making u, v edge states' ) #NEWLINE# #UNINDENT# ( _ux , _vx , _uy , _vy ) = incomp_interface_f . states ( myg . qx , myg . qy , myg . ng , myg . dx , myg . dy , self . dt , u . d , v . d , ldelta_ux , ldelta_vx , ldelta_uy , ldelta_vy , gradp_x . d , gradp_y . d , u_MAC . d , v_MAC . d ) #NEWLINE# u_xint = patch . ArrayIndexer ( d = _ux , grid = myg ) #NEWLINE# v_xint = patch . ArrayIndexer ( d = _vx , grid = myg ) #NEWLINE# u_yint = patch . ArrayIndexer ( d = _uy , grid = myg ) #NEWLINE# v_yint = patch . ArrayIndexer ( d = _vy , grid = myg ) #NEWLINE# if ( self . verbose > 0 ) : #NEWLINE# #INDENT# print ( '  doing provisional update of u, v' ) #NEWLINE# #UNINDENT# advect_x = myg . scratch_array ( ) #NEWLINE# advect_y = myg . scratch_array ( ) #NEWLINE# advect_x . v ( ) [ : , : ] = ( ( ( ( 0.5 * ( u_MAC . v ( ) + u_MAC . ip ( 1 ) ) ) * ( u_xint . ip ( 1 ) - u_xint . v ( ) ) ) / myg . dx ) + ( ( ( 0.5 * ( v_MAC . v ( ) + v_MAC . jp ( 1 ) ) ) * ( u_yint . jp ( 1 ) - u_yint . v ( ) ) ) / myg . dy ) ) #NEWLINE# advect_y . v ( ) [ : , : ] = ( ( ( ( 0.5 * ( u_MAC . v ( ) + u_MAC . ip ( 1 ) ) ) * ( v_xint . ip ( 1 ) - v_xint . v ( ) ) ) / myg . dx ) + ( ( ( 0.5 * ( v_MAC . v ( ) + v_MAC . jp ( 1 ) ) ) * ( v_yint . jp ( 1 ) - v_yint . v ( ) ) ) / myg . dy ) ) #NEWLINE# proj_type = self . rp . get_param ( 'incompressible.proj_type' ) #NEWLINE# if ( proj_type == 1 ) : #NEWLINE# #INDENT# u . d [ : , : ] -= ( ( self . dt * advect_x . d [ : , : ] ) + ( self . dt * gradp_x . d [ : , : ] ) ) #NEWLINE# v . d [ : , : ] -= ( ( self . dt * advect_y . d [ : , : ] ) + ( self . dt * gradp_y . d [ : , : ] ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( proj_type == 2 ) : #NEWLINE# #INDENT# u . d [ : , : ] -= ( self . dt * advect_x . d [ : , : ] ) #NEWLINE# v . d [ : , : ] -= ( self . dt * advect_y . d [ : , : ] ) #NEWLINE# #UNINDENT# #UNINDENT# self . cc_data . fill_BC ( 'x-velocity' ) #NEWLINE# self . cc_data . fill_BC ( 'y-velocity' ) #NEWLINE# if ( self . verbose > 0 ) : #NEWLINE# #INDENT# print ( '  final projection' ) #NEWLINE# #UNINDENT# mg = MG . CellCenterMG2d ( myg . nx , myg . ny , xl_BC_type = 'periodic' , xr_BC_type = 'periodic' , yl_BC_type = 'periodic' , yr_BC_type = 'periodic' , xmin = myg . xmin , xmax = myg . xmax , ymin = myg . ymin , ymax = myg . ymax , verbose = 0 ) #NEWLINE# divU . v ( ) [ : , : ] = ( ( ( 0.5 * ( u . ip ( 1 ) - u . ip ( ( - 1 ) ) ) ) / myg . dx ) + ( ( 0.5 * ( v . jp ( 1 ) - v . jp ( ( - 1 ) ) ) ) / myg . dy ) ) #NEWLINE# mg . init_RHS ( ( divU . d / self . dt ) ) #NEWLINE# phiGuess = mg . soln_grid . scratch_array ( ) #NEWLINE# phiGuess . v ( buf = 1 ) [ : , : ] = phi . v ( buf = 1 ) #NEWLINE# mg . init_solution ( phiGuess . d ) #NEWLINE# mg . solve ( rtol = 1e-12 ) #NEWLINE# phi . d [ : , : ] = mg . get_solution ( grid = myg ) . d #NEWLINE# ( gradphi_x , gradphi_y ) = mg . get_solution_gradient ( grid = myg ) #NEWLINE# u . d [ : , : ] -= ( self . dt * gradphi_x . d ) #NEWLINE# v . d [ : , : ] -= ( self . dt * gradphi_y . d ) #NEWLINE# if ( proj_type == 1 ) : #NEWLINE# #INDENT# gradp_x . d [ : , : ] += gradphi_x . d [ : , : ] #NEWLINE# gradp_y . d [ : , : ] += gradphi_y . d [ : , : ] #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( proj_type == 2 ) : #NEWLINE# #INDENT# gradp_x . d [ : , : ] = gradphi_x . d [ : , : ] #NEWLINE# gradp_y . d [ : , : ] = gradphi_y . d [ : , : ] #NEWLINE# #UNINDENT# #UNINDENT# self . cc_data . fill_BC ( 'x-velocity' ) #NEWLINE# self . cc_data . fill_BC ( 'y-velocity' ) #NEWLINE# if ( not self . in_preevolve ) : #NEWLINE# #INDENT# self . cc_data . t += self . dt #NEWLINE# self . n += 1