diff --git a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py index b3e6bffd..abf040f8 100644 --- a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py +++ b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py @@ -4,6 +4,411 @@ import arrayfire as af import numpy as np from numpy.fft import fftfreq +from bolt.lib.nonlinear_solver.FVM_solver.reconstruct import reconstruct +import params + +class poisson_eqn_3D(object): + """ + This user class is an application context for the problem at hand; + It contains some parameters and frames the matrix system depending on + the system state. The Poisson object is used by the + compute_electrostatic_fields function in computing the electrostatic fields + using the PETSc's SNES solver methods + """ + def __init__(self, nonlinear_solver_obj): + self.da_3D = nonlinear_solver_obj._da_snes + self.da_2D = nonlinear_solver_obj._da_f + self.obj = nonlinear_solver_obj + self.glob_phi = self.da_3D.createGlobalVec() + self.local_phi = self.da_3D.createLocalVec() # phi with ghost zones + + self.glob_residual = self.da_3D.createGlobalVec() + + self.N_ghost = self.obj.N_ghost + self.N_ghost_poisson = self.obj.N_ghost_poisson + self.dq1 = self.obj.dq1 + self.dq2 = self.obj.dq2 + self.dq3 = self.obj.dq3 + + ((i_q1_2D_start, i_q2_2D_start), + (N_q1_2D_local, N_q2_2D_local) + ) = self.da_2D.getCorners() + + ((i_q1_3D_start, i_q2_3D_start, i_q3_3D_start), + (N_q1_3D_local, N_q2_3D_local, N_q3_3D_local) + ) = self.da_3D.getCorners() + + self.i_q1_3D_start = i_q1_3D_start + self.i_q2_3D_start = i_q2_3D_start + self.i_q3_3D_start = i_q3_3D_start + + self.i_q1_3D_end = i_q1_3D_start + N_q1_3D_local - 1 + self.i_q2_3D_end = i_q2_3D_start + N_q2_3D_local - 1 + self.i_q3_3D_end = i_q3_3D_start + N_q3_3D_local - 1 + + self.N_q1_2D_local = N_q1_2D_local + self.N_q2_2D_local = N_q2_2D_local + + self.N_q1_3D_local = N_q1_3D_local + self.N_q2_3D_local = N_q2_3D_local + self.N_q3_3D_local = N_q3_3D_local + + location_in_q3 = self.obj.location_in_q3 + N_g = self.N_ghost # Ghost zones of Boltzmann solver + N_gp = self.N_ghost_poisson # Ghost zones of Poisson solver + + self.density_np = np.zeros([N_q2_2D_local + 2*N_g, + N_q1_2D_local + 2*N_g + ] + ) + # Cell centers in 3D + i_q1_3D = ( (i_q1_3D_start + 0.5) + + np.arange(-N_gp, N_q1_3D_local + N_gp) + ) + + i_q2_3D = ( (i_q2_3D_start + 0.5) + + np.arange(-N_gp, N_q2_3D_local + N_gp) + ) + + i_q3_3D = ( (i_q3_3D_start + 0.5) + + np.arange(-N_gp, N_q3_3D_local + N_gp) + ) + + q1_2D_start = self.obj.q1_start + q1_2D_end = self.obj.q1_end + q2_2D_start = self.obj.q2_start + q2_2D_end = self.obj.q2_end + q3_3D_start = self.obj.q3_3D_start + + # TODO: Code below is duplication of _calculate_q_center in nonlinear_solver.py + i_q1_2D = ( (i_q1_2D_start + 0.5) + + np.arange(-N_g, N_q1_2D_local + N_g) + ) + + i_q2_2D = ( (i_q2_2D_start + 0.5) + + np.arange(-N_g, N_q2_2D_local + N_g) + ) + + q1_2D = q1_2D_start + i_q1_2D * self.dq1 + q2_2D = q2_2D_start + i_q2_2D * self.dq2 + + length_q1_2d = (q1_2D_end - q1_2D_start) + length_q2_2d = (q2_2D_end - q2_2D_start) + + length_multiples_q1 = self.obj.length_multiples_q1 + length_multiples_q2 = self.obj.length_multiples_q2 + + q1_3D = q1_2D_start - length_multiples_q1*length_q1_2d + i_q1_3D * self.dq1 + q2_3D = q2_2D_start - length_multiples_q2*length_q2_2d + i_q2_3D * self.dq2 + q3_3D = q3_3D_start + i_q3_3D * self.dq3 + + self.q1_3D = q1_3D + self.q2_3D = q2_3D + self.q3_3D = q3_3D + + glob_epsilon = self.da_3D.createGlobalVec() + local_epsilon = self.da_3D.createLocalVec() + + epsilon_array = local_epsilon.getArray(readonly=0) + epsilon_array = epsilon_array.reshape([N_q3_3D_local + 2*N_gp, \ + N_q2_3D_local + 2*N_gp, \ + N_q1_3D_local + 2*N_gp + ] + ) + epsilon_array[:] = params.epsilon0 + epsilon_array[q3_3D= q1_2D[N_g] - tol) \ + & (q1_3D_data_structure <= q1_2D[-1-N_g] + tol) \ + & (q2_3D_data_structure >= q2_2D[N_g] - tol) \ + & (q2_3D_data_structure <= q2_2D[-1-N_g] + tol) \ + & (q3_3D_data_structure >= location_in_q3) \ + & (q3_3D_data_structure < location_in_q3 + self.dq3) + + self.cond_2D = (q1_2D_data_structure >= q1_2D[N_g]) \ + & (q1_2D_data_structure <= q1_2D[-1-N_g]) \ + & (q2_2D_data_structure >= q2_2D[N_g]) \ + & (q2_2D_data_structure <= q2_2D[-1-N_g]) + + self.cond_3D_phi = (q1_3D_data_structure >= q1_2D[0] - tol) \ + & (q1_3D_data_structure <= q1_2D[-1] + tol) \ + & (q2_3D_data_structure >= q2_2D[0] - tol) \ + & (q2_3D_data_structure <= q2_2D[-1] + tol) \ + & (q3_3D_data_structure >= location_in_q3) \ + & (q3_3D_data_structure < location_in_q3 + self.dq3) + + contact_start = params.contact_start + contact_end = params.contact_end + + self.cond_left_contact = \ + (q1_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] < q1_2D[N_g] - tol) \ + & (q2_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] > contact_start-tol) \ + & (q2_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] < contact_end + tol) \ + & (q3_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] > location_in_q3) \ + & (q3_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] < location_in_q3 + self.dq3) + + self.cond_right_contact = \ + (q1_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] > q1_2D[-1-N_g] + tol) \ + & (q2_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] > contact_start-tol) \ + & (q2_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] < contact_end + tol) \ + & (q3_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] > location_in_q3) \ + & (q3_3D_data_structure[N_gp:-N_gp, N_gp:-N_gp, N_gp:-N_gp] < location_in_q3 + self.dq3) + + backgate_potential = params.backgate_potential + + self.q3 = q3_3D_data_structure + z = self.q3 + #z_sample = q3_3D[self.q3_2D_in_3D_index_start] + #TODO: FIX THIS ASAP + z_sample = location_in_q3 + z_backgate = q3_3D[N_gp] + side_wall_boundaries = \ + backgate_potential*(z_sample - z)/(z_sample - z_backgate) + + self.bc = 0.*self.q3 # 3D boundary condition array + + self.bc[:] = 0. + self.bc[:N_gp, :, :] = backgate_potential # backgate + + below_sample = q3_3Daf + self.phi = af.to_array(phi_2D_local_array) + # Since rho was defined at (i + 0.5, j + 0.5) + # Electric Potential returned will also be at (i + 0.5, j + 0.5) + self.phi = af.moddims(self.phi, + self.N_q1_local + 2*N_g, + self.N_q2_local + 2*N_g + ) + params.phi = self.phi + + # Obtaining the electric field values at (i+0.5, j+0.5) by first + # reconstructing phi to the edges and then differencing them. Needed + # in case phi develops large gradients + + method_in_q = self.physical_system.params.reconstruction_method_in_q + + phi_left, phi_right = reconstruct(self, self.phi, 0, method_in_q) + E1 = -(phi_right - phi_left)/self.dq1 + + phi_bottom, phi_top = reconstruct(self, self.phi, 1, method_in_q) + E2 = -(phi_top - phi_bottom)/self.dq2 + + af.eval(E1, E2) + + E1 = af.moddims(E1, + 1, + self.N_q1_local + 2*N_g, + self.N_q2_local + 2*N_g + ) + + E2 = af.moddims(E2, + 1, + self.N_q1_local + 2*N_g, + self.N_q2_local + 2*N_g + ) + + self.cell_centered_EM_fields[0, :] = E1 + self.cell_centered_EM_fields[1, :] = E2 + self.cell_centered_EM_fields[2, :] = 0. # TODO: worry about this later + + if(self.performance_test_flag == True): + af.sync() + toc = af.time() + self.time_fieldsolver += toc - tic + + return + def fft_poisson(self, f=None): """ diff --git a/bolt/lib/nonlinear_solver/EM_fields_solver/fields_step.py b/bolt/lib/nonlinear_solver/EM_fields_solver/fields_step.py index 95ed59f9..9e325bff 100644 --- a/bolt/lib/nonlinear_solver/EM_fields_solver/fields_step.py +++ b/bolt/lib/nonlinear_solver/EM_fields_solver/fields_step.py @@ -3,7 +3,7 @@ import arrayfire as af -from .electrostatic import fft_poisson +from .electrostatic import fft_poisson, compute_electrostatic_fields from .fdtd_explicit import fdtd, fdtd_grid_to_ck_grid from .. import interpolation_routines @@ -12,12 +12,17 @@ def fields_step(self, dt): if(self.performance_test_flag == True): tic = af.time() - if (self.physical_system.params.fields_solver == 'fft'): - fft_poisson(self) + if (self.physical_system.params.fields_type == 'electrostatic'): + if (self.physical_system.params.fields_solver == 'fft'): + fft_poisson(self) + elif (self.physical_system.params.fields_solver == 'SNES'): + #compute_electrostatic_fields(self) + pass + self._communicate_fields() self._apply_bcs_fields() - elif (self.physical_system.params.fields_solver == 'fdtd'): + elif (self.physical_system.params.fields_type == 'electrodynamic'): # Will return a flattened array containing the values of # J1,2,3 in 2D space: self.J1 = self.physical_system.params.charge_electron \ diff --git a/bolt/lib/nonlinear_solver/FVM_solver/df_dt_fvm.py b/bolt/lib/nonlinear_solver/FVM_solver/df_dt_fvm.py index fdd97147..e43b3c89 100644 --- a/bolt/lib/nonlinear_solver/FVM_solver/df_dt_fvm.py +++ b/bolt/lib/nonlinear_solver/FVM_solver/df_dt_fvm.py @@ -1,9 +1,11 @@ import arrayfire as af from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic import fft_poisson +from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic import compute_electrostatic_fields # Importing Riemann solver used in calculating fluxes: -from .riemann_solver import riemann_solver +from .riemann_solver import riemann_solver, upwind_flux from .reconstruct import reconstruct +import params # Equation to solve: # df/dt + d(C_q1 * f)/dq1 + d(C_q2 * f)/dq2 = C[f] @@ -75,11 +77,25 @@ def df_dt_fvm(f, self, at_n = True): if( self.physical_system.params.solver_method_in_p == 'FVM' and self.physical_system.params.charge_electron != 0 ): - if(self.physical_system.params.fields_solver == 'fft'): + if(self.physical_system.params.fields_type == 'electrostatic'): - fft_poisson(self, f) - self._communicate_fields() - self._apply_bcs_fields() + if (params.time_step%params.electrostatic_solver_every_nth_step==0): + + if(self.physical_system.params.fields_solver == 'fft'): + + fft_poisson(self, f) + + elif(self.physical_system.params.fields_solver == 'SNES'): + compute_electrostatic_fields(self) + pass + + E1 = self.cell_centered_EM_fields[0] + E2 = self.cell_centered_EM_fields[1] + E3 = self.cell_centered_EM_fields[2] + + B1 = self.cell_centered_EM_fields[3] + B2 = self.cell_centered_EM_fields[4] + B3 = self.cell_centered_EM_fields[5] # This is taken care of by the timestepper that is utilized # when FDTD is to be used with FVM in p-space @@ -91,39 +107,39 @@ def df_dt_fvm(f, self, at_n = True): invalid/not-implemented' ) - if( self.physical_system.params.fields_solver == 'fdtd' - and at_n == True - ): - - E1 = self.cell_centered_EM_fields_at_n[0] - E2 = self.cell_centered_EM_fields_at_n[1] - E3 = self.cell_centered_EM_fields_at_n[2] - - B1 = self.cell_centered_EM_fields_at_n[3] - B2 = self.cell_centered_EM_fields_at_n[4] - B3 = self.cell_centered_EM_fields_at_n[5] - - elif( self.physical_system.params.fields_solver == 'fdtd' - and at_n != False - ): - - E1 = self.cell_centered_EM_fields_at_n_plus_half[0] - E2 = self.cell_centered_EM_fields_at_n_plus_half[1] - E3 = self.cell_centered_EM_fields_at_n_plus_half[2] - - B1 = self.cell_centered_EM_fields_at_n_plus_half[3] - B2 = self.cell_centered_EM_fields_at_n_plus_half[4] - B3 = self.cell_centered_EM_fields_at_n_plus_half[5] - - else: - - E1 = self.cell_centered_EM_fields[0] - E2 = self.cell_centered_EM_fields[1] - E3 = self.cell_centered_EM_fields[2] - - B1 = self.cell_centered_EM_fields[3] - B2 = self.cell_centered_EM_fields[4] - B3 = self.cell_centered_EM_fields[5] +# if( self.physical_system.params.fields_solver == 'fdtd' +# and at_n == True +# ): +# +# E1 = self.cell_centered_EM_fields_at_n[0] +# E2 = self.cell_centered_EM_fields_at_n[1] +# E3 = self.cell_centered_EM_fields_at_n[2] +# +# B1 = self.cell_centered_EM_fields_at_n[3] +# B2 = self.cell_centered_EM_fields_at_n[4] +# B3 = self.cell_centered_EM_fields_at_n[5] +# +# elif( self.physical_system.params.fields_solver == 'fdtd' +# and at_n != False +# ): +# +# E1 = self.cell_centered_EM_fields_at_n_plus_half[0] +# E2 = self.cell_centered_EM_fields_at_n_plus_half[1] +# E3 = self.cell_centered_EM_fields_at_n_plus_half[2] +# +# B1 = self.cell_centered_EM_fields_at_n_plus_half[3] +# B2 = self.cell_centered_EM_fields_at_n_plus_half[4] +# B3 = self.cell_centered_EM_fields_at_n_plus_half[5] +# +# else: +# +# E1 = self.cell_centered_EM_fields[0] +# E2 = self.cell_centered_EM_fields[1] +# E3 = self.cell_centered_EM_fields[2] +# +# B1 = self.cell_centered_EM_fields[3] +# B2 = self.cell_centered_EM_fields[4] +# B3 = self.cell_centered_EM_fields[5] (A_p1, A_p2, A_p3) = af.broadcast(self._A_p, self.q1_center, self.q2_center, self.p1, self.p2, self.p3, @@ -133,26 +149,42 @@ def df_dt_fvm(f, self, at_n = True): # Variation of p1 is along axis 0: left_plus_eps_flux_p1, right_minus_eps_flux_p1 = \ - reconstruct(self, self._convert_to_p_expanded(af.broadcast(multiply, A_p1, f)), 0, method_in_p) + reconstruct(self, + self._convert_to_p_expanded(af.broadcast(multiply, A_p1, f)), + 0, method_in_p + ) # Variation of p2 is along axis 1: bot_plus_eps_flux_p2, top_minus_eps_flux_p2 = \ - reconstruct(self, self._convert_to_p_expanded(af.broadcast(multiply, A_p2, f)), 1, method_in_p) - # Variation of p3 is along axis 2: - back_plus_eps_flux_p3, front_minus_eps_flux_p3 = \ - reconstruct(self, self._convert_to_p_expanded(af.broadcast(multiply, A_p3, f)), 2, method_in_p) + reconstruct(self, + self._convert_to_p_expanded(af.broadcast(multiply, A_p2, f)), + 1, method_in_p + ) +# # Variation of p3 is along axis 2: +# back_plus_eps_flux_p3, front_minus_eps_flux_p3 = \ +# reconstruct(self, self._convert_to_p_expanded(af.broadcast(multiply, A_p3, f)), 2, method_in_p) left_minus_eps_flux_p1 = af.shift(right_minus_eps_flux_p1, 1) bot_minus_eps_flux_p2 = af.shift(top_minus_eps_flux_p2, 0, 1) - back_minus_eps_flux_p3 = af.shift(front_minus_eps_flux_p3, 0, 0, 1) +# back_minus_eps_flux_p3 = af.shift(front_minus_eps_flux_p3, 0, 0, 1) # Obtaining the fluxes by face-averaging: - left_flux_p1 = 0.5 * (left_minus_eps_flux_p1 + left_plus_eps_flux_p1) - bot_flux_p2 = 0.5 * (bot_minus_eps_flux_p2 + bot_plus_eps_flux_p2) - back_flux_p3 = 0.5 * (back_minus_eps_flux_p3 + back_plus_eps_flux_p3) +# left_flux_p1 = 0.5 * (left_minus_eps_flux_p1 + left_plus_eps_flux_p1) +# bot_flux_p2 = 0.5 * (bot_minus_eps_flux_p2 + bot_plus_eps_flux_p2) + # back_flux_p3 = 0.5 * (back_minus_eps_flux_p3 + back_plus_eps_flux_p3) + + left_flux_p1 = upwind_flux(left_minus_eps_flux_p1, \ + left_plus_eps_flux_p1, \ + self._convert_to_p_expanded(A_p1) + ) + + bot_flux_p2 = upwind_flux(bot_minus_eps_flux_p2, \ + bot_plus_eps_flux_p2, \ + self._convert_to_p_expanded(A_p2) + ) right_flux_p1 = af.shift(left_flux_p1, -1) top_flux_p2 = af.shift(bot_flux_p2, 0, -1) - front_flux_p3 = af.shift(back_flux_p3, 0, 0, -1) +# front_flux_p3 = af.shift(back_flux_p3, 0, 0, -1) left_flux_p1 = self._convert_to_q_expanded(left_flux_p1) right_flux_p1 = self._convert_to_q_expanded(right_flux_p1) @@ -160,13 +192,12 @@ def df_dt_fvm(f, self, at_n = True): bot_flux_p2 = self._convert_to_q_expanded(bot_flux_p2) top_flux_p2 = self._convert_to_q_expanded(top_flux_p2) - back_flux_p3 = self._convert_to_q_expanded(back_flux_p3) - front_flux_p3 = self._convert_to_q_expanded(front_flux_p3) +# back_flux_p3 = self._convert_to_q_expanded(back_flux_p3) +# front_flux_p3 = self._convert_to_q_expanded(front_flux_p3) df_dt += - (right_flux_p1 - left_flux_p1)/self.dp1 \ - (top_flux_p2 - bot_flux_p2 )/self.dp2 \ - - (front_flux_p3 - back_flux_p3)/self.dp3 - + #- (front_flux_p3 - back_flux_p3)/self.dp3 af.eval(df_dt) return(df_dt) diff --git a/bolt/lib/nonlinear_solver/FVM_solver/reconstruction_methods/weno5.py b/bolt/lib/nonlinear_solver/FVM_solver/reconstruction_methods/weno5.py index 11814c23..28bf2b2a 100644 --- a/bolt/lib/nonlinear_solver/FVM_solver/reconstruction_methods/weno5.py +++ b/bolt/lib/nonlinear_solver/FVM_solver/reconstruction_methods/weno5.py @@ -75,5 +75,6 @@ def reconstruct_weno5(input_array, axis): # Reconstruction: left_value = (w1l * u1l + w2l * u2l + w3l * u3l) / denl; right_value = (w1r * u1r + w2r * u2r + w3r * u3r) / denr; - + + af.eval(left_value, right_value) return(left_value, right_value) diff --git a/bolt/lib/nonlinear_solver/FVM_solver/timestep_df_dt.py b/bolt/lib/nonlinear_solver/FVM_solver/timestep_df_dt.py index ff58735f..387c576d 100644 --- a/bolt/lib/nonlinear_solver/FVM_solver/timestep_df_dt.py +++ b/bolt/lib/nonlinear_solver/FVM_solver/timestep_df_dt.py @@ -11,59 +11,62 @@ def fvm_timestep_RK2(self, dt): + self._communicate_f() + self._apply_bcs_f() + f_initial = self.f self.f = self.f + df_dt_fvm(self.f, self, True) * (dt / 2) self._communicate_f() self._apply_bcs_f() - if( self.physical_system.params.charge_electron != 0 - and self.physical_system.params.fields_solver == 'fdtd' - ): - - # Will return a flattened array containing the values of - # J1,2,3 in 2D space: - self.J1 = self.physical_system.params.charge_electron \ - * self.compute_moments('mom_p1_bulk') # (i + 1/2, j + 1/2) - self.J2 = self.physical_system.params.charge_electron \ - * self.compute_moments('mom_p2_bulk') # (i + 1/2, j + 1/2) - self.J3 = self.physical_system.params.charge_electron \ - * self.compute_moments('mom_p3_bulk') # (i + 1/2, j + 1/2) - - # Obtaining the values for current density on the Yee-Grid: - self.J1 = 0.5 * (self.J1 + af.shift(self.J1, 0, 0, 1)) # (i + 1/2, j) - self.J2 = 0.5 * (self.J2 + af.shift(self.J2, 0, 1, 0)) # (i, j + 1/2) - - self.J3 = 0.25 * ( self.J3 + af.shift(self.J3, 0, 1, 0) + - + af.shift(self.J3, 0, 0, 1) - + af.shift(self.J3, 0, 1, 1) - ) # (i, j) - - # Here: - # cell_centered_EM_fields[:3] is at n - # cell_centered_EM_fields[3:] is at n+1/2 - # cell_centered_EM_fields_at_n_plus_half[3:] is at n-1/2 - - self.cell_centered_EM_fields_at_n[:3] = self.cell_centered_EM_fields[:3] - self.cell_centered_EM_fields_at_n[3:] = \ - 0.5 * ( self.cell_centered_EM_fields_at_n_plus_half[3:] - + self.cell_centered_EM_fields[3:] - ) - - - self.cell_centered_EM_fields_at_n_plus_half[3:] = self.cell_centered_EM_fields[3:] - - fdtd(self, dt) - fdtd_grid_to_ck_grid(self) - - # Here - # cell_centered_EM_fields[:3] is at n+1 - # cell_centered_EM_fields[3:] is at n+3/2 - - self.cell_centered_EM_fields_at_n_plus_half[:3] = \ - 0.5 * ( self.cell_centered_EM_fields_at_n_plus_half[:3] - + self.cell_centered_EM_fields[:3] - ) +# if( self.physical_system.params.charge_electron != 0 +# and self.physical_system.params.fields_solver == 'fdtd' +# ): +# +# # Will return a flattened array containing the values of +# # J1,2,3 in 2D space: +# self.J1 = self.physical_system.params.charge_electron \ +# * self.compute_moments('mom_p1_bulk') # (i + 1/2, j + 1/2) +# self.J2 = self.physical_system.params.charge_electron \ +# * self.compute_moments('mom_p2_bulk') # (i + 1/2, j + 1/2) +# self.J3 = self.physical_system.params.charge_electron \ +# * self.compute_moments('mom_p3_bulk') # (i + 1/2, j + 1/2) +# +# # Obtaining the values for current density on the Yee-Grid: +# self.J1 = 0.5 * (self.J1 + af.shift(self.J1, 0, 0, 1)) # (i + 1/2, j) +# self.J2 = 0.5 * (self.J2 + af.shift(self.J2, 0, 1, 0)) # (i, j + 1/2) +# +# self.J3 = 0.25 * ( self.J3 + af.shift(self.J3, 0, 1, 0) + +# + af.shift(self.J3, 0, 0, 1) +# + af.shift(self.J3, 0, 1, 1) +# ) # (i, j) +# +# # Here: +# # cell_centered_EM_fields[:3] is at n +# # cell_centered_EM_fields[3:] is at n+1/2 +# # cell_centered_EM_fields_at_n_plus_half[3:] is at n-1/2 +# +# self.cell_centered_EM_fields_at_n[:3] = self.cell_centered_EM_fields[:3] +# self.cell_centered_EM_fields_at_n[3:] = \ +# 0.5 * ( self.cell_centered_EM_fields_at_n_plus_half[3:] +# + self.cell_centered_EM_fields[3:] +# ) +# +# +# self.cell_centered_EM_fields_at_n_plus_half[3:] = self.cell_centered_EM_fields[3:] +# +# fdtd(self, dt) +# fdtd_grid_to_ck_grid(self) +# +# # Here +# # cell_centered_EM_fields[:3] is at n+1 +# # cell_centered_EM_fields[3:] is at n+3/2 +# +# self.cell_centered_EM_fields_at_n_plus_half[:3] = \ +# 0.5 * ( self.cell_centered_EM_fields_at_n_plus_half[:3] +# + self.cell_centered_EM_fields[:3] +# ) self.f = f_initial + df_dt_fvm(self.f, self, False) * dt diff --git a/bolt/lib/nonlinear_solver/apply_boundary_conditions.py b/bolt/lib/nonlinear_solver/apply_boundary_conditions.py index f2054d8f..ae5d0cfa 100644 --- a/bolt/lib/nonlinear_solver/apply_boundary_conditions.py +++ b/bolt/lib/nonlinear_solver/apply_boundary_conditions.py @@ -30,8 +30,12 @@ def apply_dirichlet_bcs_f(self, boundary): self.physical_system.params ) + # NOTE: This is not necessary since the Riemann solver ensures that the + # information outflow characteristics in the ghost zones do not affect + # the numerical solution inside in the physical domain. + # # Only changing inflowing characteristics: - f_left = af.select(A_q1>0, f_left, self.f) + #f_left = af.select(A_q1>0, f_left, self.f) self.f[:, :N_g] = f_left[:, :N_g] @@ -43,7 +47,7 @@ def apply_dirichlet_bcs_f(self, boundary): ) # Only changing inflowing characteristics: - f_right = af.select(A_q1<0, f_right, self.f) + #f_right = af.select(A_q1<0, f_right, self.f) self.f[:, -N_g:] = f_right[:, -N_g:] @@ -55,7 +59,7 @@ def apply_dirichlet_bcs_f(self, boundary): ) # Only changing inflowing characteristics: - f_bottom = af.select(A_q2>0, f_bottom, self.f) + #f_bottom = af.select(A_q2>0, f_bottom, self.f) self.f[:, :, :N_g] = f_bottom[:, :, :N_g] @@ -67,7 +71,7 @@ def apply_dirichlet_bcs_f(self, boundary): ) # Only changing inflowing characteristics: - f_top = af.select(A_q2<0, f_top, self.f) + #f_top = af.select(A_q2<0, f_top, self.f) self.f[:, :, -N_g:] = f_top[:, :, -N_g:] @@ -266,146 +270,146 @@ def apply_bcs_f(self): def apply_dirichlet_bcs_fields(self, boundary): - N_g = self.N_ghost - - # These arguments are defined since they are required by all the function calls: - # So the functions can be called instead using function(*args) - args = (self.q1_center, self.q2_center, self.physical_system.params) - - if(boundary == 'left'): - E1 = self.boundary_conditions.\ - E1_left(self.cell_centered_EM_fields[0],*args)[:, :N_g] - - E2 = self.boundary_conditions.\ - E2_left(self.cell_centered_EM_fields[1],*args)[:, :N_g] - - E3 = self.boundary_conditions.\ - E3_left(self.cell_centered_EM_fields[2],*args)[:, :N_g] - - B1 = self.boundary_conditions.\ - B1_left(self.cell_centered_EM_fields[3],*args)[:, :N_g] - - B2 = self.boundary_conditions.\ - B2_left(self.cell_centered_EM_fields[4],*args)[:, :N_g] - - B3 = self.boundary_conditions.\ - B3_left(self.cell_centered_EM_fields[5],*args)[:, :N_g] - - self.cell_centered_EM_fields[:, :N_g] = \ - af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) - - elif(boundary == 'right'): - E1 = self.boundary_conditions.\ - E1_right(self.cell_centered_EM_fields[0],*args)[:, -N_g:] - - E2 = self.boundary_conditions.\ - E2_right(self.cell_centered_EM_fields[1],*args)[:, -N_g:] - - E3 = self.boundary_conditions.\ - E3_right(self.cell_centered_EM_fields[2],*args)[:, -N_g:] - - B1 = self.boundary_conditions.\ - B1_right(self.cell_centered_EM_fields[3],*args)[:, -N_g:] - - B2 = self.boundary_conditions.\ - B2_right(self.cell_centered_EM_fields[4],*args)[:, -N_g:] - - B3 = self.boundary_conditions.\ - B3_right(self.cell_centered_EM_fields[5],*args)[:, -N_g:] - - self.cell_centered_EM_fields[:, -N_g:] = \ - af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) - - elif(boundary == 'bottom'): - E1 = self.boundary_conditions.\ - E1_bottom(self.cell_centered_EM_fields[0],*args)[:, :, :N_g] - - E2 = self.boundary_conditions.\ - E2_bottom(self.cell_centered_EM_fields[1],*args)[:, :, :N_g] - - E3 = self.boundary_conditions.\ - E3_bottom(self.cell_centered_EM_fields[2],*args)[:, :, :N_g] - - B1 = self.boundary_conditions.\ - B1_bottom(self.cell_centered_EM_fields[3],*args)[:, :, :N_g] - - B2 = self.boundary_conditions.\ - B2_bottom(self.cell_centered_EM_fields[4],*args)[:, :, :N_g] - - B3 = self.boundary_conditions.\ - B3_bottom(self.cell_centered_EM_fields[5],*args)[:, :, :N_g] - - self.cell_centered_EM_fields[:, :, :N_g] = \ - af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) - - elif(boundary == 'top'): - E1 = self.boundary_conditions.\ - E1_top(self.cell_centered_EM_fields[0],*args)[:, :, -N_g:] - - E2 = self.boundary_conditions.\ - E2_top(self.cell_centered_EM_fields[1],*args)[:, :, -N_g:] - - E3 = self.boundary_conditions.\ - E3_top(self.cell_centered_EM_fields[2],*args)[:, :, -N_g:] - - B1 = self.boundary_conditions.\ - B1_top(self.cell_centered_EM_fields[3],*args)[:, :, -N_g:] - - B2 = self.boundary_conditions.\ - B2_top(self.cell_centered_EM_fields[4],*args)[:, :, -N_g:] - - B3 = self.boundary_conditions.\ - B3_top(self.cell_centered_EM_fields[5],*args)[:, :, -N_g:] - - self.cell_centered_EM_fields[:, :, -N_g:] = \ - af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) - - else: - raise Exception('Invalid choice for boundary') +# N_g = self.N_ghost +# +# # These arguments are defined since they are required by all the function calls: +# # So the functions can be called instead using function(*args) +# args = (self.q1_center, self.q2_center, self.physical_system.params) +# +# if(boundary == 'left'): +# E1 = self.boundary_conditions.\ +# E1_left(self.cell_centered_EM_fields[0],*args)[:, :N_g] +# +# E2 = self.boundary_conditions.\ +# E2_left(self.cell_centered_EM_fields[1],*args)[:, :N_g] +# +# E3 = self.boundary_conditions.\ +# E3_left(self.cell_centered_EM_fields[2],*args)[:, :N_g] +# +# B1 = self.boundary_conditions.\ +# B1_left(self.cell_centered_EM_fields[3],*args)[:, :N_g] +# +# B2 = self.boundary_conditions.\ +# B2_left(self.cell_centered_EM_fields[4],*args)[:, :N_g] +# +# B3 = self.boundary_conditions.\ +# B3_left(self.cell_centered_EM_fields[5],*args)[:, :N_g] +# +# self.cell_centered_EM_fields[:, :N_g] = \ +# af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) +# +# elif(boundary == 'right'): +# E1 = self.boundary_conditions.\ +# E1_right(self.cell_centered_EM_fields[0],*args)[:, -N_g:] +# +# E2 = self.boundary_conditions.\ +# E2_right(self.cell_centered_EM_fields[1],*args)[:, -N_g:] +# +# E3 = self.boundary_conditions.\ +# E3_right(self.cell_centered_EM_fields[2],*args)[:, -N_g:] +# +# B1 = self.boundary_conditions.\ +# B1_right(self.cell_centered_EM_fields[3],*args)[:, -N_g:] +# +# B2 = self.boundary_conditions.\ +# B2_right(self.cell_centered_EM_fields[4],*args)[:, -N_g:] +# +# B3 = self.boundary_conditions.\ +# B3_right(self.cell_centered_EM_fields[5],*args)[:, -N_g:] +# +# self.cell_centered_EM_fields[:, -N_g:] = \ +# af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) +# +# elif(boundary == 'bottom'): +# E1 = self.boundary_conditions.\ +# E1_bottom(self.cell_centered_EM_fields[0],*args)[:, :, :N_g] +# +# E2 = self.boundary_conditions.\ +# E2_bottom(self.cell_centered_EM_fields[1],*args)[:, :, :N_g] +# +# E3 = self.boundary_conditions.\ +# E3_bottom(self.cell_centered_EM_fields[2],*args)[:, :, :N_g] +# +# B1 = self.boundary_conditions.\ +# B1_bottom(self.cell_centered_EM_fields[3],*args)[:, :, :N_g] +# +# B2 = self.boundary_conditions.\ +# B2_bottom(self.cell_centered_EM_fields[4],*args)[:, :, :N_g] +# +# B3 = self.boundary_conditions.\ +# B3_bottom(self.cell_centered_EM_fields[5],*args)[:, :, :N_g] +# +# self.cell_centered_EM_fields[:, :, :N_g] = \ +# af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) +# +# elif(boundary == 'top'): +# E1 = self.boundary_conditions.\ +# E1_top(self.cell_centered_EM_fields[0],*args)[:, :, -N_g:] +# +# E2 = self.boundary_conditions.\ +# E2_top(self.cell_centered_EM_fields[1],*args)[:, :, -N_g:] +# +# E3 = self.boundary_conditions.\ +# E3_top(self.cell_centered_EM_fields[2],*args)[:, :, -N_g:] +# +# B1 = self.boundary_conditions.\ +# B1_top(self.cell_centered_EM_fields[3],*args)[:, :, -N_g:] +# +# B2 = self.boundary_conditions.\ +# B2_top(self.cell_centered_EM_fields[4],*args)[:, :, -N_g:] +# +# B3 = self.boundary_conditions.\ +# B3_top(self.cell_centered_EM_fields[5],*args)[:, :, -N_g:] +# +# self.cell_centered_EM_fields[:, :, -N_g:] = \ +# af.join(0, E1, E2, E3, af.join(0, B1, B2, B3)) +# +# else: +# raise Exception('Invalid choice for boundary') return def apply_mirror_bcs_fields(self, boundary): - N_g = self.N_ghost - - if(boundary == 'left'): - # x-0-x-0-x-0-|-0-x-0-x-0-x-.... - # 0 1 2 3 4 5 - # For mirror boundary conditions: - # 0 = 5; 1 = 4; 2 = 3; - self.cell_centered_EM_fields[:, :N_g] = \ - af.flip(self.cell_centered_EM_fields[:, N_g:2 * N_g], 1) - - elif(boundary == 'right'): - # ...-x-0-x-0-x-0-|-0-x-0-x-0-x - # -6 -5 -4 -3 -2 -1 - # For mirror boundary conditions: - # -1 = -6; -2 = -5; -3 = -4; - - self.cell_centered_EM_fields[:, -N_g:] = \ - af.flip(self.cell_centered_EM_fields[:, -2 * N_g:-N_g], 1) - - elif(boundary == 'bottom'): - # x-0-x-0-x-0-|-0-x-0-x-0-x-.... - # 0 1 2 3 4 5 - # For mirror boundary conditions: - # 0 = 5; 1 = 4; 2 = 3; - - self.cell_centered_EM_fields[:, :, :N_g] = \ - af.flip(self.cell_centered_EM_fields[:, :, N_g:2 * N_g], 2) - - elif(boundary == 'top'): - # ...-x-0-x-0-x-0-|-0-x-0-x-0-x - # -6 -5 -4 -3 -2 -1 - # For mirror boundary conditions: - # -1 = -6; -2 = -5; -3 = -4; - - self.cell_centered_EM_fields[:, :, -N_g:] = \ - af.flip(self.cell_centered_EM_fields[:, :, -2 * N_g:-N_g], 2) - - else: - raise Exception('Invalid choice for boundary') +# N_g = self.N_ghost +# +# if(boundary == 'left'): +# # x-0-x-0-x-0-|-0-x-0-x-0-x-.... +# # 0 1 2 3 4 5 +# # For mirror boundary conditions: +# # 0 = 5; 1 = 4; 2 = 3; +# self.cell_centered_EM_fields[:, :N_g] = \ +# af.flip(self.cell_centered_EM_fields[:, N_g:2 * N_g], 1) +# +# elif(boundary == 'right'): +# # ...-x-0-x-0-x-0-|-0-x-0-x-0-x +# # -6 -5 -4 -3 -2 -1 +# # For mirror boundary conditions: +# # -1 = -6; -2 = -5; -3 = -4; +# +# self.cell_centered_EM_fields[:, -N_g:] = \ +# af.flip(self.cell_centered_EM_fields[:, -2 * N_g:-N_g], 1) +# +# elif(boundary == 'bottom'): +# # x-0-x-0-x-0-|-0-x-0-x-0-x-.... +# # 0 1 2 3 4 5 +# # For mirror boundary conditions: +# # 0 = 5; 1 = 4; 2 = 3; +# +# self.cell_centered_EM_fields[:, :, :N_g] = \ +# af.flip(self.cell_centered_EM_fields[:, :, N_g:2 * N_g], 2) +# +# elif(boundary == 'top'): +# # ...-x-0-x-0-x-0-|-0-x-0-x-0-x +# # -6 -5 -4 -3 -2 -1 +# # For mirror boundary conditions: +# # -1 = -6; -2 = -5; -3 = -4; +# +# self.cell_centered_EM_fields[:, :, -N_g:] = \ +# af.flip(self.cell_centered_EM_fields[:, :, -2 * N_g:-N_g], 2) +# +# else: +# raise Exception('Invalid choice for boundary') return diff --git a/bolt/lib/nonlinear_solver/file_io/dump.py b/bolt/lib/nonlinear_solver/file_io/dump.py index 9b67706b..0ec6484a 100644 --- a/bolt/lib/nonlinear_solver/file_io/dump.py +++ b/bolt/lib/nonlinear_solver/file_io/dump.py @@ -5,6 +5,37 @@ import numpy as np import arrayfire as af +def dump_aux_arrays(self, arrays, name, file_name): + + if (self.dump_aux_arrays_initial_call): + self._da_aux_arrays = PETSc.DMDA().create([self.N_q1, self.N_q2], + dof = len(arrays), + proc_sizes = (PETSc.DECIDE, + PETSc.DECIDE + ), + comm = self._comm + ) + + self._glob_aux = self._da_aux_arrays.createGlobalVec() + self._glob_aux_array = self._glob_aux.getArray() + + self.dump_aux_arrays_initial_call = 0 + + N_g = self.N_ghost + + for i in range(len(arrays)): + if (i==0): + array_to_dump = arrays[0][:, N_g:-N_g, N_g:-N_g] + else: + array_to_dump = af.join(0, array_to_dump, + arrays[i][:, N_g:-N_g, N_g:-N_g] + ) + + af.flat(array_to_dump).to_ndarray(self._glob_aux_array) + PETSc.Object.setName(self._glob_aux, name) + viewer = PETSc.Viewer().createHDF5(file_name + '.h5', 'w', comm=self._comm) + viewer(self._glob_aux) + def dump_moments(self, file_name): """ This function is used to dump variables to a file for later usage. diff --git a/bolt/lib/nonlinear_solver/file_io/load.py b/bolt/lib/nonlinear_solver/file_io/load.py index ae433b2e..0b246993 100644 --- a/bolt/lib/nonlinear_solver/file_io/load.py +++ b/bolt/lib/nonlinear_solver/file_io/load.py @@ -24,7 +24,7 @@ def load_distribution_function(self, file_name): The above statemant will load the distribution function data stored in the file distribution_function.h5 into self.f """ - viewer = PETSc.Viewer().createHDF5(file_name + '.h5', + viewer = PETSc.Viewer().createHDF5(file_name, PETSc.Viewer.Mode.READ, comm=self._comm ) @@ -33,7 +33,7 @@ def load_distribution_function(self, file_name): N_g = self.N_ghost self.f[:, N_g:-N_g, N_g:-N_g] = af.moddims(af.to_array(self._glob_f_array), self.N_p1 * self.N_p2 * self.N_p3, - self.N_q1, self.N_q2 + self.N_q1_local, self.N_q2_local ) return diff --git a/bolt/lib/nonlinear_solver/interpolation_routines.py b/bolt/lib/nonlinear_solver/interpolation_routines.py index f029ac1c..23d41706 100644 --- a/bolt/lib/nonlinear_solver/interpolation_routines.py +++ b/bolt/lib/nonlinear_solver/interpolation_routines.py @@ -2,7 +2,6 @@ # -*- coding: utf-8 -*- import arrayfire as af -import numpy as np def f_interp_2d(self, dt): diff --git a/bolt/lib/nonlinear_solver/nonlinear_solver.py b/bolt/lib/nonlinear_solver/nonlinear_solver.py index 226083d2..98b0e570 100644 --- a/bolt/lib/nonlinear_solver/nonlinear_solver.py +++ b/bolt/lib/nonlinear_solver/nonlinear_solver.py @@ -53,7 +53,7 @@ from .utils.print_with_indent import indent from .utils.performance_timings import print_table from .compute_moments import compute_moments as compute_moments_imported -from .EM_fields_solver.electrostatic import fft_poisson +from .EM_fields_solver.electrostatic import fft_poisson, poisson_eqn_3D class nonlinear_solver(object): """ @@ -208,21 +208,53 @@ def __init__(self, physical_system, performance_test_flag = False): comm = self._comm ) - # Additionally, a DMDA object also needs to be created for - # the KSP/SNES solver with a DOF of 1. This is used to solve for - # the electrostatic case: - - self._da_ksp = PETSc.DMDA().create([self.N_q1, self.N_q2], - stencil_width = self.N_ghost, - boundary_type = (petsc_bc_in_q1, - petsc_bc_in_q2 - ), - proc_sizes = (PETSc.DECIDE, - PETSc.DECIDE - ), - stencil_type = 1, - comm = self._comm - ) + # Additionally, a DA object also needs to be created for the SNES solver + # with a DOF of 1: + # TODO: Remove the following hardcoded values + self.length_multiples_q1 = .5 + self.length_multiples_q2 = .25 + self.dq3 = self.dq1 + self.location_in_q3 = 0.3 + self.q3_3D_start = 0.; self.q3_3D_end = 1.3 + + self.N_q1_poisson = (2*self.length_multiples_q1+1)*self.N_q1 + self.N_q2_poisson = (2*self.length_multiples_q2+1)*self.N_q2 + self.N_q3_poisson = (int)((self.q3_3D_end - self.q3_3D_start) / self.dq3) + self.N_ghost_poisson = self.N_ghost + + PETSc.Sys.Print("dq3 = ", self.dq3, "N_q3 = ", self.N_q3_poisson) + self._da_snes = PETSc.DMDA().create([self.N_q1_poisson, + self.N_q2_poisson, + self.N_q3_poisson], + stencil_width = self.N_ghost_poisson, + boundary_type = (petsc_bc_in_q1, + petsc_bc_in_q2, + 'periodic' + ), + proc_sizes = (PETSc.DECIDE, + PETSc.DECIDE, + 1 + ), + stencil_type = 0, # Star stencil + dof = 1, + comm = self._comm + ) + self.snes = PETSc.SNES().create() + self.poisson = poisson_eqn_3D(self) + self.snes.setFunction(self.poisson.compute_residual, + self.poisson.glob_residual + ) + + self.snes.setDM(self._da_snes) + self.snes.setFromOptions() + + # Obtaining the left-bottom corner coordinates + # (lowest values of the canonical coordinates in the local zone) + ((i_q1_start, i_q2_start), (N_q1_local, N_q2_local)) = self._da_f.getCorners() + self.i_q1_start = i_q1_start + self.i_q2_start = i_q2_start + self.N_q1_local = N_q1_local + self.N_q2_local = N_q2_local # This DA is used by the FileIO routine dump_moments(): self._da_dump_moments = PETSc.DMDA().create([self.N_q1, self.N_q2], @@ -235,6 +267,8 @@ def __init__(self, physical_system, performance_test_flag = False): ), comm = self._comm ) + # For dumping aux arrays: + self.dump_aux_arrays_initial_call = 1 # Creation of the local and global vectors from the DA: # This is for the distribution function @@ -245,7 +279,7 @@ def __init__(self, physical_system, performance_test_flag = False): # the communication routines for EM fields self._glob_fields = self._da_fields.createGlobalVec() self._local_fields = self._da_fields.createLocalVec() - + # The following vector is used to dump the data to file: self._glob_moments = self._da_dump_moments.createGlobalVec() @@ -451,7 +485,6 @@ def _calculate_p_center(self): p1_center, p3_center ) - # Flattening the arrays: p1_center = af.flat(af.to_array(p1_center)) p2_center = af.flat(af.to_array(p2_center)) @@ -586,15 +619,16 @@ def _initialize(self, params): _apply_bcs_f = apply_boundary_conditions.apply_bcs_f _apply_bcs_fields = apply_boundary_conditions.apply_bcs_fields - strang_timestep = timestep.strang_step - lie_timestep = timestep.lie_step - swss_timestep = timestep.swss_step - jia_timestep = timestep.jia_step + strang_timestep = timestep.strang_step + lie_timestep = timestep.lie_step + swss_timestep = timestep.swss_step + jia_timestep = timestep.jia_step compute_moments = compute_moments_imported dump_distribution_function = dump.dump_distribution_function dump_moments = dump.dump_moments + dump_aux_arrays = dump.dump_aux_arrays load_distribution_function = load.load_distribution_function print_performance_timings = print_table diff --git a/bolt/lib/nonlinear_solver/tests/test_3D_poisson_solver_2D_density.py b/bolt/lib/nonlinear_solver/tests/test_3D_poisson_solver_2D_density.py new file mode 100644 index 00000000..07decc7c --- /dev/null +++ b/bolt/lib/nonlinear_solver/tests/test_3D_poisson_solver_2D_density.py @@ -0,0 +1,372 @@ + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import petsc4py, sys +petsc4py.init(sys.argv) +from petsc4py import PETSc +import arrayfire as af +import numpy as np +import pylab as pl + +pl.rcParams['figure.figsize'] = 20, 7.5 +pl.rcParams['figure.dpi'] = 150 +pl.rcParams['image.cmap'] = 'jet' +pl.rcParams['lines.linewidth'] = 1.5 +pl.rcParams['font.family'] = 'serif' +pl.rcParams['font.weight'] = 'bold' +pl.rcParams['font.size'] = 20 +pl.rcParams['font.sans-serif'] = 'serif' +pl.rcParams['text.usetex'] = True +pl.rcParams['axes.linewidth'] = 1.5 +pl.rcParams['axes.titlesize'] = 'medium' +pl.rcParams['axes.labelsize'] = 'medium' + +pl.rcParams['xtick.major.size'] = 8 +pl.rcParams['xtick.minor.size'] = 4 +pl.rcParams['xtick.major.pad'] = 8 +pl.rcParams['xtick.minor.pad'] = 8 +pl.rcParams['xtick.color'] = 'k' +pl.rcParams['xtick.labelsize'] = 'medium' +pl.rcParams['xtick.direction'] = 'in' + +pl.rcParams['ytick.major.size'] = 8 +pl.rcParams['ytick.minor.size'] = 4 +pl.rcParams['ytick.major.pad'] = 8 +pl.rcParams['ytick.minor.pad'] = 8 +pl.rcParams['ytick.color'] = 'k' +pl.rcParams['ytick.labelsize'] = 'medium' +pl.rcParams['ytick.direction'] = 'in' + +comm = PETSc.COMM_WORLD.tompi4py() + +N_q1_density = 33 +N_q2_density = 33 + +N_ghost = 1 + +q1_2D_start = -.5; q1_2D_end = 0.5 +q2_2D_start = -.5; q2_2D_end = 0.5 +location_in_q3 = 1. + +q3_3D_start = 0.; q3_3D_end = 2. + +da_2D = PETSc.DMDA().create([N_q1_density, + N_q2_density], + stencil_width = N_ghost, + boundary_type = ('periodic', + 'periodic' + ), + stencil_type = 1, + dof = 1, + comm = comm + ) + +glob_density = da_2D.createGlobalVec() +local_density = da_2D.createLocalVec() + + +((i_q1_2D_start, i_q2_2D_start), + (N_q1_2D_local, N_q2_2D_local) +) = \ + da_2D.getCorners() + +dq1_2D = (q1_2D_end - q1_2D_start) / N_q1_density +dq2_2D = (q2_2D_end - q2_2D_start) / N_q2_density + +i_q1_2D = ( (i_q1_2D_start + 0.5) + + np.arange(-N_ghost, N_q1_2D_local + N_ghost) + ) + +i_q2_2D = ( (i_q2_2D_start + 0.5) + + np.arange(-N_ghost, N_q2_2D_local + N_ghost) + ) + +q1_2D = q1_2D_start + i_q1_2D * dq1_2D +q2_2D = q2_2D_start + i_q2_2D * dq2_2D + +dq1_3D = dq1_2D +dq2_3D = dq2_2D +dq3_3D = dq1_2D + +length_multiples_q1 = 1 +length_multiples_q2 = 1 +N_q1_poisson = (2*length_multiples_q1+1)*N_q1_density +N_q2_poisson = (2*length_multiples_q2+1)*N_q2_density +N_q3_poisson = (int)((q3_3D_end - q3_3D_start) / dq1_3D) + +da_3D = PETSc.DMDA().create([N_q1_poisson, + N_q2_poisson, + N_q3_poisson], + stencil_width = N_ghost, + boundary_type = ('periodic', + 'periodic', + 'periodic' + ), + stencil_type = 1, + dof = 1, + comm = comm + ) + +glob_phi = da_3D.createGlobalVec() +local_phi = da_3D.createLocalVec() +glob_residual = da_3D.createGlobalVec() + +glob_epsilon = da_3D.createGlobalVec() +local_epsilon = da_3D.createLocalVec() + +((i_q1_3D_start, i_q2_3D_start, i_q3_3D_start), + (N_q1_3D_local, N_q2_3D_local, N_q3_3D_local) +) = \ + da_3D.getCorners() + +i_q1_3D = ( (i_q1_3D_start + 0.5) + + np.arange(-N_ghost, N_q1_3D_local + N_ghost) + ) + +i_q2_3D = ( (i_q2_3D_start + 0.5) + + np.arange(-N_ghost, N_q2_3D_local + N_ghost) + ) + +i_q3_3D = ( (i_q3_3D_start + 0.5) + + np.arange(-N_ghost, N_q3_3D_local + N_ghost) + ) + +length_q1_2d = (q1_2D_end - q1_2D_start) +length_q2_2d = (q2_2D_end - q2_2D_start) + +q1_3D = q1_2D_start - length_multiples_q1*length_q1_2d + i_q1_3D * dq1_3D +q2_3D = q2_2D_start - length_multiples_q2*length_q2_2d + i_q2_3D * dq2_3D +q3_3D = q3_3D_start + i_q3_3D * dq3_3D + +glob_density_array = glob_density.getArray(readonly=0) +glob_density_array = glob_density_array.reshape([N_q2_2D_local, \ + N_q1_2D_local], \ + ) +glob_density_array[:] = -335. + +epsilon_array = local_epsilon.getArray(readonly=0) +epsilon_array = epsilon_array.reshape([N_q3_3D_local + 2*N_ghost, \ + N_q2_3D_local + 2*N_ghost, \ + N_q1_3D_local + 2*N_ghost + ] + ) +epsilon_array[:] = 1. + + +print("rank = ", comm.rank) +print("N_q1_poisson = ", N_q1_poisson) +print("N_q2_poisson = ", N_q2_poisson) +print("N_q3_poisson = ", N_q3_poisson) + +offset = 2*N_ghost-1 +print("i_q1_local_3D = [", i_q1_3D_start, ",", i_q1_3D_start+N_q1_3D_local+offset, "]") +print("i_q2_local_3D = [", i_q2_3D_start, ",", i_q2_3D_start+N_q2_3D_local+offset, "]") +print("i_q3_local_3D = [", i_q3_3D_start, ",", i_q3_3D_start+N_q3_3D_local+offset, "]") +print(" ") +print("i_q1_local_2D = [", i_q1_2D_start, ",", i_q1_2D_start+N_q2_2D_local+offset, "]") +print("i_q2_local_2D = [", i_q2_2D_start, ",", i_q2_2D_start+N_q2_2D_local+offset, "]") + +# Figure out the coordinates of the 3D phi cube of the current rank +print("q1_3D_start = ", q1_3D[N_ghost]) +print("q2_3D_start = ", q2_3D[N_ghost]) +print("q3_3D_start = ", q3_3D[N_ghost]) +print(" ") +print("q1_2D_start = ", q1_2D[N_ghost]) +print("q2_2D_start = ", q2_2D[N_ghost]) + +q3_2D_in_3D_index_start = np.where(q3_3D > location_in_q3 - dq3_3D)[0][0] +q3_2D_in_3D_index_end = np.where(q3_3D < location_in_q3 + dq3_3D)[0][-1] + +q1_2D_in_3D_index_start = np.where(abs(q1_3D - q1_2D[0+N_ghost] ) < 1e-10)[0][0] +q1_2D_in_3D_index_end = np.where(abs(q1_3D - q1_2D[-1-N_ghost]) < 1e-10)[0][0] +q2_2D_in_3D_index_start = np.where(abs(q2_3D - q2_2D[0+N_ghost] ) < 1e-10)[0][0] +q2_2D_in_3D_index_end = np.where(abs(q2_3D - q2_2D[-1-N_ghost]) < 1e-10)[0][0] + +print("q1_2D_in_3D_index_start = ", q1_2D_in_3D_index_start, "q1_3D_start = ", q1_3D[q1_2D_in_3D_index_start]) +print("q1_2D_in_3D_index_end = ", q1_2D_in_3D_index_end, "q1_3D_end = ", q1_3D[q1_2D_in_3D_index_end]) +print("q2_2D_in_3D_index_start = ", q2_2D_in_3D_index_start, "q2_3D_start = ", q2_3D[q2_2D_in_3D_index_start]) +print("q2_2D_in_3D_index_end = ", q2_2D_in_3D_index_end, "q2_3D_end = ", q2_3D[q2_2D_in_3D_index_end]) +print("q3_2D_in_3D_index_start = ", q3_2D_in_3D_index_start, "q3_3D_start = ", q3_3D[q3_2D_in_3D_index_start]) +print("q3_2D_in_3D_index_end = ", q3_2D_in_3D_index_end, "q3_3D_end = ", q3_3D[q3_2D_in_3D_index_end]) + +print(" ") +print("After ghost zone offset:") +print("q1_3D_start = ", q1_3D[q1_2D_in_3D_index_start]) +print("q1_3D_end = ", q1_3D[q1_2D_in_3D_index_end]) +print("q2_3D_start = ", q2_3D[q2_2D_in_3D_index_start]) +print("q2_3D_end = ", q2_3D[q2_2D_in_3D_index_end]) + +epsilon_array[:q3_2D_in_3D_index_start, + :, : + ] = 10. +q3_3D_data_structure = 0.*epsilon_array +for j in range(q3_3D_data_structure.shape[1]): + for i in range(q3_3D_data_structure.shape[2]): + q3_3D_data_structure[:, j, i] = q3_3D + +print("z = ", q3_3D_data_structure[q3_2D_in_3D_index_start, 0, 0]) +class poisson_eqn(object): + + def __init__(self): + self.local_phi = local_phi + self.residual_counter = 0 + + def compute_residual(self, snes, phi, residual): + self.residual_counter += 1 +# print("residual iter = ", self.residual_counter) + da_3D.globalToLocal(phi, local_phi) + + N_g = N_ghost + + phi_array = local_phi.getArray(readonly=0) + phi_array = phi_array.reshape([N_q3_3D_local + 2*N_g, \ + N_q2_3D_local + 2*N_g, \ + N_q1_3D_local + 2*N_g + ] + ) + + residual_array = residual.getArray(readonly=0) + residual_array = residual_array.reshape([N_q3_3D_local, \ + N_q2_3D_local, \ + N_q1_3D_local + ] + ) + + phi_array[:N_ghost, :, :] = 0. + phi_array[N_q3_3D_local+N_ghost:, :, :] = 0. + phi_array[:, :N_ghost, :] = 0. + phi_array[:, N_q2_3D_local+N_ghost:, :] = 0. + phi_array[:, :, :N_ghost] = 0. + phi_array[:, :, N_q1_3D_local+N_ghost:] = 0. + + z = q3_3D_data_structure + z_sample = q3_3D[q3_2D_in_3D_index_start] + z_backgate = q3_3D[0] + side_wall_boundaries = (z_sample - z)/(z_sample - z_backgate) + + phi_array[:q3_2D_in_3D_index_start, :N_ghost, :] = \ + side_wall_boundaries[:q3_2D_in_3D_index_start, :N_ghost, :] + + phi_array[:q3_2D_in_3D_index_start, N_q2_3D_local+N_ghost:, :] = \ + side_wall_boundaries[:q3_2D_in_3D_index_start, N_q2_3D_local+N_ghost:, :] + + phi_array[:q3_2D_in_3D_index_start, :, :N_ghost] = \ + side_wall_boundaries[:q3_2D_in_3D_index_start, :, :N_ghost] + + phi_array[:q3_2D_in_3D_index_start, :, N_q1_3D_local+N_ghost:] = \ + side_wall_boundaries[:q3_2D_in_3D_index_start, :, N_q1_3D_local+N_ghost:] + + #Backgate +# phi_array[:N_ghost, +# q2_2D_in_3D_index_start:q2_2D_in_3D_index_end+1, +# q1_2D_in_3D_index_start:q1_2D_in_3D_index_end+1 +# ] = 1. + phi_array[:N_ghost, :, :] = 1. + + phi_plus_x = np.roll(phi_array, shift=-1, axis=2) # (i+3/2, j+1/2, k+1/2) + phi_minus_x = np.roll(phi_array, shift=1, axis=2) # (i-1/2, j+1/2, k+1/2) + phi_plus_y = np.roll(phi_array, shift=-1, axis=1) # (i+1/2, j+3/2, k+1/2) + phi_minus_y = np.roll(phi_array, shift=1, axis=1) # (i+1/2, j-1/2, k+1/2) + phi_plus_z = np.roll(phi_array, shift=-1, axis=0) # (i+1/2, j+1/2, k+3/2) + phi_minus_z = np.roll(phi_array, shift=1, axis=0) # (i+1/2, j+1/2, k+3/2) + + eps_left_edge = epsilon_array # (i, j+1/2, k+1/2) + eps_right_edge = np.roll(epsilon_array, shift=-1, axis=2) # (i+1, j+1/2, k+1/2) + + eps_bot_edge = epsilon_array # (i+1/2, j, k+1/2) + eps_top_edge = np.roll(epsilon_array, shift=-1, axis=1) # (i+1/2, j+1, k+1/2) + + eps_back_edge = epsilon_array # (i+1/2, j+1/2, k) + eps_front_edge = np.roll(epsilon_array, shift=-1, axis=0) # (i+1/2, j+1/2, k+1) + + D_left_edge = eps_left_edge * (phi_array - phi_minus_x)/dq1_3D + D_right_edge = eps_right_edge * (phi_plus_x - phi_array) /dq1_3D + + D_bot_edge = eps_bot_edge * (phi_array - phi_minus_y)/dq2_3D + D_top_edge = eps_top_edge * (phi_plus_y - phi_array )/dq2_3D + + D_back_edge = eps_back_edge * (phi_array - phi_minus_z)/dq3_3D + D_front_edge = eps_front_edge * (phi_plus_z - phi_array )/dq3_3D + + laplacian_phi = (D_right_edge - D_left_edge) /dq1_3D \ + + (D_top_edge - D_bot_edge) /dq2_3D \ + + (D_front_edge - D_back_edge) /dq3_3D + +# d2phi_dx2 = (phi_minus_x - 2.*phi_array + phi_plus_x)/dq1_3D**2. +# d2phi_dy2 = (phi_minus_y - 2.*phi_array + phi_plus_y)/dq2_3D**2. +# d2phi_dz2 = (phi_minus_z - 2.*phi_array + phi_plus_z)/dq3_3D**2. +# +# laplacian_phi = d2phi_dx2 + d2phi_dy2 + d2phi_dz2 + + laplacian_phi[q3_2D_in_3D_index_start, + q2_2D_in_3D_index_start:q2_2D_in_3D_index_end+1, + q1_2D_in_3D_index_start:q1_2D_in_3D_index_end+1 + ] \ + += glob_density_array + + residual_array[:, :, :] = \ + laplacian_phi[N_g:-N_g, N_g:-N_g, N_g:-N_g] + + #Side contacts + mid_point_q2_index = \ + (int)((q2_2D_in_3D_index_start + q2_2D_in_3D_index_end)/2) + + residual_array[q3_2D_in_3D_index_start-N_g, + mid_point_q2_index-5-N_g:mid_point_q2_index+5+1-N_g, + :q1_2D_in_3D_index_start-N_g + ] = \ + phi_array[q3_2D_in_3D_index_start, + mid_point_q2_index-5:mid_point_q2_index+5+1, + N_g:q1_2D_in_3D_index_start + ] - 0.1 + + residual_array[q3_2D_in_3D_index_start-N_g, + mid_point_q2_index-5-N_g:mid_point_q2_index+5+1-N_g, + q1_2D_in_3D_index_end+1-N_g: + ] = \ + phi_array[q3_2D_in_3D_index_start, + mid_point_q2_index-5:mid_point_q2_index+5+1, + q1_2D_in_3D_index_end+1:-N_g + ] + 0.1*0. + + return + +snes = PETSc.SNES().create() +pde = poisson_eqn() +snes.setFunction(pde.compute_residual, glob_residual) + +snes.setDM(da_3D) +snes.setFromOptions() + +snes.solve(None, glob_phi) +phi_array = glob_phi.getArray() +phi_array = phi_array.reshape([N_q3_3D_local, \ + N_q2_3D_local, \ + N_q1_3D_local] + ) +pl.subplot(121) +pl.contourf( + phi_array[q3_2D_in_3D_index_start, :, :], 100, cmap='jet' + ) +pl.colorbar() +pl.title('Top View') +pl.xlabel('$x$') +pl.ylabel('$y$') +pl.gca().set_aspect('equal') + +pl.subplot(122) +pl.contourf(phi_array[:, N_q2_poisson/2, :], 100, cmap='jet') +pl.title('Side View') +pl.xlabel('$x$') +pl.ylabel('$z$') +pl.colorbar() +pl.gca().set_aspect('equal') +pl.show() + +#for n in range(10): +# +# print("====== n = ", n, "======") +# glob_density_array[:] = n +# snes.solve(None, glob_phi) +# pde.residual_counter = 0 diff --git a/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py b/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py index 650f5ca9..cd3e74ea 100644 --- a/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py +++ b/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py @@ -1,123 +1,194 @@ -# #!/usr/bin/env python3 -# # -*- coding: utf-8 -*- +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- -##NOT WORKING IN THIS BRANCH +""" +In this test we check that the 2D Poisson solver +works as intended. For this purpose, we assign +a density distribution for which the analytical +solution for electrostatic fields may be computed. -# """ -# In this test we check that the 2D Poisson solver -# works as intended. For this purpose, we assign -# a density distribution for which the analytical -# solution for electrostatic fields may be computed. +This solution is then checked against the solution +given by the KSP solver +""" -# This solution is then checked against the solution -# given by the KSP solver -# """ +import numpy as np +import arrayfire as af +import sys, petsc4py; petsc4py.init(sys.argv) +from petsc4py import PETSc +import pylab as pl -# import numpy as np -# import arrayfire as af -# import sys, petsc4py; petsc4py.init(sys.argv) -# from petsc4py import PETSc -# from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic \ -# import compute_electrostatic_fields +from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic \ + import compute_electrostatic_fields -# def compute_moments_sinusoidal(self, *args): -# return (1 + af.sin(2 * np.pi * self.q1 + 4 * np.pi * self.q2)) +def compute_moments_sinusoidal(self, *args): + rho = 1. + 0.1*af.sin(2 * np.pi * self.q1 + 4 * np.pi * self.q2) -# def compute_moments_gaussian(self, *args): -# q2_minus = 0.25 -# q2_plus = 0.75 + rho_mean = af.mean(rho) + print("rho_mean = ", rho_mean) -# regulator = 40 # larger value makes the transition sharper + return(rho - rho_mean) -# rho = 1 + 0.5 * ( af.tanh(( self.q2 - q2_minus)*regulator) -# - af.tanh(( self.q2 - q2_plus )*regulator) -# ) +def compute_moments_gaussian(self, *args): + q2_minus = 0.25 + q2_plus = 0.75 -# return(rho) + regulator = 40 # larger value makes the transition sharper -# class test(object): -# def __init__(self, N): -# # Creating an object with necessary attributes: -# self.physical_system = type('obj', (object, ), -# {'params': type('obj', (object, ), -# {'charge_electron': -1}) -# } -# ) - -# self.q1_start = 0 -# self.q2_start = 0 - -# self.q1_end = 1 -# self.q2_end = 1 + rho = 1 + 0.5 * ( af.tanh(( self.q2 - q2_minus)*regulator) + - af.tanh(( self.q2 - q2_plus )*regulator) + ) -# self.N_q1 = N -# self.N_q2 = N +# rho = af.exp(-0.*(self.q1 - 0.5)**2./0.01 - (self.q2 - 0.5)**2./0.01) +# rho = (self.q2 - 0.5)**2. + (self.q1 - 0.5)**2. -# self.dq1 = (self.q1_end - self.q1_start) / self.N_q1 -# self.dq2 = (self.q2_end - self.q2_start) / self.N_q2 +# rho = 1. + 0.1*af.sin(2*np.pi*self.q2) -# self.N_ghost = np.random.randint(3, 5) +# sigma = 0.1 +# rho = af.exp(-(self.q1)**2./(2.*sigma**2.) -(self.q2)**2./(2.*sigma**2.)) \ +# * 1./ sigma**2. / (2. * np.pi) -# self.q1 = self.q1_start \ -# + (0.5 + np.arange(-self.N_ghost, -# self.N_q1 + self.N_ghost -# ) -# ) * self.dq1 - -# self.q2 = self.q2_start \ -# + (0.5 + np.arange(-self.N_ghost, -# self.N_q2 + self.N_ghost -# ) -# ) * self.dq2 - -# self.q2, self.q1 = np.meshgrid(self.q2, self.q1) -# self.q2, self.q1 = af.to_array(self.q2), af.to_array(self.q1) + N_g = self.N_ghost + net_charge = af.sum(rho[N_g:-N_g, N_g:-N_g]) * self.dq1 * self.dq2 -# self.E1 = af.constant(0, self.q1.shape[0], self.q1.shape[1], -# dtype=af.Dtype.f64 -# ) + total_volume = (self.q1_end - self.q1_start) \ + * (self.q2_end - self.q2_start) -# self.E2 = self.E1.copy() -# self.E3 = self.E1.copy() + rho_zero_net_charge = rho - net_charge/total_volume -# self.B1 = self.E1.copy() -# self.B2 = self.E1.copy() -# self.B3 = self.E1.copy() + print("Initial net charge = ", net_charge) -# self._comm = PETSc.COMM_WORLD.tompi4py() + return(rho_zero_net_charge) -# self._da_ksp = PETSc.DMDA().create([self.N_q1, self.N_q2], -# stencil_width=self.N_ghost, -# boundary_type=('periodic', -# 'periodic'), -# stencil_type=1, -# ) +class test(object): + def __init__(self, N_q1, N_q2): + # Creating an object with necessary attributes: + self.physical_system = type('obj', (object, ), + {'params': type('obj', (object, ), + {'charge_electron': -1}) + } + ) -# self.performance_test_flag = False + self.q1_start = 0. + self.q2_start = 0. -# compute_moments = compute_moments_sinusoidal + self.q1_end = 1. + self.q2_end = 1. -# def test_compute_electrostatic_fields_1(): + self.N_q1 = N_q1 + self.N_q2 = N_q2 -# error_E1 = np.zeros(5) -# error_E2 = np.zeros(5) + self.dq1 = (self.q1_end - self.q1_start) / self.N_q1 + self.dq2 = (self.q2_end - self.q2_start) / self.N_q2 -# N = 2**np.arange(5, 10) + #self.N_ghost = np.random.randint(3, 5) + self.N_ghost = 3 -# for i in range(N.size): -# obj = test(N[i]) -# compute_electrostatic_fields(obj) - -# E1_expected = (0.1 / np.pi) \ -# * af.cos( 2 * np.pi * obj.q1 -# + 4 * np.pi * obj.q2 -# ) - -# E2_expected = (0.2 / np.pi) \ + self.q1 = self.q1_start \ + + (0.5 + np.arange(-self.N_ghost, + self.N_q1 + self.N_ghost + ) + ) * self.dq1 + + self.q2 = self.q2_start \ + + (0.5 + np.arange(-self.N_ghost, + self.N_q2 + self.N_ghost + ) + ) * self.dq2 + + self.q2, self.q1 = np.meshgrid(self.q2, self.q1) + self.q2, self.q1 = af.to_array(self.q2), af.to_array(self.q1) + + self.E1 = af.constant(0, self.q1.shape[0], self.q1.shape[1], + dtype=af.Dtype.f64 + ) + + self.E2 = self.E1.copy() + self.E3 = self.E1.copy() + + self.B1 = self.E1.copy() + self.B2 = self.E1.copy() + self.B3 = self.E1.copy() + + self._comm = PETSc.COMM_WORLD.tompi4py() + + self._da_snes = PETSc.DMDA().create([self.N_q1, self.N_q2], + stencil_width=self.N_ghost, + boundary_type=('periodic', + 'periodic'), + stencil_type=1, + dof = 1 + ) + self.glob_residual = self._da_snes.createGlobalVec() + self.glob_phi = self._da_snes.createGlobalVec() + self.glob_phi.set(0.) + + self.performance_test_flag = False + + compute_moments = compute_moments_sinusoidal + +def test_compute_electrostatic_fields_1(): + obj = test(70, 70) + compute_electrostatic_fields(obj) + +# N = 7*np.array([2, 4, 6, 8, 10, 12]) +# error_E1 = np.zeros(N.size) +# error_E2 = np.zeros(N.size) +# +# for i in range(N.size): +# obj = test(N[i], N[i]) +# compute_electrostatic_fields(obj) +# +# E1_expected = (0.1 / np.pi) \ # * af.cos( 2 * np.pi * obj.q1 # + 4 * np.pi * obj.q2 # ) +# +# E2_expected = (0.2 / np.pi) \ +# * af.cos( 2 * np.pi * obj.q1 +# + 4 * np.pi * obj.q2 +# ) +# +# N_g = obj.N_ghost +# +# error_E1[i] = af.sum(af.abs( obj.E1[N_g:-N_g, N_g:-N_g] +# - E1_expected[N_g:-N_g, N_g:-N_g] +# ) +# ) / (obj.E1[N_g:-N_g, N_g:-N_g].elements()) +# +# error_E2[i] = af.sum(af.abs( obj.E2[N_g:-N_g, N_g:-N_g] +# - E2_expected[N_g:-N_g, N_g:-N_g] +# ) +# ) / (obj.E2[N_g:-N_g, N_g:-N_g].elements()) +# +# print("Error E1 = ", error_E1) +# print("Error E2 = ", error_E2) +# +# poly_E1 = np.polyfit(np.log10(N), np.log10(error_E1), 1) +# poly_E2 = np.polyfit(np.log10(N), np.log10(error_E2), 1) +# +# assert (abs(poly_E1[0] + 2) < 0.2) +# assert (abs(poly_E2[0] + 2) < 0.2) + +# def test_compute_electrostatic_fields_2(): + +# N = 7*np.array([2, 4, 6, 8, 10, 12]) +# N = 7*np.array([12]) +# error_E1 = np.zeros(N.size) +# error_E2 = np.zeros(N.size) +# +# for i in range(N.size): +# obj = test(N[i], N[i]) +# compute_electrostatic_fields(obj) +# +# E1_expected = 0 * obj.q1 +# +# q2_minus = 0.25 +# q2_plus = 0.75 + +# E2_expected = -0.5/40 * ( af.log(af.cosh(( obj.q2 - 0.25)*40)) +# - af.log(af.cosh(( obj.q2 - 0.75)*40)) +# ) # N_g = obj.N_ghost @@ -131,45 +202,14 @@ # ) # ) / (obj.E2[N_g:-N_g, N_g:-N_g].elements()) -# poly_E1 = np.polyfit(np.log10(N), np.log10(error_E1), 1) -# poly_E2 = np.polyfit(np.log10(N), np.log10(error_E2), 1) - -# assert (abs(poly_E1[0] + 2) < 0.2) -# assert (abs(poly_E2[0] + 2) < 0.2) - -# # def test_compute_electrostatic_fields_2(): - -# # error_E1 = np.zeros(5) -# # error_E2 = np.zeros(5) - -# # N = 2**np.arange(10, 11) - -# # for i in range(N.size): -# # obj = test(N[i]) -# # compute_electrostatic_fields(obj) - -# # E1_expected = 0 - -# # E2_expected = -0.5/40 * ( af.log(af.cosh(( obj.q2 - 0.25)*40)) -# # - af.log(af.cosh(( obj.q2 - 0.75)*40)) -# # ) - -# # N_g = obj.N_ghost - -# # error_E1[i] = af.sum(af.abs( obj.E1[N_g:-N_g, N_g:-N_g] -# # - E1_expected[N_g:-N_g, N_g:-N_g] -# # ) -# # ) / (obj.E1[N_g:-N_g, N_g:-N_g].elements()) - -# # error_E2[i] = af.sum(af.abs( obj.E2[N_g:-N_g, N_g:-N_g] -# # - E2_expected[N_g:-N_g, N_g:-N_g] -# # ) -# # ) / (obj.E2[N_g:-N_g, N_g:-N_g].elements()) - -# # poly_E1 = np.polyfit(np.log10(N), np.log10(error_E1), 1) -# # poly_E2 = np.polyfit(np.log10(N), np.log10(error_E2), 1) +# print("Error E1 = ", error_E1) +# print("Error E2 = ", error_E2) +# +# poly_E1 = np.polyfit(np.log10(N), np.log10(error_E1), 1) +# poly_E2 = np.polyfit(np.log10(N), np.log10(error_E2), 1) +# +# assert (abs(poly_E1[0] + 2) < 0.2) +# assert (abs(poly_E2[0] + 2) < 0.2) -# # assert (abs(poly_E1[0] + 2) < 0.2) -# # assert (abs(poly_E2[0] + 2) < 0.2) -# test_compute_electrostatic_fields_1() +test_compute_electrostatic_fields_1() diff --git a/bolt/lib/nonlinear_solver/timestep.py b/bolt/lib/nonlinear_solver/timestep.py index f624eac7..7321b5a4 100644 --- a/bolt/lib/nonlinear_solver/timestep.py +++ b/bolt/lib/nonlinear_solver/timestep.py @@ -18,20 +18,7 @@ # Defining the operators: # When using FVM: def op_fvm_q(self, dt): - - self._communicate_f() - self._apply_bcs_f() - - if(self.performance_test_flag == True): - tic = af.time() - - fvm_timestep_RK2(self, dt) - if(self.performance_test_flag == True): - af.sync() - toc = af.time() - self.time_fvm_solver += toc - tic - # Solving for tau = 0 systems if(af.any_true(self.physical_system.params.tau(self.q1_center, self.q2_center, self.p1, self.p2, self.p3 @@ -47,12 +34,29 @@ def op_fvm_q(self, dt): self.physical_system.params, True ) + + fvm_timestep_RK2(self, dt) if(self.performance_test_flag == True): af.sync() toc = af.time() self.time_sourcets += toc - tic + af.eval(self.f) + return + + + if(self.performance_test_flag == True): + tic = af.time() + + fvm_timestep_RK2(self, dt) + + if(self.performance_test_flag == True): + af.sync() + toc = af.time() + self.time_fvm_solver += toc - tic + + af.eval(self.f) return @@ -91,7 +95,7 @@ def op_solve_src(self, dt): self.compute_moments, self.physical_system.params ) - + if(self.performance_test_flag == True): af.sync() toc = af.time() diff --git a/bolt/lib/physical_system.py b/bolt/lib/physical_system.py index 2a63335e..40817d25 100644 --- a/bolt/lib/physical_system.py +++ b/bolt/lib/physical_system.py @@ -98,21 +98,21 @@ def __init__(self, raise TypeError('Expected source_or_sink to be of type function') # Checking for the types of the methods in advection_term: - attributes = [a for a in dir(advection_term) if not a.startswith('_')] - for i in range(len(attributes)): - if(isinstance(getattr(advection_term, attributes[i]), - types.FunctionType - ) is False - ): - raise TypeError('Expected attributes of advection_term \ - to be of type function' - ) +# attributes = [a for a in dir(advection_term) if not a.startswith('_')] +# for i in range(len(attributes)): +# if(isinstance(getattr(advection_term, attributes[i]), +# types.FunctionType +# ) is False +# ): +# raise TypeError('Expected attributes of advection_term \ +# to be of type function' +# ) # Checking for the data-type in moment_defs: if(not isinstance(moment_defs.moment_exponents, dict) or not isinstance(moment_defs.moment_coeffs, dict) ): - raise TypeError('Expected attributes of boundary_conditions \ + raise TypeError('Expected attributes of moment definitions \ to be of type dict' ) diff --git a/bolt/src/electronic_boltzmann/advection_terms.py b/bolt/src/electronic_boltzmann/advection_terms.py new file mode 100644 index 00000000..694512d0 --- /dev/null +++ b/bolt/src/electronic_boltzmann/advection_terms.py @@ -0,0 +1,37 @@ +#import numpy as np +import arrayfire as af + +#@af.broadcast +def A_q(q1, q2, p1, p2, p3, params): + """Return the terms A_q1, A_q2.""" + + A_q1, A_q2 = params.vel_band + + return (A_q1, A_q2) + +def C_q(q1, q2, p1, p2, p3, params): + """Return the terms A_q1, A_q2.""" + + C_q1, C_q2 = params.vel_band + + return (C_q1, C_q2) + +# This can then be called inside A_p if needed: +# F1 = (params.char....)(E1 + ....) + T1(q1, q2, p1, p2, p3) + +def A_p(q1, q2, p1, p2, p3, + E1, E2, E3, B1, B2, B3, + params + ): + """Return the terms A_p1, A_p2 and A_p3.""" + e = params.charge_electron + c = params.speed_of_light + B3_mean = params.B3_mean + + v1, v2 = params.vel_band + + dp1_dt = -e*(E1 + v2*B3_mean/c) # p1 = hcross * k1 + dp2_dt = -e*(E2 - v1*B3_mean/c) # p2 = hcross * k2 + dp3_dt = 0.*p1 + + return (dp1_dt, dp2_dt, dp3_dt) diff --git a/bolt/src/electronic_boltzmann/collision_operator.py b/bolt/src/electronic_boltzmann/collision_operator.py new file mode 100644 index 00000000..36f50869 --- /dev/null +++ b/bolt/src/electronic_boltzmann/collision_operator.py @@ -0,0 +1,357 @@ +"""Contains the function which returns the Source/Sink term.""" + +from petsc4py import PETSc +import numpy as np +import arrayfire as af +from .matrix_inverse import inverse_4x4_matrix +import domain + +@af.broadcast +def f0_defect_constant_T(f, p1, p2, p3, params): + + mu = params.mu + T = params.T + + for n in range(params.collision_nonlinear_iters): + + E_upper = params.E_band + k = params.boltzmann_constant + + tmp = ((E_upper - mu)/(k*T)) + denominator = (k*T**2.*(af.exp(tmp) + 2. + af.exp(-tmp)) ) + + # TODO: Multiply with the integral measure dp1 * dp2 + a00 = af.sum(T / denominator, 0) + + fermi_dirac = 1./(af.exp( (E_upper - mu)/(k*T) ) + 1.) + af.eval(fermi_dirac) + + zeroth_moment = f - fermi_dirac + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + + N_g = domain.N_ghost + error_mass_conservation = af.max(af.abs(eqn_mass_conservation)[0, N_g:-N_g, N_g:-N_g]) + + print(" rank = ", params.rank, + "||residual_defect|| = ", error_mass_conservation + ) + + res = eqn_mass_conservation + dres_dmu = -a00 + + delta_mu = -res/dres_dmu + + mu = mu + delta_mu + + af.eval(mu) + + # Solved for mu. Now store in params + params.mu = mu + + # Print final residual + fermi_dirac = 1./(af.exp( (E_upper - mu)/(k*T) ) + 1.) + af.eval(fermi_dirac) + + zeroth_moment = f - fermi_dirac + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + + N_g = domain.N_ghost + error_mass_conservation = af.max(af.abs(eqn_mass_conservation)[0, N_g:-N_g, N_g:-N_g]) + + print(" rank = ", params.rank, + "||residual_defect|| = ", error_mass_conservation + ) + print(" rank = ", params.rank, + "mu = ", af.mean(params.mu[0, N_g:-N_g, N_g:-N_g]), + "T = ", af.mean(params.T[0, N_g:-N_g, N_g:-N_g]) + ) + PETSc.Sys.Print(" ------------------") + + return(fermi_dirac) + + +# Using af.broadcast, since p1, p2, p3 are of size (1, 1, Np1*Np2*Np3) +# All moment quantities are of shape (Nq1, Nq2) +# By wrapping with af.broadcast, we can perform batched operations +# on arrays of different sizes. +@af.broadcast +def f0_defect(f, p1, p2, p3, params): + + # Initial guess + mu = params.mu + T = params.T + + for n in range(params.collision_nonlinear_iters): + + E_upper = params.E_band + k = params.boltzmann_constant + + tmp = ((E_upper - mu)/(k*T)) + denominator = (k*T**2.*(af.exp(tmp) + 2. + af.exp(-tmp)) ) + + # TODO: Multiply with the integral measure dp1 * dp2 + a00 = af.sum(T / denominator, 0) + a01 = af.sum((E_upper - mu) / denominator, 0) + a10 = af.sum(E_upper*T / denominator, 0) + a11 = af.sum(E_upper*(E_upper - mu) / denominator, 0) + + # Solve Ax = b + # where A == Jacobian, + # x == delta guess (correction to guess), + # b = -residual + + fermi_dirac = 1./(af.exp( (E_upper - mu)/(k*T) ) + 1.) + af.eval(fermi_dirac) + + zeroth_moment = f - fermi_dirac + second_moment = E_upper*(f - fermi_dirac) + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + eqn_energy_conservation = af.sum(second_moment, 0) + + N_g = domain.N_ghost + error_mass_conservation = af.max(af.abs(eqn_mass_conservation)[0, N_g:-N_g, N_g:-N_g]) + error_energy_conservation = af.max(af.abs(eqn_energy_conservation)[0, N_g:-N_g, N_g:-N_g]) + + residual = [eqn_mass_conservation, eqn_energy_conservation] + error_norm = np.max([af.max(af.abs(residual[0])), + af.max(af.abs(residual[1]))] + ) + print(" ||residual_defect|| = ", error_norm) + +# if (error_norm < 1e-9): +# params.mu = mu +# params.T = T +# return(fermi_dirac) + + b0 = eqn_mass_conservation + b1 = eqn_energy_conservation + + det = a01*a10 - a00*a11 + delta_mu = -(a11*b0 - a01*b1)/det + delta_T = (a10*b0 - a00*b1)/det + + mu = mu + delta_mu + T = T + delta_T + + af.eval(mu, T) + + # Solved for (mu, T). Now store in params + params.mu = mu + params.T = T + + # Print final residual + fermi_dirac = 1./(af.exp( (E_upper - mu)/(k*T) ) + 1.) + af.eval(fermi_dirac) + + zeroth_moment = f - fermi_dirac + second_moment = E_upper*(f - fermi_dirac) + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + eqn_energy_conservation = af.sum(second_moment, 0) + + residual = [eqn_mass_conservation, eqn_energy_conservation] + error_norm = np.max([af.max(af.abs(residual[0])), + af.max(af.abs(residual[1]))] + ) + print(" ||residual_defect|| = ", error_norm) + print(" mu = ", af.mean(params.mu[0, N_g:-N_g, N_g:-N_g]), + "T = ", af.mean(params.T[0, N_g:-N_g, N_g:-N_g]) + ) + print(" ------------------") + + return(fermi_dirac) + +@af.broadcast +def f0_ee(f, p1, p2, p3, params): + + # Initial guess + mu_ee = params.mu_ee + T_ee = params.T_ee + vel_drift_x = params.vel_drift_x + vel_drift_y = params.vel_drift_y + + for n in range(params.collision_nonlinear_iters): + + E_upper = params.E_band + k = params.boltzmann_constant + + tmp1 = (E_upper - mu_ee - p1*vel_drift_x - p2*vel_drift_y) + tmp = (tmp1/(k*T_ee)) + denominator = (k*T_ee**2.*(af.exp(tmp) + 2. + af.exp(-tmp)) ) + + a_0 = T_ee / denominator + a_1 = tmp1 / denominator + a_2 = T_ee * p1 / denominator + a_3 = T_ee * p2 / denominator + + af.eval(a_0, a_1, a_2, a_3) + + # TODO: Multiply with the integral measure dp1 * dp2 + a_00 = af.sum(a_0, 0) + a_01 = af.sum(a_1, 0) + a_02 = af.sum(a_2, 0) + a_03 = af.sum(a_3, 0) + + a_10 = af.sum(E_upper * a_0, 0) + a_11 = af.sum(E_upper * a_1, 0) + a_12 = af.sum(E_upper * a_2, 0) + a_13 = af.sum(E_upper * a_3, 0) + + a_20 = af.sum(p1 * a_0, 0) + a_21 = af.sum(p1 * a_1, 0) + a_22 = af.sum(p1 * a_2, 0) + a_23 = af.sum(p1 * a_3, 0) + + a_30 = af.sum(p2 * a_0, 0) + a_31 = af.sum(p2 * a_1, 0) + a_32 = af.sum(p2 * a_2, 0) + a_33 = af.sum(p2 * a_3, 0) + + A = [ [a_00, a_01, a_02, a_03], \ + [a_10, a_11, a_12, a_13], \ + [a_20, a_21, a_22, a_23], \ + [a_30, a_31, a_32, a_33] \ + ] + + fermi_dirac = 1./(af.exp( ( E_upper - mu_ee + - vel_drift_x*p1 - vel_drift_y*p2 + )/(k*T_ee) + ) + 1. + ) + af.eval(fermi_dirac) + + zeroth_moment = (f - fermi_dirac) + second_moment = E_upper*(f - fermi_dirac) + first_moment_x = p1*(f - fermi_dirac) + first_moment_y = p2*(f - fermi_dirac) + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + eqn_energy_conservation = af.sum(second_moment, 0) + eqn_mom_x_conservation = af.sum(first_moment_x, 0) + eqn_mom_y_conservation = af.sum(first_moment_y, 0) + + residual = [eqn_mass_conservation, \ + eqn_energy_conservation, \ + eqn_mom_x_conservation, \ + eqn_mom_y_conservation] + + error_norm = np.max([af.max(af.abs(residual[0])), + af.max(af.abs(residual[1])), + af.max(af.abs(residual[2])), + af.max(af.abs(residual[3])) + ] + ) + print(" rank = ", params.rank, + "||residual_ee|| = ", error_norm + ) + +# if (error_norm < 1e-13): +# params.mu_ee = mu_ee +# params.T_ee = T_ee +# params.vel_drift_x = vel_drift_x +# params.vel_drift_y = vel_drift_y +# return(fermi_dirac) + + b_0 = eqn_mass_conservation + b_1 = eqn_energy_conservation + b_2 = eqn_mom_x_conservation + b_3 = eqn_mom_y_conservation + b = [b_0, b_1, b_2, b_3] + + # Solve Ax = b + # where A == Jacobian, + # x == delta guess (correction to guess), + # b = -residual + + A_inv = inverse_4x4_matrix(A) + + x_0 = A_inv[0][0]*b[0] + A_inv[0][1]*b[1] + A_inv[0][2]*b[2] + A_inv[0][3]*b[3] + x_1 = A_inv[1][0]*b[0] + A_inv[1][1]*b[1] + A_inv[1][2]*b[2] + A_inv[1][3]*b[3] + x_2 = A_inv[2][0]*b[0] + A_inv[2][1]*b[1] + A_inv[2][2]*b[2] + A_inv[2][3]*b[3] + x_3 = A_inv[3][0]*b[0] + A_inv[3][1]*b[1] + A_inv[3][2]*b[2] + A_inv[3][3]*b[3] + + delta_mu = x_0 + delta_T = x_1 + delta_vx = x_2 + delta_vy = x_3 + + mu_ee = mu_ee + delta_mu + T_ee = T_ee + delta_T + vel_drift_x = vel_drift_x + delta_vx + vel_drift_y = vel_drift_y + delta_vy + + af.eval(mu_ee, T_ee, vel_drift_x, vel_drift_y) + + # Solved for (mu_ee, T_ee, vel_drift_x, vel_drift_y). Now store in params + params.mu_ee = mu_ee + params.T_ee = T_ee + params.vel_drift_x = vel_drift_x + params.vel_drift_y = vel_drift_y + + fermi_dirac = 1./(af.exp( ( E_upper - mu_ee + - vel_drift_x*p1 - vel_drift_y*p2 + )/(k*T_ee) + ) + 1. + ) + af.eval(fermi_dirac) + + zeroth_moment = f - fermi_dirac + second_moment = E_upper*(f - fermi_dirac) + first_moment_x = p1*(f - fermi_dirac) + first_moment_y = p2*(f - fermi_dirac) + + eqn_mass_conservation = af.sum(zeroth_moment, 0) + eqn_energy_conservation = af.sum(second_moment, 0) + eqn_mom_x_conservation = af.sum(first_moment_x, 0) + eqn_mom_y_conservation = af.sum(first_moment_y, 0) + + residual = [eqn_mass_conservation, \ + eqn_energy_conservation, \ + eqn_mom_x_conservation, \ + eqn_mom_y_conservation + ] + + error_norm = np.max([af.max(af.abs(residual[0])), + af.max(af.abs(residual[1])), + af.max(af.abs(residual[2])), + af.max(af.abs(residual[3])) + ] + ) + print(" rank = ", params.rank, + "||residual_ee|| = ", error_norm + ) + N_g = domain.N_ghost + print(" rank = ", params.rank, + "mu_ee = ", af.mean(params.mu_ee[0, N_g:-N_g, N_g:-N_g]), + "T_ee = ", af.mean(params.T_ee[0, N_g:-N_g, N_g:-N_g]), + " = ", af.mean(params.vel_drift_x[0, N_g:-N_g, N_g:-N_g]), + " = ", af.mean(params.vel_drift_y[0, N_g:-N_g, N_g:-N_g]) + ) + PETSc.Sys.Print(" ------------------") + + return(fermi_dirac) + +def RTA(f, q1, q2, p1, p2, p3, moments, params, flag = False): + """Return BGK operator -(f-f0)/tau.""" + + if(af.any_true(params.tau_defect(q1, q2, p1, p2, p3) == 0)): + if (flag == False): + return(0.*f) + + f = f0_defect_constant_T(f, p1, p2, p3, params) + + return(f) + + C_f = -( f - f0_defect_constant_T(f, p1, p2, p3, params) \ + ) / params.tau_defect(q1, q2, p1, p2, p3) \ + -( f - f0_ee(f, p1, p2, p3, params) + ) / params.tau_ee(q1, q2, p1, p2, p3) + # When (f - f0) is NaN. Dividing by np.inf doesn't give 0 + # TODO: WORKAROUND + + af.eval(C_f) + return(C_f) + diff --git a/bolt/src/electronic_boltzmann/matrix_inverse.py b/bolt/src/electronic_boltzmann/matrix_inverse.py new file mode 100644 index 00000000..4b987180 --- /dev/null +++ b/bolt/src/electronic_boltzmann/matrix_inverse.py @@ -0,0 +1,189 @@ +import numpy as np +import arrayfire as af + +def inverse_4x4_matrix(A): +# TO TEST: +# A_test = np.random.rand(4, 4) +# A_inv_test = np.linalg.inv(A_test) +# A_inv = np.array(inverse_4x4_matrix(A_test)) +# print("err = ", np.max(np.abs(A_inv - A_inv_test))) + + + det = \ + A[0][0]*A[1][1]*A[2][2]*A[3][3] \ + + A[0][0]*A[1][2]*A[2][3]*A[3][1] \ + + A[0][0]*A[1][3]*A[2][1]*A[3][2] \ + + A[0][1]*A[1][0]*A[2][3]*A[3][2] \ + + A[0][1]*A[1][2]*A[2][0]*A[3][3] \ + + A[0][1]*A[1][3]*A[2][2]*A[3][0] \ + + A[0][2]*A[1][0]*A[2][1]*A[3][3] \ + + A[0][2]*A[1][1]*A[2][3]*A[3][0] \ + + A[0][2]*A[1][3]*A[2][0]*A[3][1] \ + + A[0][3]*A[1][0]*A[2][2]*A[3][1] \ + + A[0][3]*A[1][1]*A[2][0]*A[3][2] \ + + A[0][3]*A[1][2]*A[2][1]*A[3][0] \ + - A[0][0]*A[1][1]*A[2][3]*A[3][2] \ + - A[0][0]*A[1][2]*A[2][1]*A[3][3] \ + - A[0][0]*A[1][3]*A[2][2]*A[3][1] \ + - A[0][1]*A[1][0]*A[2][2]*A[3][3] \ + - A[0][1]*A[1][2]*A[2][3]*A[3][0] \ + - A[0][1]*A[1][3]*A[2][0]*A[3][2] \ + - A[0][2]*A[1][0]*A[2][3]*A[3][1] \ + - A[0][2]*A[1][1]*A[2][0]*A[3][3] \ + - A[0][2]*A[1][3]*A[2][1]*A[3][0] \ + - A[0][3]*A[1][0]*A[2][1]*A[3][2] \ + - A[0][3]*A[1][1]*A[2][2]*A[3][0] \ + - A[0][3]*A[1][2]*A[2][0]*A[3][1] + + af.eval(det) + + A_inv = [[0, 0, 0, 0], + [0, 0, 0, 0], + [0, 0, 0, 0], + [0, 0, 0, 0] + ] + + A_inv[0][0] = \ + ( A[1][1]*A[2][2]*A[3][3] + + A[1][2]*A[2][3]*A[3][1] + + A[1][3]*A[2][1]*A[3][2] + - A[1][1]*A[2][3]*A[3][2] + - A[1][2]*A[2][1]*A[3][3] + - A[1][3]*A[2][2]*A[3][1])/det + + A_inv[0][1] = \ + ( A[0][1]*A[2][3]*A[3][2] + + A[0][2]*A[2][1]*A[3][3] + + A[0][3]*A[2][2]*A[3][1] + - A[0][1]*A[2][2]*A[3][3] + - A[0][2]*A[2][3]*A[3][1] + - A[0][3]*A[2][1]*A[3][2])/det + + A_inv[0][2] = \ + ( A[0][1]*A[1][2]*A[3][3] + + A[0][2]*A[1][3]*A[3][1] + + A[0][3]*A[1][1]*A[3][2] + - A[0][1]*A[1][3]*A[3][2] + - A[0][2]*A[1][1]*A[3][3] + - A[0][3]*A[1][2]*A[3][1])/det + + A_inv[0][3] = \ + ( A[0][1]*A[1][3]*A[2][2] + + A[0][2]*A[1][1]*A[2][3] + + A[0][3]*A[1][2]*A[2][1] + - A[0][1]*A[1][2]*A[2][3] + - A[0][2]*A[1][3]*A[2][1] + - A[0][3]*A[1][1]*A[2][2])/det + + # b21 + A_inv[1][0] = \ + ( A[1][0]*A[2][3]*A[3][2] + + A[1][2]*A[2][0]*A[3][3] + + A[1][3]*A[2][2]*A[3][0] + - A[1][0]*A[2][2]*A[3][3] + - A[1][2]*A[2][3]*A[3][0] + - A[1][3]*A[2][0]*A[3][2])/det + + A_inv[1][1] = \ + ( A[0][0]*A[2][2]*A[3][3] + + A[0][2]*A[2][3]*A[3][0] + + A[0][3]*A[2][0]*A[3][2] + - A[0][0]*A[2][3]*A[3][2] + - A[0][2]*A[2][0]*A[3][3] + - A[0][3]*A[2][2]*A[3][0])/det + + A_inv[1][2] = \ + ( A[0][0]*A[1][3]*A[3][2] + + A[0][2]*A[1][0]*A[3][3] + + A[0][3]*A[1][2]*A[3][0] + - A[0][0]*A[1][2]*A[3][3] + - A[0][2]*A[1][3]*A[3][0] + - A[0][3]*A[1][0]*A[3][2])/det + + A_inv[1][3] = \ + ( A[0][0]*A[1][2]*A[2][3] + + A[0][2]*A[1][3]*A[2][0] + + A[0][3]*A[1][0]*A[2][2] + - A[0][0]*A[1][3]*A[2][2] + - A[0][2]*A[1][0]*A[2][3] + - A[0][3]*A[1][2]*A[2][0])/det + + # b31 + A_inv[2][0] = \ + ( A[1][0]*A[2][1]*A[3][3] + + A[1][1]*A[2][3]*A[3][0] + + A[1][3]*A[2][0]*A[3][1] + - A[1][0]*A[2][3]*A[3][1] + - A[1][1]*A[2][0]*A[3][3] + - A[1][3]*A[2][1]*A[3][0])/det + + # b32 + A_inv[2][1] = \ + ( A[0][0]*A[2][3]*A[3][1] + + A[0][1]*A[2][0]*A[3][3] + + A[0][3]*A[2][1]*A[3][0] + - A[0][0]*A[2][1]*A[3][3] + - A[0][1]*A[2][3]*A[3][0] + - A[0][3]*A[2][0]*A[3][1])/det + + A_inv[2][2] = \ + ( A[0][0]*A[1][1]*A[3][3] + + A[0][1]*A[1][3]*A[3][0] + + A[0][3]*A[1][0]*A[3][1] + - A[0][0]*A[1][3]*A[3][1] + - A[0][1]*A[1][0]*A[3][3] + - A[0][3]*A[1][1]*A[3][0])/det + + A_inv[2][3] = \ + ( A[0][0]*A[1][3]*A[2][1] + + A[0][1]*A[1][0]*A[2][3] + + A[0][3]*A[1][1]*A[2][0] + - A[0][0]*A[1][1]*A[2][3] + - A[0][1]*A[1][3]*A[2][0] + - A[0][3]*A[1][0]*A[2][1])/det + + # b41 + A_inv[3][0] = \ + ( A[1][0]*A[2][2]*A[3][1] + + A[1][1]*A[2][0]*A[3][2] + + A[1][2]*A[2][1]*A[3][0] + - A[1][0]*A[2][1]*A[3][2] + - A[1][1]*A[2][2]*A[3][0] + - A[1][2]*A[2][0]*A[3][1])/det + + # b42 + A_inv[3][1] = \ + ( A[0][0]*A[2][1]*A[3][2] + + A[0][1]*A[2][2]*A[3][0] + + A[0][2]*A[2][0]*A[3][1] + - A[0][0]*A[2][2]*A[3][1] + - A[0][1]*A[2][0]*A[3][2] + - A[0][2]*A[2][1]*A[3][0])/det + + # b43 + A_inv[3][2] = \ + ( A[0][0]*A[1][2]*A[3][1] + + A[0][1]*A[1][0]*A[3][2] + + A[0][2]*A[1][1]*A[3][0] + - A[0][0]*A[1][1]*A[3][2] + - A[0][1]*A[1][2]*A[3][0] + - A[0][2]*A[1][0]*A[3][1])/det + + A_inv[3][3] = \ + ( A[0][0]*A[1][1]*A[2][2] + + A[0][1]*A[1][2]*A[2][0] + + A[0][2]*A[1][0]*A[2][1] + - A[0][0]*A[1][2]*A[2][1] + - A[0][1]*A[1][0]*A[2][2] + - A[0][2]*A[1][1]*A[2][0])/det + + arrays_to_be_evaled = \ + [A_inv[0][0], A_inv[0][1], A_inv[0][2], A_inv[0][3], \ + A_inv[1][0], A_inv[1][1], A_inv[1][2], A_inv[1][3], \ + A_inv[2][0], A_inv[2][1], A_inv[2][2], A_inv[2][3], \ + A_inv[3][0], A_inv[3][1], A_inv[3][2], A_inv[3][3] \ + ] + + af.eval(*arrays_to_be_evaled) + + return(A_inv) diff --git a/bolt/src/electronic_boltzmann/moment_defs.py b/bolt/src/electronic_boltzmann/moment_defs.py new file mode 100644 index 00000000..9ac8c96e --- /dev/null +++ b/bolt/src/electronic_boltzmann/moment_defs.py @@ -0,0 +1,12 @@ +import numpy as np +import params + +moment_exponents = dict(density = [0, 0, 0], + j_x = [1, 0, 0], + j_y = [0, 1, 0] + ) + +moment_coeffs = dict(density = [4./(2.*np.pi*params.h_bar)**2., 0, 0], + j_x = [4./(2.*np.pi*params.h_bar)**2., 0, 0], + j_y = [0, 4./(2.*np.pi*params.h_bar)**2., 0] + ) diff --git a/example_problems/electronic_boltzmann/graphene/boundary_conditions.py b/example_problems/electronic_boltzmann/graphene/boundary_conditions.py new file mode 100644 index 00000000..fc1a9927 --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/boundary_conditions.py @@ -0,0 +1,88 @@ +import numpy as np +import arrayfire as af +import domain + +in_q1_left = 'mirror+dirichlet' +in_q1_right = 'mirror+dirichlet' +in_q2_bottom = 'mirror' +in_q2_top = 'mirror' + +@af.broadcast +def f_left(f, q1, q2, p1, p2, p3, params): + + k = params.boltzmann_constant + E_upper = params.E_band + T = params.initial_temperature + mu = params.initial_mu + + t = params.current_time + omega = 2. * np.pi * params.AC_freq + vel_drift_x_in = params.vel_drift_x_in * np.sin(omega*t) + + fermi_dirac_in = (1./(af.exp( (E_upper - vel_drift_x_in*p1 - mu)/(k*T) ) + 1.) + ) + + if (params.contact_geometry="straight"): + # Contacts on either side of the device + + q2_contact_start = params.contact_start + q2_contact_end = params.contact_end + + cond = ((q2 >= q2_contact_start) & \ + (q2 <= q2_contact_end) \ + ) + + f_left = cond*fermi_dirac_in + (1 - cond)*f + + elif (params.contact_geometry="turn_around"): + # Contacts on the same side of the device + + vel_drift_x_out = -params.vel_drift_x_in * np.sin(omega*t) + + fermi_dirac_out = (1./(af.exp( (E_upper - vel_drift_x_out*p1 - mu)/(k*T) ) + 1.) + ) + + # TODO: set these parameters in params.py + cond_in = ((q2 >= 3.5) & (q2 <= 4.5)) + cond_out = ((q2 >= 5.5) & (q2 <= 6.5)) + + f_left = cond_in*fermi_dirac_in + cond_out*fermi_dirac_out \ + + (1 - cond_in)*(1 - cond_out)*f + + af.eval(f_left) + return(f_left) + +@af.broadcast +def f_right(f, q1, q2, p1, p2, p3, params): + + k = params.boltzmann_constant + E_upper = params.E_band + T = params.initial_temperature + mu = params.initial_mu + + t = params.current_time + omega = 2. * np.pi * params.AC_freq + vel_drift_x_out = params.vel_drift_x_out * np.sin(omega*t) + + fermi_dirac_out = (1./(af.exp( (E_upper - vel_drift_x_out*p1 - mu)/(k*T) ) + 1.) + ) + + if (params.contact_geometry="straight"): + # Contacts on either side of the device + + q2_contact_start = params.contact_start + q2_contact_end = params.contact_end + + cond = ((q2 >= q2_contact_start) & \ + (q2 <= q2_contact_end) \ + ) + + f_right = cond*fermi_dirac_out + (1 - cond)*f + + elif (params.contact_geometry="turn_around"): + # Contacts on the same side of the device + + f_right = f + + af.eval(f_right) + return(f_right) diff --git a/example_problems/electronic_boltzmann/graphene/domain.py b/example_problems/electronic_boltzmann/graphene/domain.py new file mode 100644 index 00000000..36ff4b4c --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/domain.py @@ -0,0 +1,21 @@ +q1_start = 0. +q1_end = 5. +N_q1 = 60 + +q2_start = 0. +q2_end = 10. +N_q2 = 120 + +p1_start = -0.04 +p1_end = 0.04 +N_p1 = 32 + +p2_start = -0.04 +p2_end = 0.04 +N_p2 = 32 + +p3_start = -0.5 +p3_end = 0.5 +N_p3 = 1 + +N_ghost = 2 diff --git a/example_problems/electronic_boltzmann/graphene/initialize.py b/example_problems/electronic_boltzmann/graphene/initialize.py new file mode 100644 index 00000000..aef71ebe --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/initialize.py @@ -0,0 +1,56 @@ +""" +Functions which are used in assigning the I.C's to +the system. +""" + +import arrayfire as af +import numpy as np +from petsc4py import PETSc + +def initialize_f(q1, q2, p1, p2, p3, params): + + PETSc.Sys.Print("Initializing f") + k = params.boltzmann_constant + + params.mu = 0.*q1 + params.initial_mu + params.T = 0.*q1 + params.initial_temperature + params.vel_drift_x = 0.*q1 + params.vel_drift_y = 0.*q1 + params.phi = 0.*q1 + + params.mu_ee = params.mu.copy() + params.T_ee = params.T.copy() + params.vel_drift_x = 0.*q1 + 0e-3 + params.vel_drift_y = 0.*q1 + 0e-3 + + params.E_band = params.band_energy(p1, p2) + params.vel_band = params.band_velocity(p1, p2) + + E_upper = params.E_band + params.charge_electron*params.phi + + f = (1./(af.exp( (E_upper - params.vel_drift_x*p1 + - params.vel_drift_y*p2 + - params.mu + )/(k*params.T) + ) + 1. + )) + + af.eval(f) + return(f) + + +def initialize_E(q1, q2, params): + + E1 = 0.*q1 + E2 = 0.*q1 + E3 = 0.*q1 + + return(E1, E2, E3) + +def initialize_B(q1, q2, params): + + B1 = 0.*q1 + B2 = 0.*q1 + B3 = 0.*q1 + + return(B1, B2, B3) diff --git a/example_problems/electronic_boltzmann/graphene/job_script b/example_problems/electronic_boltzmann/graphene/job_script new file mode 100644 index 00000000..8cb37850 --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/job_script @@ -0,0 +1,4 @@ +#!/bin/bash +#SBATCH -p gpu -n 2 --gres=gpu:2 -c1 --hint=nomultithread -t 2-00:00 +# (Note: use one MPI process per gpu, update both gres and -n together; max 6) +time mpirun python main.py -snes_monitor -snes_max_it 1 -snes_lag_jacobian_persists TRUE -snes_lag_jacobian 1000000 -snes_atol 1e-50 -snes_rtol 1e-50 > output.txt diff --git a/example_problems/electronic_boltzmann/graphene/main.py b/example_problems/electronic_boltzmann/graphene/main.py new file mode 100644 index 00000000..3c243bed --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/main.py @@ -0,0 +1,376 @@ +import sys +import arrayfire as af +import numpy as np +import pylab as pl +import h5py +import petsc4py, sys; petsc4py.init(sys.argv) +from petsc4py import PETSc + +from bolt.lib.physical_system import physical_system + +from bolt.lib.nonlinear_solver.nonlinear_solver \ + import nonlinear_solver +from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic \ + import compute_electrostatic_fields + +import domain +import boundary_conditions +import params +import initialize + +import bolt.src.electronic_boltzmann.advection_terms as advection_terms + +import bolt.src.electronic_boltzmann.collision_operator \ + as collision_operator + +import bolt.src.electronic_boltzmann.moment_defs as moment_defs + +#pl.rcParams['figure.figsize'] = 12, 7.5 +#pl.rcParams['figure.dpi'] = 150 +#pl.rcParams['image.cmap'] = 'jet' +#pl.rcParams['lines.linewidth'] = 1.5 +#pl.rcParams['font.family'] = 'serif' +#pl.rcParams['font.weight'] = 'bold' +#pl.rcParams['font.size'] = 20 +#pl.rcParams['font.sans-serif'] = 'serif' +#pl.rcParams['text.usetex'] = False +#pl.rcParams['axes.linewidth'] = 1.5 +#pl.rcParams['axes.titlesize'] = 'medium' +#pl.rcParams['axes.labelsize'] = 'medium' +# +#pl.rcParams['xtick.major.size'] = 8 +#pl.rcParams['xtick.minor.size'] = 4 +#pl.rcParams['xtick.major.pad'] = 8 +#pl.rcParams['xtick.minor.pad'] = 8 +#pl.rcParams['xtick.color'] = 'k' +#pl.rcParams['xtick.labelsize'] = 'medium' +#pl.rcParams['xtick.direction'] = 'in' +# +#pl.rcParams['ytick.major.size'] = 8 +#pl.rcParams['ytick.minor.size'] = 4 +#pl.rcParams['ytick.major.pad'] = 8 +#pl.rcParams['ytick.minor.pad'] = 8 +#pl.rcParams['ytick.color'] = 'k' +#pl.rcParams['ytick.labelsize'] = 'medium' +#pl.rcParams['ytick.direction'] = 'in' + + +# Defining the physical system to be solved: +system = physical_system(domain, + boundary_conditions, + params, + initialize, + advection_terms, + collision_operator.RTA, + moment_defs + ) + +# Time parameters: +dt = params.dt +t_final = params.t_final +params.current_time = t0 = 0.0 +params.time_step = time_step = 0 +dump_counter = 0 +dump_time_array = [] + +N_g = domain.N_ghost + +# Declaring a nonlinear system object which will evolve the defined physical system: +nls = nonlinear_solver(system) +params.rank = nls._comm.rank + +if (params.restart): + nls.load_distribution_function(params.restart_file) + + viewer = PETSc.Viewer().createHDF5(params.phi_restart_file, + PETSc.Viewer.Mode.READ, + comm=nls._comm + ) + PETSc.Object.setName(nls.poisson.glob_phi, 'phi') + nls.poisson.glob_phi.load(viewer) + + params.solve_for_equilibrium = 0 + compute_electrostatic_fields(nls) + + params.mu = params.global_chem_potential + params.charge_electron*params.phi + params.mu_analytic = params.global_chem_potential + params.charge_electron*params.phi + params.mu = af.moddims(params.mu, + 1, + nls.N_q1_local + 2*N_g, + nls.N_q2_local + 2*N_g + ) +else: + if (params.solve_for_equilibrium): + PETSc.Sys.Print("=============================") + PETSc.Sys.Print("Solving for equilibrium state") + PETSc.Sys.Print("=============================") + params.solve_for_equilibrium = 1 + compute_electrostatic_fields(nls) + params.solve_for_equilibrium = 0 + + params.mu = params.global_chem_potential + params.charge_electron*params.phi + params.mu = af.moddims(params.mu, + 1, + nls.N_q1_local + 2*N_g, + nls.N_q2_local + 2*N_g + ) + + nls.f = params.fermi_dirac(params.mu, params.E_band) + + PETSc.Sys.Print("=====================") + PETSc.Sys.Print("Solution diagnostics:") + PETSc.Sys.Print("=====================") + + density = nls.compute_moments('density') + print("rank = ", params.rank, "\n", + " = ", af.mean(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(mu) = ", af.max(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + #" = ", af.mean(params.phi[N_g:-N_g, N_g:-N_g]), "\n", + " = ", af.mean(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(n) = ", af.max(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " |E1| = ", af.mean(af.abs(nls.cell_centered_EM_fields[0, N_g:-N_g, N_g:-N_g])), + "\n", + " |E2| = ", af.mean(af.abs(nls.cell_centered_EM_fields[1, N_g:-N_g, N_g:-N_g])) + ) + + PETSc.Sys.Print("===============") + PETSc.Sys.Print("Dumping data...") + PETSc.Sys.Print("===============") + + nls.dump_moments('dumps/moments_eqbm') + nls.dump_distribution_function('dumps/f_eqbm') + + # Dump EM fields + af.flat(nls.cell_centered_EM_fields[:, N_g:-N_g, N_g:-N_g]).to_ndarray(nls._glob_fields_array) + PETSc.Object.setName(nls._glob_fields, 'fields') + viewer = PETSc.Viewer().createHDF5('dumps/fields_eqbm.h5', + 'w', comm=nls._comm + ) + viewer(nls._glob_fields) + + # Dump eqbm potential + PETSc.Object.setName(nls.poisson.glob_phi, 'phi') + viewer = PETSc.Viewer().createHDF5('dumps/phi_eqbm.h5', + 'w', comm=nls._comm + ) + viewer(nls.poisson.glob_phi) + + sys.exit("Terminating execution") + else: + + #compute_electrostatic_fields(nls) + pass + +density = nls.compute_moments('density') +print("rank = ", params.rank, "\n", + " = ", af.mean(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(mu) = ", af.max(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + #" = ", af.mean(params.phi[N_g:-N_g, N_g:-N_g]), "\n", + " = ", af.mean(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(n) = ", af.max(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " |E1| = ", af.mean(af.abs(nls.cell_centered_EM_fields[0, N_g:-N_g, N_g:-N_g])), + "\n", + " |E2| = ", af.mean(af.abs(nls.cell_centered_EM_fields[1, N_g:-N_g, N_g:-N_g])) + ) + +nls.f = af.select(nls.f < 1e-20, 1e-20, nls.f) +while t0 < t_final: + + # Refine to machine error + if (time_step==0): + params.collision_nonlinear_iters = 10 + else: + params.collision_nonlinear_iters = params.collision_operator_nonlinear_iters + + dump_steps = params.dump_steps + # Uncomment if need to dump more frequently during a desired time interval + #if (params.current_time > 149. and params.current_time < 154): + # dump_steps = 1 + #else: + # dump_steps = params.dump_steps + if (time_step%dump_steps==0): + file_number = '%06d'%dump_counter + dump_counter= dump_counter + 1 + dump_time_array.append(params.current_time) + PETSc.Sys.Print("=====================================================") + PETSc.Sys.Print("Dumping data at time step =", time_step, + ", file number =", file_number + ) + PETSc.Sys.Print("=====================================================") + nls.dump_moments('dumps/moments_' + file_number) + #nls.dump_distribution_function('dumps/f_' + file_number) + nls.dump_aux_arrays([params.mu, + params.mu_ee, + params.T_ee, + params.vel_drift_x, params.vel_drift_y], + 'lagrange_multipliers', + 'dumps/lagrange_multipliers_' + file_number + ) + + # Dump EM fields + af.flat(nls.cell_centered_EM_fields[:, N_g:-N_g, N_g:-N_g]).to_ndarray(nls._glob_fields_array) + PETSc.Object.setName(nls._glob_fields, 'fields') + viewer = PETSc.Viewer().createHDF5('dumps/fields_' + + '%06d'%(time_step/dump_steps) + '.h5', + 'w', comm=nls._comm + ) + viewer(nls._glob_fields) + + + dt_force_constraint = 0. +# dt_force_constraint = \ +# 0.5 * np.min(nls.dp1, nls.dp2) \ +# / np.max((af.max(nls.cell_centered_EM_fields[0]), +# af.max(nls.cell_centered_EM_fields[1]) +# ) +# ) + + + PETSc.Sys.Print("Time step =", time_step, ", Time =", t0) + + mean_density = af.mean(density[0, N_g:-N_g, N_g:-N_g]) + density_pert = density - mean_density + +# if (time_step%1==0): +# chemical_potential = np.array(params.mu)[0, :, :] \ +# - params.charge_electron*np.array(params.phi)[:, :] +# pl.contourf(np.array(nls.q1_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.q2_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(params.mu)[0, N_g:-N_g, N_g:-N_g], \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/mu_' + '%06d'%time_step + '.png' ) +# pl.clf() +# +# pl.contourf(np.array(nls.q1_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.q2_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(density)[0, N_g:-N_g, N_g:-N_g], \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/density_' + '%06d'%time_step + '.png' ) +# pl.clf() +# +# pl.contourf(np.array(nls.q1_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.q2_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(params.phi)[N_g:-N_g, N_g:-N_g], \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/phi_' + '%06d'%time_step + '.png' ) +# pl.clf() +# +# pl.contourf(np.array(nls.q1_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.q2_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.cell_centered_EM_fields)[0, N_g:-N_g, N_g:-N_g], \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/E1_' + '%06d'%time_step + '.png' ) +# pl.clf() +# +# pl.contourf(np.array(nls.q1_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.q2_center)[0, N_g:-N_g, N_g:-N_g], \ +# np.array(nls.cell_centered_EM_fields)[1, N_g:-N_g, N_g:-N_g], \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/E2_' + '%06d'%time_step + '.png' ) +# pl.clf() +# +# f_at_desired_q = af.moddims(nls.f[:, N_g, N_g + 0.*nls.N_q2/2], +# nls.N_p1, nls.N_p2 +# ) +# p1 = af.moddims(nls.p1, nls.N_p1, nls.N_p2) +# p2 = af.moddims(nls.p2, nls.N_p1, nls.N_p2) +# pl.contourf(np.array(p1), \ +# np.array(p2), \ +# np.array((f_at_desired_q)), \ +# 100, cmap='bwr' +# ) +# pl.title('Time = ' + "%.2f"%(t0) ) +# pl.xlabel('$x$') +# pl.ylabel('$y$') +# pl.colorbar() +# pl.gca().set_aspect('equal') +# pl.savefig('/tmp/f_' + '%06d'%time_step + '.png' ) +# pl.clf() + + nls.strang_timestep(dt) + t0 = t0 + dt + time_step = time_step + 1 + params.time_step = time_step + params.current_time = t0 + + # Floors + nls.f = af.select(nls.f < 1e-20, 1e-20, nls.f) + + density = nls.compute_moments('density') + print("rank = ", params.rank, "\n", + " = ", af.mean(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(mu) = ", af.max(params.mu[0, N_g:-N_g, N_g:-N_g]), "\n", + #" = ", af.mean(params.phi[N_g:-N_g, N_g:-N_g]), "\n", + " = ", af.mean(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " max(n) = ", af.max(density[0, N_g:-N_g, N_g:-N_g]), "\n", + " |E1| = ", af.mean(af.abs(nls.cell_centered_EM_fields[0, N_g:-N_g, N_g:-N_g])), + "\n", + " |E2| = ", af.mean(af.abs(nls.cell_centered_EM_fields[1, N_g:-N_g, N_g:-N_g])) + ) + PETSc.Sys.Print("--------------------\n") + +if (params.rank==0): + np.savetxt("dump_time_array.txt", dump_time_array) + +#phi_array = nls.poisson.glob_phi.getArray() +#phi_array = phi_array.reshape([nls.poisson.N_q3_3D_local, \ +# nls.poisson.N_q2_3D_local, \ +# nls.poisson.N_q1_3D_local] +# ) +# +#q3_index = np.where(nls.poisson.q3_3D[N_g:-N_g] >= nls.location_in_q3)[0][0] +#pl.plot(nls.poisson.q1_3D[N_g:-N_g], +# phi_array[q3_index, nls.poisson.N_q2_3D_local/2, :], 'o-' +# ) +#pl.show() + +#pl.rcParams['figure.figsize'] = 20, 7.5 +#pl.subplot(121) +#N_g = domain.N_ghost +#pl.contourf( +# phi_array[nls.N_q3_poisson/2, :, :], 100, cmap='jet' +# ) +#pl.colorbar() +#pl.title('Top View') +#pl.xlabel('$x$') +#pl.ylabel('$y$') +#pl.gca().set_aspect('equal') +# +#pl.subplot(122) +#pl.contourf(phi_array[:, nls.N_q2_poisson/2, :], 100, cmap='jet') +#pl.title('Side View') +#pl.xlabel('$x$') +#pl.ylabel('$z$') +#pl.colorbar() +#pl.gca().set_aspect('equal') +#pl.show() diff --git a/example_problems/electronic_boltzmann/graphene/params.py b/example_problems/electronic_boltzmann/graphene/params.py new file mode 100644 index 00000000..02c08f3b --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/params.py @@ -0,0 +1,134 @@ +import numpy as np +import arrayfire as af + +# Can be defined as 'electrostatic', 'user-defined'. +# The initial conditions need to be specified under initialize +# Ensure that the initial conditions specified satisfy +# Maxwell's constraint equations +fields_initialize = 'user-defined' + +# Can be defined as 'electrostatic' and 'fdtd' +fields_type = 'electrostatic' +fields_solver = 'SNES' + +# Can be defined as 'strang' and 'lie' +time_splitting = 'strang' + +# Method in q-space +solver_method_in_q = 'FVM' +solver_method_in_p = 'FVM' + +reconstruction_method_in_q = 'minmod' +reconstruction_method_in_p = 'minmod' + +riemann_solver = 'upwind-flux' + +# Restart(Set to zero for no-restart): +restart = 0 +restart_file = '/home/mani/work/quazar_research/bolt/example_problems/electronic_boltzmann/graphene/dumps/f_eqbm.h5' +phi_restart_file = '/home/mani/work/quazar_research/bolt/example_problems/electronic_boltzmann/graphene/dumps/phi_eqbm.h5' +electrostatic_solver_every_nth_step = 1000000 +solve_for_equilibrium = 0 + + +# File-writing Parameters: +dump_steps = 10 + +# Time parameters: +dt = 0.025/2 # ps +t_final = 200 # ps + +# Dimensionality considered in velocity space: +p_dim = 2 + +# Number of devices(GPUs/Accelerators) on each node: +num_devices = 6 + +# Constants: +mass_particle = 0.910938356 # x 1e-30 kg +h_bar = 1.0545718e-4 # x aJ ps +boltzmann_constant = 1 +charge_electron = 0.*-0.160217662 # x aC +speed_of_light = 300. # x [um/ps] +fermi_velocity = speed_of_light/300 +epsilon0 = 8.854187817 # x [aC^2 / (aJ um) ] + +epsilon_relative = 3.9 # SiO2 +backgate_potential = -10 # V +global_chem_potential = 0.03 +contact_start = 4.5 # um +contact_end = 5.5 # um +contact_geometry = "straight" # Contacts on either side of the device + # For contacts on the same side, use + # contact_geometry = "turn_around" + +initial_temperature = 12e-4 +initial_mu = 0.015 +vel_drift_x_in = 1e-4*fermi_velocity +vel_drift_x_out = 1e-4*fermi_velocity +AC_freq = 1./100 # ps^-1 + +B3_mean = 1. # T + +# Spatial quantities (will be initialized to shape = [q1, q2] in initalize.py) +mu = None # chemical potential used in the e-ph operator +T = None # Electron temperature used in the e-ph operator +mu_ee = None # chemical potential used in the e-e operator +T_ee = None # Electron temperature used in the e-e operator +vel_drift_x = None +vel_drift_y = None +phi = None # Electric potential in the plane of graphene sheet + +# Momentum quantities (will be initialized to shape = [p1*p2*p3] in initialize.py) +E_band = None +vel_band = None + +collision_operator_nonlinear_iters = 2 + +# Variation of collisional-timescale parameter through phase space: +@af.broadcast +def tau_defect(q1, q2, p1, p2, p3): + return(1. * q1**0 * p1**0) + +@af.broadcast +def tau_ee(q1, q2, p1, p2, p3): + return(0.2 * q1**0 * p1**0) + +def tau(q1, q2, p1, p2, p3): + return(tau_defect(q1, q2, p1, p2, p3)) + +def band_energy(p_x, p_y): + + p = af.sqrt(p_x**2. + p_y**2.) + + E_upper = p*fermi_velocity + + af.eval(E_upper) + return(E_upper) + +def band_velocity(p_x, p_y): + + p = af.sqrt(p_x**2. + p_y**2.) + p_hat = [p_x / (p + 1e-20), p_y / (p + 1e-20)] + + v_f = fermi_velocity + + upper_band_velocity = [ v_f * p_hat[0], v_f * p_hat[1]] + + af.eval(upper_band_velocity[0], upper_band_velocity[1]) + return(upper_band_velocity) + +@af.broadcast +def fermi_dirac(mu, E_band): + + k = boltzmann_constant + T = initial_temperature + + f = (1./(af.exp( (E_band - mu + )/(k*T) + ) + 1. + ) + ) + + af.eval(f) + return(f) diff --git a/example_problems/electronic_boltzmann/graphene/post.py b/example_problems/electronic_boltzmann/graphene/post.py new file mode 100644 index 00000000..a297c6f0 --- /dev/null +++ b/example_problems/electronic_boltzmann/graphene/post.py @@ -0,0 +1,422 @@ +import arrayfire as af +import numpy as np +from scipy.signal import correlate +import glob +import h5py +import matplotlib +import matplotlib.gridspec as gridspec +import matplotlib.patches as patches +matplotlib.use('agg') +import pylab as pl +import yt +yt.enable_parallelism() + +import petsc4py, sys; petsc4py.init(sys.argv) +from petsc4py import PETSc + +from bolt.lib.physical_system import physical_system + +from bolt.lib.nonlinear_solver.nonlinear_solver \ + import nonlinear_solver +from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic \ + import compute_electrostatic_fields + +import domain +import boundary_conditions +import params +import initialize + +import bolt.src.electronic_boltzmann.advection_terms as advection_terms + +import bolt.src.electronic_boltzmann.collision_operator \ + as collision_operator + +import bolt.src.electronic_boltzmann.moment_defs as moment_defs + +# Optimized plot parameters to make beautiful plots: +#pl.rcParams['figure.figsize'] = 8, 7.5 +pl.rcParams['figure.figsize'] = 8, 8 +#pl.rcParams['figure.figsize'] = 17, 9.5 +pl.rcParams['figure.dpi'] = 100 +pl.rcParams['image.cmap'] = 'jet' +pl.rcParams['lines.linewidth'] = 1.5 +#pl.rcParams['lines.linewidth'] = 3 +pl.rcParams['font.family'] = 'serif' +pl.rcParams['font.weight'] = 'bold' +pl.rcParams['font.size'] = 25 +pl.rcParams['font.sans-serif'] = 'serif' +pl.rcParams['text.usetex'] = True +pl.rcParams['axes.linewidth'] = 1.5 +pl.rcParams['axes.titlesize'] = 'medium' +pl.rcParams['axes.labelsize'] = 'medium' + +pl.rcParams['xtick.major.size'] = 8 +pl.rcParams['xtick.minor.size'] = 4 +pl.rcParams['xtick.major.pad'] = 8 +pl.rcParams['xtick.minor.pad'] = 8 +pl.rcParams['xtick.color'] = 'k' +pl.rcParams['xtick.labelsize'] = 'medium' +pl.rcParams['xtick.direction'] = 'in' + +pl.rcParams['ytick.major.size'] = 8 +pl.rcParams['ytick.minor.size'] = 4 +pl.rcParams['ytick.major.pad'] = 8 +pl.rcParams['ytick.minor.pad'] = 8 +pl.rcParams['ytick.color'] = 'k' +pl.rcParams['ytick.labelsize'] = 'medium' +pl.rcParams['ytick.direction'] = 'in' + +N_q1 = domain.N_q1 +N_q2 = domain.N_q2 +#N_q1 = 120 +#N_q2 = 240 + +q1 = domain.q1_start + (0.5 + np.arange(N_q1)) * (domain.q1_end - domain.q1_start)/N_q1 +q2 = domain.q2_start + (0.5 + np.arange(N_q2)) * (domain.q2_end - domain.q2_start)/N_q2 + +q2_meshgrid, q1_meshgrid = np.meshgrid(q2, q1) + +source_start = params.contact_start +source_end = params.contact_end + +drain_start = params.contact_start +drain_end = params.contact_end + +#source_start = 3.5; source_end = 4.5 +#drain_start = 5.5; drain_end = 6.5 + +source_indices = (q2 > source_start) & (q2 < source_end) +drain_indices = (q2 > drain_start) & (q2 < drain_end ) + +sensor_1_left_start = 8.5 # um +sensor_1_left_end = 9.5 # um + +sensor_1_right_start = 8.5 # um +sensor_1_right_end = 9.5 # um + +# Left needs to be near source, right sensor near drain +#sensor_1_left_start = 1.5 # um +#sensor_1_left_end = 2.5 # um + +#sensor_1_right_start = 7.5 # um +#sensor_1_right_end = 8.5 # um + +sensor_1_left_indices = (q2 > sensor_1_left_start ) & (q2 < sensor_1_left_end) +sensor_1_right_indices = (q2 > sensor_1_right_start) & (q2 < sensor_1_right_end) + +sensor_2_left_start = 6.5 # um +sensor_2_left_end = 7.5 # um + +sensor_2_right_start = 6.5 # um +sensor_2_right_end = 7.5 # um + +sensor_2_left_indices = (q2 > sensor_2_left_start ) & (q2 < sensor_2_left_end) +sensor_2_right_indices = (q2 > sensor_2_right_start) & (q2 < sensor_2_right_end) + +#dump_index = 0 +#h5f = h5py.File('dumps/moments_' + '%06d'%(dump_index) + '.h5', 'r') +#moments = np.swapaxes(h5f['moments'][:], 0, 1) +#h5f.close() +# +#density = moments[:, :, 0] +#j_x = moments[:, :, 1] +#j_y = moments[:, :, 2] +#pl.contourf(q1, q2, density, 100) +##pl.title('Time = ' + "%.2f"%(t0)) +#pl.axes().set_aspect('equal') +#pl.xlabel(r'$x$') +#pl.ylabel(r'$y$') +#pl.colorbar() +#pl.savefig('images/density' + '.png') +#pl.clf() +# +#h5f = h5py.File('dumps/lagrange_multipliers_' + '%06d'%(dump_index) + '.h5', 'r') +#lagrange_multipliers = np.swapaxes(h5f['lagrange_multipliers'][:], 0, 1) +#h5f.close() +# +#print("lagrange_multipliers.shape = ", lagrange_multipliers.shape) +#mu = lagrange_multipliers[:, :, 0] +#mu_ee = lagrange_multipliers[:, :, 1] +#T_ee = lagrange_multipliers[:, :, 2] +#vel_drift_x = lagrange_multipliers[:, :, 3] +#vel_drift_y = lagrange_multipliers[:, :, 4] +#j_x_prime = density*vel_drift_x +#print("err = ", np.mean(np.abs(j_x_prime - j_x))) + +#filepath = \ +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/geom_1/DC/tau_D_50_tau_ee_0.2' +#dump_file= np.sort(glob.glob(filepath+'/moment*.h5'))[-1] +# +#h5f = h5py.File(dump_file, 'r') +#moments = np.swapaxes(h5f['moments'][:], 0, 1) +#h5f.close() +# +#density = moments[:, :, 0] +#np.savetxt('paper_plots/density_tau_D_50_tau_ee_0.2.txt', density) +#np.savetxt('paper_plots/q2_DC_tau_D_50_tau_ee_0.2.txt', q2) +##pl.plot(q2[q2>source_end], density[0, q2>source_end]-np.mean(density)) +# +#filepath = \ +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/geom_1/DC/tau_D_5_tau_ee_0.2' +#dump_file= np.sort(glob.glob(filepath+'/moment*.h5'))[-1] +# +# +#h5f = h5py.File(dump_file, 'r') +#moments = np.swapaxes(h5f['moments'][:], 0, 1) +#h5f.close() +# +#density = moments[:, :, 0] +#np.savetxt('paper_plots/density_tau_D_5_tau_ee_0.2.txt', density) +#np.savetxt('paper_plots/q2_DC_tau_D_5_tau_ee_0.2.txt', q2) +# +#filepath = \ +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/geom_1/DC/tau_D_10_tau_ee_0.2' +#dump_file= np.sort(glob.glob(filepath+'/moment*.h5'))[-1] +# +#N_q1 = 120 +#N_q2 = 240 +# +#q1 = domain.q1_start + (0.5 + np.arange(N_q1)) * (domain.q1_end - domain.q1_start)/N_q1 +#q2 = domain.q2_start + (0.5 + np.arange(N_q2)) * (domain.q2_end - domain.q2_start)/N_q2 +# +#h5f = h5py.File(dump_file, 'r') +#moments = np.swapaxes(h5f['moments'][:], 0, 1) +#h5f.close() +# +#density = moments[:, :, 0] +#np.savetxt('paper_plots/density_tau_D_10_tau_ee_0.2.txt', density) +#np.savetxt('paper_plots/q2_DC_tau_D_10_tau_ee_0.2.txt', q2) + + +#pl.plot(q2[q2>source_end], density[0, q2>source_end]-np.mean(density)) +#pl.axhline(0, color='black', linestyle='--') +#pl.legend(['$\\tau_{ee}=0.2$ ps, $\\tau_{e-ph}=50$ ps', +# '$\\tau_{ee}=0.2$ ps, $\\tau_{e-ph}=10$ ps'], loc='lower right') +#pl.xlabel(r'$x\;(\mu \mathrm{m})$') +#pl.ylabel(r'$R\; (\mathrm{a.u.})$') +#pl.xlim(xmin=(source_end+0.1)) +#pl.savefig('paper_plots/DC.png') + +filepath = \ +'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/dumps_tau_D_1_tau_ee_0.2_movie' +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/geom_1/55_GHz/tau_D_5_tau_ee_0.2' +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/dumps_AC_10_Ghz_tau_D_10_tau_ee_1_geom_2/' +#'/home/mchandra/bolt/example_problems/electronic_boltzmann/graphene/dumps_tau_D_2_tau_ee_1_AC/' +moment_files = np.sort(glob.glob(filepath+'/moment*.h5')) +lagrange_multiplier_files = \ + np.sort(glob.glob(filepath+'/lagrange_multipliers*.h5')) + +dt = params.dt +dump_interval = params.dump_steps + +sensor_1_signal_array = [] +print("Reading sensor signal...") +for file_number, dump_file in enumerate(moment_files): + + h5f = h5py.File(dump_file, 'r') + moments = np.swapaxes(h5f['moments'][:], 0, 1) + h5f.close() + + density = moments[:, :, 0] + + source = np.mean(density[0, source_indices]) + drain = np.mean(density[-1, drain_indices]) + + sensor_1_left = np.mean(density[0, sensor_1_left_indices] ) + sensor_1_right = np.mean(density[-1, sensor_1_right_indices]) + + sensor_1_signal = sensor_1_left - sensor_1_right + + sensor_1_signal_array.append(sensor_1_signal) + +time_array = np.loadtxt("dump_time_array.txt") +AC_freq = 1./100 +input_signal_array = np.sin(2.*np.pi*AC_freq*time_array) +sensor_1_signal_array = np.array(sensor_1_signal_array) +half_time = (int)(time_array.size/2) +sensor_normalized = \ + sensor_1_signal_array/np.max(np.abs(sensor_1_signal_array[half_time:])) + +pl.rcParams['figure.figsize'] = 10, 8 +for file_number, dump_file in yt.parallel_objects(enumerate(moment_files)): + + print("file number = ", file_number, "of ", moment_files.size) + + h5f = h5py.File(dump_file, 'r') + moments = np.swapaxes(h5f['moments'][:], 0, 1) + h5f.close() + + gs = gridspec.GridSpec(3, 2) + pl.subplot(gs[:, 0]) + + density = moments[:, :, 0] + j_x = moments[:, :, 1] + j_y = moments[:, :, 2] + pl.contourf(q1_meshgrid, q2_meshgrid, density, 100, cmap='bwr') + pl.title(r'Time = ' + "%.2f"%(time_array[file_number]) + " ps") + #pl.colorbar() + + h5f = h5py.File(lagrange_multiplier_files[file_number], 'r') + lagrange_multipliers = h5f['lagrange_multipliers'][:] + h5f.close() + + mu = lagrange_multipliers[:, :, 0] + mu_ee = lagrange_multipliers[:, :, 1] + T_ee = lagrange_multipliers[:, :, 2] + vel_drift_x = lagrange_multipliers[:, :, 3] + vel_drift_y = lagrange_multipliers[:, :, 4] + +# pl.streamplot(q1[(int)(N_q1/2):], q2, +# vel_drift_x[:, (int)(N_q1/2):], vel_drift_y[:, (int)(N_q1/2):], +# density=2, color='blue', +# linewidth=0.7, arrowsize=1 +# ) + pl.streamplot(q1, q2, + vel_drift_x, vel_drift_y, + density=2, color='blue', + linewidth=0.7, arrowsize=1 + ) +# pl.streamplot(q1, q2, +# vel_drift_x, vel_drift_y, +# density=3, color='blue', +# linewidth=0.8, arrowsize=1.1 +# ) + pl.xlim([domain.q1_start, domain.q1_end]) + pl.ylim([domain.q2_start, domain.q2_end]) + #pl.ylim([0, 5]) + pl.gca().set_aspect('equal') + pl.xlabel(r'$x\;(\mu \mathrm{m})$') + pl.ylabel(r'$y\;(\mu \mathrm{m})$') + + pl.gca().annotate("+", xy=(-0.07, .9), xycoords=("axes fraction"), + ha="center", va="center", size=30, + bbox=dict(fc="white")) + pl.gca().annotate("-", xy=(1.05, .9), xycoords=("axes fraction"), + ha="center", va="center", size=30, + bbox=dict(fc="white", pad=6.5)) + + + pl.subplot(gs[1, 1]) + + pl.plot(time_array, input_signal_array) + pl.plot(time_array, sensor_normalized) + pl.axhline(0, color='black', linestyle='--') + pl.axvline(time_array[file_number], color='black', alpha=0.75) + pl.legend(['Source $I(t)$', 'Measured $V(t)$'], loc=(0.04, 1.125)) + pl.xlabel(r'Time (ps)') + pl.xlim([100, 200]) + pl.ylim([-1.1, 1.1]) + + + pl.suptitle('$\\tau_\mathrm{mc} = 0.2$ ps, $\\tau_\mathrm{mr} = 1$ ps') + #pl.tight_layout() + pl.savefig('images/dump_' + '%06d'%file_number + '.png') + #pl.savefig('paper_plots/DC.png') + pl.clf() + + +#time_array = [] +#input_signal_array = [] +#sensor_1_signal_array = [] +#sensor_2_signal_array = [] +#for file_number, dump_file in enumerate(moment_files): +# +# print("file number = ", file_number, "of ", moment_files.size) +# +# h5f = h5py.File(dump_file, 'r') +# moments = np.swapaxes(h5f['moments'][:], 0, 1) +# h5f.close() +# +# density = moments[:, :, 0] +# +# source = np.mean(density[0, source_indices]) +# drain = np.mean(density[-1, drain_indices]) +# +# sensor_1_left = np.mean(density[0, sensor_1_left_indices] ) +# sensor_1_right = np.mean(density[-1, sensor_1_right_indices]) +# #sensor_1_right = np.mean(density[0, sensor_1_right_indices]) +# +# sensor_2_left = np.mean(density[0, sensor_2_left_indices] ) +# sensor_2_right = np.mean(density[-1, sensor_2_right_indices]) +# +# #sensor_1_left = density[0, q2>source_end] +# #sensor_1_right = density[-1, q2>source_end] +# +# input_signal = source - drain +# sensor_1_signal = sensor_1_left - sensor_1_right +# sensor_2_signal = sensor_2_left - sensor_2_right +# +# time_array.append(file_number*dt*dump_interval) +# input_signal_array.append(input_signal) +# sensor_1_signal_array.append(sensor_1_signal) +# sensor_2_signal_array.append(sensor_2_signal) +# +##pl.rcParams['figure.figsize'] = 12, 8 +## +#AC_freq = 5.5/100 +#time_array = np.array(time_array) +#input_signal_array = np.sin(2.*np.pi*AC_freq*time_array) +#sensor_1_signal_array = np.array(sensor_1_signal_array) +###np.savetxt('drive.txt', input_signal_array) +##np.savetxt('paper_plots/sensor_tau_ee_0.2_tau_D_5.txt', sensor_1_signal_array) +##np.savetxt('time.txt', time_array) +## +##print("sensor.shape = ", sensor_1_signal_array.shape) +##sensor_2_signal_array = np.array(sensor_2_signal_array) +## +#half_time = (int)(time_array.size/2) +#pl.plot(time_array, input_signal_array/np.max(input_signal_array[half_time:])) +#pl.plot(time_array, +# sensor_1_signal_array/np.max(sensor_1_signal_array[half_time:]) +# ) +##pl.plot(time_array, +## sensor_2_signal_array/np.max(sensor_2_signal_array[half_time:]) +## ) +#pl.xlabel(r"Time (ps)") +#pl.ylim([-1.1, 1.1]) +#pl.legend(['Input', 'Sensor 1', 'Sensor 2']) +#pl.savefig('paper_plots/IV_55_Ghz_tau_ee_0.2_tau_D_5.png') +##half_time = 0 +#input_normalized = input_signal_array[half_time:]/np.max(input_signal_array[half_time:]) +#sensor_normalized = sensor_1_signal_array[half_time:]/np.max(sensor_1_signal_array[half_time:]) +# +## Calculate the phase diff. Copied from: +## https://stackoverflow.com/questions/6157791/find-phase-difference-between-two-inharmonic-waves +#corr = correlate(input_normalized, sensor_normalized) +#nsamples = input_normalized.size +#time_corr = time_array[half_time:] +#dt_corr = np.linspace(-time_corr[-1] + time_corr[0], +# time_corr[-1] - time_corr[0], 2*nsamples-1) +#time_shift = dt_corr[corr.argmax()] +# +## force the phase shift to be in [-pi:pi] +#period = 1./AC_freq +#phase_diff = 2*np.pi*(((0.5 + time_shift/period) % 1.0) - 0.5) +#print("density.shape = ", density.shape) +#print("Phase diff = ", phase_diff) + +#phase_vs_x = [] +#for i in range(sensor_1_signal_array[0, :].size): +# signal = sensor_1_signal_array[:, i] +# corr = correlate(input_signal_array, signal) +# nsamples = input_signal_array.size +# half_time = 0 +# time_corr = time_array[half_time:] +# dt_corr = np.linspace(-time_corr[-1] + time_corr[0], +# time_corr[-1] - time_corr[0], +# 2*nsamples-1 +# ) +# time_shift = dt_corr[corr.argmax()] +# +# # force the phase shift to be in [-pi:pi] +# period = 1./params.AC_freq +# phase_diff = 2*np.pi*(((0.5 + time_shift/period) % 1.0) - 0.5) +# phase_vs_x.append(phase_diff) +# +#phase_vs_x = np.array(phase_vs_x) +#print("phase_vs_x.shape = ", phase_vs_x.shape) +#np.savetxt("paper_plots/phase_vs_x_tau_ee_0.2_tau_D_1.txt", phase_vs_x) +#np.savetxt("paper_plots/q2_tau_ee_0.2_tau_D_1.txt", q2) +#print("density.shape = ", density.shape) + +