Poisson solver debugging [WIP]

bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py

@@ -1,6 +1,8 @@
 #!/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
@@ -74,7 +76,8 @@ def compute_electrostatic_fields(self, performance_test_flag = False):
     ksp = PETSc.KSP().create()
 
     ksp.setOperators(A)
-    ksp.setType('cg')
+    ksp.setFromOptions()
+    # ksp.setType('cg')
 
     pc = ksp.getPC()
     pc.setType('none')
@@ -152,7 +155,7 @@ def fft_poisson(self, performance_test_flag = False):
     else:
         N_g = self.N_ghost
         rho =   self.physical_system.params.charge_electron \
-              * self.compute_moments('density')[N_g:-N_g, N_g:-N_g]
+              * (self.compute_moments('density')[N_g:-N_g, N_g:-N_g])
 
         k_q1 = fftfreq(rho.shape[0], self.dq1)
         k_q2 = fftfreq(rho.shape[1], self.dq2)
@@ -167,8 +170,8 @@ def fft_poisson(self, performance_test_flag = False):
         potential_hat       = rho_hat / (4 * np.pi**2 * (k_q1**2 + k_q2**2))
         potential_hat[0, 0] = 0
 
-        E1_hat = -1j * 2 * np.pi * (k_q1) * potential_hat
-        E2_hat = -1j * 2 * np.pi * (k_q2) * potential_hat
+        E1_hat = -1j * 2 * np.pi * k_q1 * potential_hat
+        E2_hat = -1j * 2 * np.pi * k_q2 * potential_hat
 
         self.E1[N_g:-N_g, N_g:-N_g] = af.real(af.ifft2(E1_hat))
         self.E2[N_g:-N_g, N_g:-N_g] = af.real(af.ifft2(E2_hat))

bolt/lib/nonlinear_solver/nonlinear_solver.py

@@ -10,6 +10,8 @@
 # Importing dependencies:
 import arrayfire as af
 import numpy as np
+import petsc4py, sys
+petsc4py.init(sys.argv)
 from mpi4py import MPI
 from petsc4py import PETSc
 from prettytable import PrettyTable

bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py

@@ -13,11 +13,26 @@
 
 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
 
+def compute_moments_sinusoidal(self, *args):
+    return (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
+
+    regulator = 40  # larger value makes the transition sharper
+
+    rho = 1 + 0.5 * (  af.tanh(( self.q2 - q2_minus)*regulator) 
+                     - af.tanh(( self.q2 - q2_plus )*regulator)
+                    )
+
+    return(rho)
 
 class test(object):
     def __init__(self, N):
@@ -77,11 +92,9 @@ def __init__(self, N):
                                             stencil_type=1,
                                           )
 
-    def compute_moments(self, *args):
-        return (1 + af.sin(2 * np.pi * self.q1 + 4 * np.pi * self.q2))
-
+    compute_moments = compute_moments_sinusoidal
 
-def test_compute_electrostatic_fields():
+def test_compute_electrostatic_fields_1():
 
     error_E1 = np.zeros(5)
     error_E2 = np.zeros(5)
@@ -119,3 +132,40 @@ def test_compute_electrostatic_fields():
 
     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)
+
+#     assert (abs(poly_E1[0] + 2) < 0.2)
+#     assert (abs(poly_E2[0] + 2) < 0.2)
+
+test_compute_electrostatic_fields_1()

bolt/lib/nonlinear_solver/tests/test_fft_poisson.py

@@ -12,11 +12,33 @@
 
 import numpy as np
 import arrayfire as af
+af.set_backend('cpu')
+import pylab as pl
 from petsc4py import PETSc
 
 from bolt.lib.nonlinear_solver.EM_fields_solver.electrostatic import fft_poisson
 from bolt.lib.nonlinear_solver.communicate import communicate_fields
 
+def compute_moments_sinusoidal(self, *args):
+    return (1 + af.sin(2 * np.pi * self.q1 + 4 * np.pi * self.q2))
+
+def compute_moments_hyperbolic(self, *args):
+    q2_minus = 0.25
+    q2_plus  = 0.75
+
+    regulator = 80 # larger value makes the transition sharper
+
+    rho = 1 + 0.5 * (  af.tanh(( self.q2 - q2_minus)*regulator) 
+                     - af.tanh(( self.q2 - q2_plus )*regulator)
+                    )
+
+    rho[:self.N_ghost]  = rho[-2*self.N_ghost:-self.N_ghost]
+    rho[-self.N_ghost:] = rho[self.N_ghost:2*self.N_ghost]
+
+    rho[:, :self.N_ghost]  = rho[:, -2*self.N_ghost:-self.N_ghost]
+    rho[:, -self.N_ghost:] = rho[:, self.N_ghost:2*self.N_ghost]
+
+    return(rho)
 
 class test(object):
     def __init__(self):
