From c6f751d83db5d7311113b7f06c8dcdec2a384228 Mon Sep 17 00:00:00 2001 From: Shyam Sundar Sankaran Date: Fri, 22 Sep 2017 11:08:07 +0530 Subject: [PATCH 1/2] WIP Poisson solver debugging --- .../EM_fields_solver/electrostatic.py | 2 +- .../test_compute_electrostatic_fields.py | 55 +++++++++++- .../tests/test_fft_poisson.py | 86 ++++++++++++++++--- 3 files changed, 127 insertions(+), 16 deletions(-) diff --git a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py index 3eb651c2..b23ac42d 100644 --- a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py +++ b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py @@ -152,7 +152,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) 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 4cf9b2a2..e38b1198 100644 --- a/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py +++ b/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py @@ -18,6 +18,20 @@ 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 = 20 # 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 +91,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 +131,38 @@ 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(5, 10) + + for i in range(N.size): + obj = test(N[i]) + compute_electrostatic_fields(obj) + + E1_expected = 0 + + E2_expected = -0.5/20 * ( af.log(af.cosh(( obj.q2 - q2_minus)*20)) + - af.log(af.cosh(( obj.q2 - q2_plus )*20)) + ) + + 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) diff --git a/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py b/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py index a9605ff0..30098b5a 100644 --- a/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py +++ b/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 = 20 # 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): @@ -85,26 +107,68 @@ 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 - + compute_moments = compute_moments_hyperbolic def test_fft_poisson(): obj = test() 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/20 * ( af.log(af.cosh(( obj.q2 - 0.25)*20)) + - af.log(af.cosh(( obj.q2 - 0.75)*20)) + ) + + # 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.contourf(-np.array(obj.compute_moments('density') - 1)[N_g:-N_g, N_g:-N_g], 100) + pl.colorbar() + 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() \ No newline at end of file From 32a6e8d34df4e98dd1ea8ff21192b513c06be838 Mon Sep 17 00:00:00 2001 From: Shyam Sundar Sankaran Date: Tue, 26 Sep 2017 17:59:29 +0530 Subject: [PATCH 2/2] Changes made during debugging --- .../EM_fields_solver/electrostatic.py | 9 ++-- bolt/lib/nonlinear_solver/nonlinear_solver.py | 2 + .../test_compute_electrostatic_fields.py | 53 ++++++++++--------- .../tests/test_fft_poisson.py | 39 +++++++------- .../1D1V/fields_collisionless/main.py | 3 ++ 5 files changed, 59 insertions(+), 47 deletions(-) diff --git a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py index b23ac42d..3bdb9b4f 100644 --- a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py +++ b/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') @@ -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)) diff --git a/bolt/lib/nonlinear_solver/nonlinear_solver.py b/bolt/lib/nonlinear_solver/nonlinear_solver.py index 677ff0ed..618da49b 100644 --- a/bolt/lib/nonlinear_solver/nonlinear_solver.py +++ b/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 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 e38b1198..b59405d7 100644 --- a/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py +++ b/bolt/lib/nonlinear_solver/tests/test_compute_electrostatic_fields.py @@ -13,6 +13,7 @@ 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 \ @@ -25,7 +26,7 @@ def compute_moments_gaussian(self, *args): q2_minus = 0.25 q2_plus = 0.75 - regulator = 20 # larger value makes the transition sharper + 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) @@ -132,37 +133,39 @@ def test_compute_electrostatic_fields_1(): assert (abs(poly_E1[0] + 2) < 0.2) assert (abs(poly_E2[0] + 2) < 0.2) -def test_compute_electrostatic_fields_2(): +# def test_compute_electrostatic_fields_2(): - error_E1 = np.zeros(5) - error_E2 = np.zeros(5) +# error_E1 = np.zeros(5) +# error_E2 = np.zeros(5) - N = 2**np.arange(5, 10) +# N = 2**np.arange(10, 11) - for i in range(N.size): - obj = test(N[i]) - compute_electrostatic_fields(obj) +# for i in range(N.size): +# obj = test(N[i]) +# compute_electrostatic_fields(obj) - E1_expected = 0 +# E1_expected = 0 - E2_expected = -0.5/20 * ( af.log(af.cosh(( obj.q2 - q2_minus)*20)) - - af.log(af.cosh(( obj.q2 - q2_plus )*20)) - ) +# 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 +# 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_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()) +# 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) +# 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() diff --git a/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py b/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py index 30098b5a..2efc4bf7 100644 --- a/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py +++ b/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py @@ -26,7 +26,7 @@ def compute_moments_hyperbolic(self, *args): q2_minus = 0.25 q2_plus = 0.75 - regulator = 20 # larger value makes the transition sharper + 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) @@ -56,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 @@ -107,11 +107,13 @@ 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) - _communicate_fields = communicate_fields - compute_moments = compute_moments_hyperbolic + _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 @@ -124,8 +126,8 @@ def test_fft_poisson(): E1_expected = 0 - E2_expected = -0.5/20 * ( af.log(af.cosh(( obj.q2 - 0.25)*20)) - - af.log(af.cosh(( obj.q2 - 0.75)*20)) + 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 @@ -140,22 +142,21 @@ def test_fft_poisson(): N_g = obj.N_ghost - pl.contourf(-np.array(obj.compute_moments('density') - 1)[N_g:-N_g, N_g:-N_g], 100) - pl.colorbar() + 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.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.colorbar() + # pl.show() + # pl.clf() pl.contourf(np.array(E2_expected)[N_g:-N_g, N_g:-N_g], 100) pl.colorbar() diff --git a/example_problems/nonrelativistic_boltzmann/linear_modes_damping/perturbed_density/1D1V/fields_collisionless/main.py b/example_problems/nonrelativistic_boltzmann/linear_modes_damping/perturbed_density/1D1V/fields_collisionless/main.py index 241b79eb..3370a38e 100644 --- a/example_problems/nonrelativistic_boltzmann/linear_modes_damping/perturbed_density/1D1V/fields_collisionless/main.py +++ b/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