SNES Poisson solver [WIP]

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 \

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,40 +149,55 @@ 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)
 
         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)

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)

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