@@ -34,8 +56,8 @@ def __init__(self):
         self.q1_end = 1
         self.q2_end = 1
 
-        self.N_q1 = np.random.randint(24, 48)
-        self.N_q2 = np.random.randint(24, 48)
+        self.N_q1 = 1024
+        self.N_q2 = 1024
 
         self.dq1 = (self.q1_end - self.q1_start) / self.N_q1
         self.dq2 = (self.q2_end - self.q2_start) / self.N_q2
@@ -85,26 +107,69 @@ def __init__(self):
         self._local_value_fields = self._da_fields.getVecArray(self._local_fields)
         self._glob_value_fields  = self._da_fields.getVecArray(self._glob_fields)
 
-
-    def compute_moments(self, *args):
-        return (af.sin(2 * np.pi * self.q1 + 4 * np.pi * self.q2))
-
-    _communicate_fields = communicate_fields
-
+    _communicate_fields        = communicate_fields
+    compute_moments_hyperbolic = compute_moments_hyperbolic
+    compute_moments_sinusoidal = compute_moments_sinusoidal
 
 def test_fft_poisson():
     obj = test()
+    obj.compute_moments = obj.compute_moments_hyperbolic
     fft_poisson(obj)
 
-    E1_expected = (0.1 / np.pi) * af.cos(  2 * np.pi * obj.q1
-                                         + 4 * np.pi * obj.q2
-                                        )
+    # 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
+    #                                     )
+
+    E1_expected = 0
 
-    E2_expected = (0.2 / np.pi) * af.cos(  2 * np.pi * obj.q1
-                                         + 4 * np.pi * obj.q2
-                                        )
+    E2_expected = -0.5/80 * (   af.log(af.cosh(( obj.q2 - 0.25)*80)) 
+                              - af.log(af.cosh(( obj.q2 - 0.75)*80))
+                            ) 
+
+    # E1_expected = 0
+
+    # E2_expected = -0.5/5 * 0.5 * (  af.exp(-10*( obj.q2 - 0.25)**2) 
+    #                               - af.exp(-10*( obj.q2 - 0.75)**2)
+    #                              )
+
+    # E1_expected = -0.5 * af.log(1+obj.q1+obj.q2)
+
+    # E2_expected = -0.5 * af.log(1+obj.q1+obj.q2) 
+
+    N_g = obj.N_ghost
+
+    pl.plot((np.array(obj.compute_moments('density') - 1)[N_g:-N_g, N_g:-N_g]).T)
+    pl.show()
+    pl.clf()
+
+    # pl.contourf(np.array(af.convolve2_separable(af.Array([0, 1, 0]), 
+    #                                             af.Array([1, 0, -1]),
+    #                                             E2_expected
+    #                                            )
+    #                     )[N_g:-N_g, N_g:-N_g]/(2*obj.dq2), 
+    #             100
+    #            )
+
+    # pl.colorbar()
+    # pl.show()
+    # pl.clf()
+
+    pl.contourf(np.array(E2_expected)[N_g:-N_g, N_g:-N_g], 100)
+    pl.colorbar()
+    pl.show()
+    pl.clf()
+
+    pl.contourf(np.array(obj.E2)[N_g:-N_g, N_g:-N_g], 100)
+    pl.colorbar()
+    pl.show()
 
     error_E1 = af.sum(af.abs(obj.E1 - E1_expected)) / (obj.E1.elements())
     error_E2 = af.sum(af.abs(obj.E2 - E2_expected)) / (obj.E2.elements())
 
     assert (error_E1 < 1e-14 and error_E2 < 1e-14)
+
+test_fft_poisson()

example_problems/nonrelativistic_boltzmann/linear_modes_damping/perturbed_density/1D1V/fields_collisionless/main.py

@@ -1,6 +1,9 @@
 import arrayfire as af
 import numpy as np
 import pylab as pl
+import petsc4py
+import sys
+petsc4py.init(sys.argv)
 
 from bolt.lib.physical_system import physical_system