diff --git a/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py b/bolt/lib/nonlinear_solver/EM_fields_solver/electrostatic.py index 3eb651c2..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') @@ -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)) 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 4cf9b2a2..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,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() diff --git a/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py b/bolt/lib/nonlinear_solver/tests/test_fft_poisson.py index a9605ff0..2efc4bf7 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 = 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() \ No newline at end of file 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