<s> from __future__ import print_function <EOL> import sys <EOL> import numpy as np <EOL> import matplotlib . pyplot as plt <EOL> from lm_atm . problems import * <EOL> import lm_atm . LM_atm_interface_f as lm_interface_f <EOL> import mesh . reconstruction_f as reconstruction_f <EOL> import mesh . patch as patch <EOL> from simulation_null import NullSimulation , grid_setup , bc_setup <EOL> import multigrid . variable_coeff_MG as vcMG <EOL> from util import profile <EOL> class Basestate ( object ) : <EOL> def __init__ ( self , ny , ng = <NUM_LIT:0> ) : <EOL> self . ny = ny <EOL> self . ng = ng <EOL> self . qy = ny + <NUM_LIT:2> * ng <EOL> self . d = np . zeros ( ( self . qy ) , dtype = np . float64 ) <EOL> self . jlo = ng <EOL> self . jhi = ng + ny - <NUM_LIT:1> <EOL> def v ( self , buf = <NUM_LIT:0> ) : <EOL> return self . d [ self . jlo - buf : self . jhi + <NUM_LIT:1> + buf ] <EOL> def v2d ( self , buf = <NUM_LIT:0> ) : <EOL> return self . d [ np . newaxis , self . jlo - buf : self . jhi + <NUM_LIT:1> + buf ] <EOL> def v2dp ( self , shift , buf = <NUM_LIT:0> ) : <EOL> return self . d [ np . newaxis , self . jlo + shift - buf : self . jhi + <NUM_LIT:1> + shift + buf ] <EOL> def jp ( self , shift , buf = <NUM_LIT:0> ) : <EOL> return self . d [ self . jlo - buf + shift : self . jhi + <NUM_LIT:1> + buf + shift ] <EOL> class Simulation ( NullSimulation ) : <EOL> def __init__ ( self , solver_name , problem_name , rp , timers = None ) : <EOL> NullSimulation . __init__ ( self , solver_name , problem_name , rp , timers = timers ) <EOL> self . base = { } <EOL> self . aux_data = None <EOL> def initialize ( self ) : <EOL> """<STR_LIT>""" <EOL> myg = grid_setup ( self . rp , ng = <NUM_LIT:4> ) <EOL> bc_dens , bc_xodd , bc_yodd = bc_setup ( self . rp ) <EOL> my_data = patch . CellCenterData2d ( myg ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_dens ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_xodd ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_yodd ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_dens ) <EOL> bcs = [ ] <EOL> for bc in [ self . rp . get_param ( "<STR_LIT>" ) , <EOL> self . rp . get_param ( "<STR_LIT>" ) , <EOL> self . rp . get_param ( "<STR_LIT>" ) , <EOL> self . rp . get_param ( "<STR_LIT>" ) ] : <EOL> if bc == "<STR_LIT>" : <EOL> bctype = "<STR_LIT>" <EOL> elif bc in [ "<STR_LIT>" , "<STR_LIT>" ] : <EOL> bctype = "<STR_LIT>" <EOL> elif bc in [ "<STR_LIT>" ] : <EOL> bctype = "<STR_LIT>" <EOL> bcs . append ( bctype ) <EOL> bc_phi = patch . BCObject ( xlb = bcs [ <NUM_LIT:0> ] , xrb = bcs [ <NUM_LIT:1> ] , ylb = bcs [ <NUM_LIT:2> ] , yrb = bcs [ <NUM_LIT:3> ] ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_phi ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_phi ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_dens ) <EOL> my_data . register_var ( "<STR_LIT>" , bc_dens ) <EOL> my_data . create ( ) <EOL> self . cc_data = my_data <EOL> aux_data = patch . CellCenterData2d ( myg ) <EOL> aux_data . register_var ( "<STR_LIT>" , bc_dens ) <EOL> aux_data . register_var ( "<STR_LIT>" , bc_yodd ) <EOL> aux_data . create ( ) <EOL> self . aux_data = aux_data <EOL> self . base [ "<STR_LIT>" ] = Basestate ( myg . ny , ng = myg . ng ) <EOL> self . base [ "<STR_LIT>" ] = Basestate ( myg . ny , ng = myg . ng ) <EOL> exec ( self . problem_name + '<STR_LIT>' ) <EOL> gamma = self . rp . get_param ( "<STR_LIT>" ) <EOL> self . base [ "<STR_LIT>" ] = Basestate ( myg . ny , ng = myg . ng ) <EOL> self . base [ "<STR_LIT>" ] . d [ : ] = self . base [ "<STR_LIT>" ] . d ** ( <NUM_LIT:1.0> / gamma ) <EOL> self . base [ "<STR_LIT>" ] = Basestate ( myg . ny , ng = myg . ng ) <EOL> self . base [ "<STR_LIT>" ] . jp ( <NUM_LIT:1> ) [ : ] = <NUM_LIT:0.5> * ( self . base [ "<STR_LIT>" ] . v ( ) + self . base [ "<STR_LIT>" ] . jp ( <NUM_LIT:1> ) ) <EOL> self . base [ "<STR_LIT>" ] . d [ myg . jlo ] = self . base [ "<STR_LIT>" ] . d [ myg . jlo ] <EOL> self . base [ "<STR_LIT>" ] . d [ myg . jhi + <NUM_LIT:1> ] = self . base [ "<STR_LIT>" ] . d [ myg . jhi ] <EOL> def make_prime ( self , a , a0 ) : <EOL> return a - a0 . v2d ( buf = a0 . ng ) <EOL> def compute_timestep ( self ) : <EOL> """<STR_LIT>""" <EOL> myg = self . cc_data . grid <EOL> cfl = self . rp . get_param ( "<STR_LIT>" ) <EOL> u = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> v = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> xtmp = ytmp = <NUM_LIT> <EOL> if not abs ( u ) . max ( ) == <NUM_LIT:0> : <EOL> xtmp = myg . dx / abs ( u . v ( ) ) . max ( ) <EOL> if not abs ( v ) . max ( ) == <NUM_LIT:0> : <EOL> ytmp = myg . dy / abs ( v . v ( ) ) . max ( ) <EOL> dt = cfl * min ( xtmp , ytmp ) <EOL> rho = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> rho0 = self . base [ "<STR_LIT>" ] <EOL> rhoprime = self . make_prime ( rho , rho0 ) <EOL> g = self . rp . get_param ( "<STR_LIT>" ) <EOL> F_buoy = ( abs ( rhoprime * g ) . v ( ) / rho . v ( ) ) . max ( ) <EOL> dt_buoy = np . sqrt ( <NUM_LIT> * myg . dx / F_buoy ) <EOL> self . dt = min ( dt , dt_buoy ) <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" . format ( dt ) ) <EOL> def preevolve ( self ) : <EOL> """<STR_LIT>""" <EOL> self . in_preevolve = True <EOL> myg = self . cc_data . grid <EOL> rho = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> u = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> v = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> coeff = <NUM_LIT:1> / rho <EOL> beta0 = self . base [ "<STR_LIT>" ] <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) ** <NUM_LIT:2> <EOL> mg = vcMG . VarCoeffCCMG2d ( myg . nx , myg . ny , <EOL> xl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xlb , <EOL> xr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xrb , <EOL> yl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . ylb , <EOL> yr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . yrb , <EOL> xmin = myg . xmin , xmax = myg . xmax , <EOL> ymin = myg . ymin , ymax = myg . ymax , <EOL> coeffs = coeff , <EOL> coeffs_bc = self . cc_data . BCs [ "<STR_LIT>" ] , <EOL> verbose = <NUM_LIT:0> ) <EOL> div_beta_U = mg . soln_grid . scratch_array ( ) <EOL> div_beta_U . v ( ) [ : , : ] = <NUM_LIT:0.5> * beta0 . v2d ( ) * ( u . ip ( <NUM_LIT:1> ) - u . ip ( - <NUM_LIT:1> ) ) / myg . dx + <NUM_LIT:0.5> * ( beta0 . v2dp ( <NUM_LIT:1> ) * v . jp ( <NUM_LIT:1> ) - beta0 . v2dp ( - <NUM_LIT:1> ) * v . jp ( - <NUM_LIT:1> ) ) / myg . dy <EOL> mg . init_RHS ( div_beta_U . d ) <EOL> mg . solve ( rtol = <NUM_LIT> ) <EOL> phi = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> phi . d [ : , : ] = mg . get_solution ( grid = myg ) . d <EOL> gradp_x , gradp_y = mg . get_solution_gradient ( grid = myg ) <EOL> coeff = <NUM_LIT:1.0> / rho <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) <EOL> u . v ( ) [ : , : ] -= coeff . v ( ) * gradp_x . v ( ) <EOL> v . v ( ) [ : , : ] -= coeff . v ( ) * gradp_y . v ( ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> orig_data = patch . cell_center_data_clone ( self . cc_data ) <EOL> self . compute_timestep ( ) <EOL> self . evolve ( ) <EOL> new_gp_x = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> new_gp_y = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> orig_gp_x = orig_data . get_var ( "<STR_LIT>" ) <EOL> orig_gp_y = orig_data . get_var ( "<STR_LIT>" ) <EOL> orig_gp_x . d [ : , : ] = new_gp_x . d [ : , : ] <EOL> orig_gp_y . d [ : , : ] = new_gp_y . d [ : , : ] <EOL> self . cc_data = orig_data <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> self . in_preevolve = False <EOL> def evolve ( self ) : <EOL> """<STR_LIT>""" <EOL> rho = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> u = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> v = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> gradp_x = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> gradp_y = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> beta0 = self . base [ "<STR_LIT>" ] <EOL> beta0_edges = self . base [ "<STR_LIT>" ] <EOL> rho0 = self . base [ "<STR_LIT>" ] <EOL> phi = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> myg = self . cc_data . grid <EOL> limiter = self . rp . get_param ( "<STR_LIT>" ) <EOL> if limiter == <NUM_LIT:0> : limitFunc = reconstruction_f . nolimit <EOL> elif limiter == <NUM_LIT:1> : limitFunc = reconstruction_f . limit2 <EOL> else : limitFunc = reconstruction_f . limit4 <EOL> ldelta_rx = limitFunc ( <NUM_LIT:1> , rho . d , myg . qx , myg . qy , myg . ng ) <EOL> ldelta_ux = limitFunc ( <NUM_LIT:1> , u . d , myg . qx , myg . qy , myg . ng ) <EOL> ldelta_vx = limitFunc ( <NUM_LIT:1> , v . d , myg . qx , myg . qy , myg . ng ) <EOL> ldelta_ry = limitFunc ( <NUM_LIT:2> , rho . d , myg . qx , myg . qy , myg . ng ) <EOL> ldelta_uy = limitFunc ( <NUM_LIT:2> , u . d , myg . qx , myg . qy , myg . ng ) <EOL> ldelta_vy = limitFunc ( <NUM_LIT:2> , v . d , myg . qx , myg . qy , myg . ng ) <EOL> """<STR_LIT>""" <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> coeff = self . aux_data . get_var ( "<STR_LIT>" ) <EOL> coeff . v ( ) [ : , : ] = <NUM_LIT:1.0> / rho . v ( ) <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) <EOL> self . aux_data . fill_BC ( "<STR_LIT>" ) <EOL> source = self . aux_data . get_var ( "<STR_LIT>" ) <EOL> g = self . rp . get_param ( "<STR_LIT>" ) <EOL> rhoprime = self . make_prime ( rho , rho0 ) <EOL> source . v ( ) [ : , : ] = rhoprime . v ( ) * g / rho . v ( ) <EOL> self . aux_data . fill_BC ( "<STR_LIT>" ) <EOL> _um , _vm = lm_interface_f . mac_vels ( myg . qx , myg . qy , myg . ng , <EOL> myg . dx , myg . dy , self . dt , <EOL> u . d , v . d , <EOL> ldelta_ux , ldelta_vx , <EOL> ldelta_uy , ldelta_vy , <EOL> coeff . d * gradp_x . d , coeff . d * gradp_y . d , <EOL> source . d ) <EOL> u_MAC = patch . ArrayIndexer ( d = _um , grid = myg ) <EOL> v_MAC = patch . ArrayIndexer ( d = _vm , grid = myg ) <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> coeff . v ( buf = <NUM_LIT:1> ) [ : , : ] = <NUM_LIT:1.0> / rho . v ( buf = <NUM_LIT:1> ) <EOL> coeff . v ( buf = <NUM_LIT:1> ) [ : , : ] = coeff . v ( buf = <NUM_LIT:1> ) * beta0 . v2d ( buf = <NUM_LIT:1> ) ** <NUM_LIT:2> <EOL> mg = vcMG . VarCoeffCCMG2d ( myg . nx , myg . ny , <EOL> xl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xlb , <EOL> xr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xrb , <EOL> yl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . ylb , <EOL> yr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . yrb , <EOL> xmin = myg . xmin , xmax = myg . xmax , <EOL> ymin = myg . ymin , ymax = myg . ymax , <EOL> coeffs = coeff , <EOL> coeffs_bc = self . cc_data . BCs [ "<STR_LIT>" ] , <EOL> verbose = <NUM_LIT:0> ) <EOL> div_beta_U = mg . soln_grid . scratch_array ( ) <EOL> div_beta_U . v ( ) [ : , : ] = beta0 . v2d ( ) * ( u_MAC . ip ( <NUM_LIT:1> ) - u_MAC . v ( ) ) / myg . dx + ( beta0_edges . v2dp ( <NUM_LIT:1> ) * v_MAC . jp ( <NUM_LIT:1> ) - <EOL> beta0_edges . v2d ( ) * v_MAC . v ( ) ) / myg . dy <EOL> mg . init_RHS ( div_beta_U . d ) <EOL> mg . solve ( rtol = <NUM_LIT> ) <EOL> phi_MAC = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> phi_MAC . d [ : , : ] = mg . get_solution ( grid = myg ) . d <EOL> coeff = self . aux_data . get_var ( "<STR_LIT>" ) <EOL> coeff . v ( ) [ : , : ] = <NUM_LIT:1.0> / rho . v ( ) <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) <EOL> self . aux_data . fill_BC ( "<STR_LIT>" ) <EOL> coeff_x = myg . scratch_array ( ) <EOL> b = ( <NUM_LIT:3> , <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> coeff_x . v ( buf = b ) [ : , : ] = <NUM_LIT:0.5> * ( coeff . ip ( - <NUM_LIT:1> , buf = b ) + coeff . v ( buf = b ) ) <EOL> coeff_y = myg . scratch_array ( ) <EOL> b = ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:3> , <NUM_LIT:1> ) <EOL> coeff_y . v ( buf = b ) [ : , : ] = <NUM_LIT:0.5> * ( coeff . jp ( - <NUM_LIT:1> , buf = b ) + coeff . v ( buf = b ) ) <EOL> b = ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> u_MAC . v ( buf = b ) [ : , : ] -= coeff_x . v ( buf = b ) * ( phi_MAC . v ( buf = b ) - phi_MAC . ip ( - <NUM_LIT:1> , buf = b ) ) / myg . dx <EOL> b = ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> v_MAC . v ( buf = b ) [ : , : ] -= coeff_y . v ( buf = b ) * ( phi_MAC . v ( buf = b ) - phi_MAC . jp ( - <NUM_LIT:1> , buf = b ) ) / myg . dy <EOL> _rx , _ry = lm_interface_f . rho_states ( myg . qx , myg . qy , myg . ng , <EOL> myg . dx , myg . dy , self . dt , <EOL> rho . d , u_MAC . d , v_MAC . d , <EOL> ldelta_rx , ldelta_ry ) <EOL> rho_xint = patch . ArrayIndexer ( d = _rx , grid = myg ) <EOL> rho_yint = patch . ArrayIndexer ( d = _ry , grid = myg ) <EOL> rho_old = rho . copy ( ) <EOL> rho . v ( ) [ : , : ] -= self . dt * ( <EOL> ( rho_xint . ip ( <NUM_LIT:1> ) * u_MAC . ip ( <NUM_LIT:1> ) - rho_xint . v ( ) * u_MAC . v ( ) ) / myg . dx + <EOL> ( rho_yint . jp ( <NUM_LIT:1> ) * v_MAC . jp ( <NUM_LIT:1> ) - rho_yint . v ( ) * v_MAC . v ( ) ) / myg . dy ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> eint = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> gamma = self . rp . get_param ( "<STR_LIT>" ) <EOL> eint . v ( ) [ : , : ] = self . base [ "<STR_LIT>" ] . v2d ( ) / ( gamma - <NUM_LIT:1.0> ) / rho . v ( ) <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> coeff = self . aux_data . get_var ( "<STR_LIT>" ) <EOL> coeff . v ( ) [ : , : ] = <NUM_LIT> / ( rho . v ( ) + rho_old . v ( ) ) <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) <EOL> self . aux_data . fill_BC ( "<STR_LIT>" ) <EOL> _ux , _vx , _uy , _vy = lm_interface_f . states ( myg . qx , myg . qy , myg . ng , <EOL> myg . dx , myg . dy , self . dt , <EOL> u . d , v . d , <EOL> ldelta_ux , ldelta_vx , <EOL> ldelta_uy , ldelta_vy , <EOL> coeff . d * gradp_x . d , coeff . d * gradp_y . d , <EOL> source . d , <EOL> u_MAC . d , v_MAC . d ) <EOL> u_xint = patch . ArrayIndexer ( d = _ux , grid = myg ) <EOL> v_xint = patch . ArrayIndexer ( d = _vx , grid = myg ) <EOL> u_yint = patch . ArrayIndexer ( d = _uy , grid = myg ) <EOL> v_yint = patch . ArrayIndexer ( d = _vy , grid = myg ) <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> advect_x = myg . scratch_array ( ) <EOL> advect_y = myg . scratch_array ( ) <EOL> advect_x . v ( ) [ : , : ] = <NUM_LIT:0.5> * ( u_MAC . v ( ) + u_MAC . ip ( <NUM_LIT:1> ) ) * ( u_xint . ip ( <NUM_LIT:1> ) - u_xint . v ( ) ) / myg . dx + <NUM_LIT:0.5> * ( v_MAC . v ( ) + v_MAC . jp ( <NUM_LIT:1> ) ) * ( u_yint . jp ( <NUM_LIT:1> ) - u_yint . v ( ) ) / myg . dy <EOL> advect_y . v ( ) [ : , : ] = <NUM_LIT:0.5> * ( u_MAC . v ( ) + u_MAC . ip ( <NUM_LIT:1> ) ) * ( v_xint . ip ( <NUM_LIT:1> ) - v_xint . v ( ) ) / myg . dx + <NUM_LIT:0.5> * ( v_MAC . v ( ) + v_MAC . jp ( <NUM_LIT:1> ) ) * ( v_yint . jp ( <NUM_LIT:1> ) - v_yint . v ( ) ) / myg . dy <EOL> proj_type = self . rp . get_param ( "<STR_LIT>" ) <EOL> if proj_type == <NUM_LIT:1> : <EOL> u . v ( ) [ : , : ] -= ( self . dt * advect_x . v ( ) + self . dt * gradp_x . v ( ) ) <EOL> v . v ( ) [ : , : ] -= ( self . dt * advect_y . v ( ) + self . dt * gradp_y . v ( ) ) <EOL> elif proj_type == <NUM_LIT:2> : <EOL> u . v ( ) [ : , : ] -= self . dt * advect_x . v ( ) <EOL> v . v ( ) [ : , : ] -= self . dt * advect_y . v ( ) <EOL> rho_half = <NUM_LIT:0.5> * ( rho + rho_old ) <EOL> rhoprime = self . make_prime ( rho_half , rho0 ) <EOL> source . d [ : , : ] = ( rhoprime * g / rho_half ) . d <EOL> self . aux_data . fill_BC ( "<STR_LIT>" ) <EOL> v . d [ : , : ] += self . dt * source . d <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> if self . verbose > <NUM_LIT:0> : <EOL> print ( "<STR_LIT>" . format ( self . cc_data . min ( "<STR_LIT>" ) , self . cc_data . max ( "<STR_LIT>" ) ) ) <EOL> print ( "<STR_LIT>" . format ( self . cc_data . min ( "<STR_LIT>" ) , self . cc_data . max ( "<STR_LIT>" ) ) ) <EOL> print ( "<STR_LIT>" . format ( self . cc_data . min ( "<STR_LIT>" ) , self . cc_data . max ( "<STR_LIT>" ) ) ) <EOL> if self . verbose > <NUM_LIT:0> : print ( "<STR_LIT>" ) <EOL> coeff = <NUM_LIT:1.0> / rho <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) ** <NUM_LIT:2> <EOL> mg = vcMG . VarCoeffCCMG2d ( myg . nx , myg . ny , <EOL> xl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xlb , <EOL> xr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . xrb , <EOL> yl_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . ylb , <EOL> yr_BC_type = self . cc_data . BCs [ "<STR_LIT>" ] . yrb , <EOL> xmin = myg . xmin , xmax = myg . xmax , <EOL> ymin = myg . ymin , ymax = myg . ymax , <EOL> coeffs = coeff , <EOL> coeffs_bc = self . cc_data . BCs [ "<STR_LIT>" ] , <EOL> verbose = <NUM_LIT:0> ) <EOL> div_beta_U . v ( ) [ : , : ] = <NUM_LIT:0.5> * beta0 . v2d ( ) * ( u . ip ( <NUM_LIT:1> ) - u . ip ( - <NUM_LIT:1> ) ) / myg . dx + <NUM_LIT:0.5> * ( beta0 . v2dp ( <NUM_LIT:1> ) * v . jp ( <NUM_LIT:1> ) - beta0 . v2dp ( - <NUM_LIT:1> ) * v . jp ( - <NUM_LIT:1> ) ) / myg . dy <EOL> mg . init_RHS ( div_beta_U . d / self . dt ) <EOL> phiGuess = mg . soln_grid . scratch_array ( ) <EOL> phiGuess . v ( buf = <NUM_LIT:1> ) [ : , : ] = phi . v ( buf = <NUM_LIT:1> ) <EOL> mg . init_solution ( phiGuess . d ) <EOL> mg . solve ( rtol = <NUM_LIT> ) <EOL> phi . d [ : , : ] = mg . get_solution ( grid = myg ) . d <EOL> gradphi_x , gradphi_y = mg . get_solution_gradient ( grid = myg ) <EOL> coeff = <NUM_LIT:1.0> / rho <EOL> coeff . v ( ) [ : , : ] = coeff . v ( ) * beta0 . v2d ( ) <EOL> u . v ( ) [ : , : ] -= self . dt * coeff . v ( ) * gradphi_x . v ( ) <EOL> v . v ( ) [ : , : ] -= self . dt * coeff . v ( ) * gradphi_y . v ( ) <EOL> if proj_type == <NUM_LIT:1> : <EOL> gradp_x . v ( ) [ : , : ] += gradphi_x . v ( ) <EOL> gradp_y . v ( ) [ : , : ] += gradphi_y . v ( ) <EOL> elif proj_type == <NUM_LIT:2> : <EOL> gradp_x . v ( ) [ : , : ] = gradphi_x . v ( ) <EOL> gradp_y . v ( ) [ : , : ] = gradphi_y . v ( ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> self . cc_data . fill_BC ( "<STR_LIT>" ) <EOL> if not self . in_preevolve : <EOL> self . cc_data . t += self . dt <EOL> self . n += <NUM_LIT:1> <EOL> def dovis ( self ) : <EOL> """<STR_LIT>""" <EOL> plt . clf ( ) <EOL> rho = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> rho0 = self . base [ "<STR_LIT>" ] <EOL> rhoprime = self . make_prime ( rho , rho0 ) <EOL> u = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> v = self . cc_data . get_var ( "<STR_LIT>" ) <EOL> myg = self . cc_data . grid <EOL> magvel = np . sqrt ( u ** <NUM_LIT:2> + v ** <NUM_LIT:2> ) <EOL> vort = myg . scratch_array ( ) <EOL> dv = <NUM_LIT:0.5> * ( v . ip ( <NUM_LIT:1> ) - v . ip ( - <NUM_LIT:1> ) ) / myg . dx <EOL> du = <NUM_LIT:0.5> * ( u . jp ( <NUM_LIT:1> ) - u . jp ( - <NUM_LIT:1> ) ) / myg . dy <EOL> vort . v ( ) [ : , : ] = dv - du <EOL> fig , axes = plt . subplots ( nrows = <NUM_LIT:2> , ncols = <NUM_LIT:2> , num = <NUM_LIT:1> ) <EOL> plt . subplots_adjust ( hspace = <NUM_LIT> ) <EOL> fields = [ rho , magvel , vort , rhoprime ] <EOL> field_names = [ r"<STR_LIT>" , r"<STR_LIT>" , r"<STR_LIT>" , r"<STR_LIT>" ] <EOL> for n in range ( len ( fields ) ) : <EOL> ax = axes . flat [ n ] <EOL> f = fields [ n ] <EOL> img = ax . imshow ( np . transpose ( f . v ( ) ) , <EOL> interpolation = "<STR_LIT>" , origin = "<STR_LIT>" , <EOL> extent = [ myg . xmin , myg . xmax , myg . ymin , myg . ymax ] ) <EOL> ax . set_xlabel ( "<STR_LIT:x>" ) <EOL> ax . set_ylabel ( "<STR_LIT:y>" ) <EOL> ax . set_title ( field_names [ n ] ) <EOL> plt . figtext ( <NUM_LIT> , <NUM_LIT> , "<STR_LIT>" % self . cc_data . t ) <EOL> plt . draw ( ) </s>
