Commit bfe4fe54 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'fix-fd-operator-padding' into 'master'

Use correct padding for finite-difference operator

Closes #232

See merge request !600
parents 718c0087 7a1f3126
Pipeline #98865814 passed with stage
in 3 minutes and 25 seconds
......@@ -45,21 +45,18 @@ class FDOperator:
if sum([offset != 0 for offset in offset_c]) >= 2:
cfd = False
maxoffset_c = [max([offset_c[c] for offset_c in offset_pc])
for c in range(3)]
mp = maxoffset_c[0]
if maxoffset_c[1] != mp or maxoffset_c[2] != mp:
mp = max(maxoffset_c)
mp = np.abs(offset_pc).max() # padding
n_c = gd.n_c
M_c = n_c + 2 * mp
stride_c = np.array([M_c[1] * M_c[2], M_c[2], 1])
offset_p = np.dot(offset_pc, stride_c)
coef_p = np.ascontiguousarray(coef_p, float)
neighbor_cd = gd.neighbor_cd
assert coef_p.ndim == 1
assert coef_p.shape == offset_p.shape
assert dtype in [float, complex]
self.dtype = dtype
self.shape = tuple(n_c)
......@@ -69,7 +66,7 @@ class FDOperator:
comm = None
assert neighbor_cd.flags.c_contiguous and offset_p.flags.c_contiguous
self.mp = mp # padding
self.mp = mp
self.gd = gd
self.npoints = len(coef_p)
self.coef_p = coef_p
......
from __future__ import print_function
from gpaw.fd_operators import Gradient
from gpaw.fd_operators import Gradient, FDOperator
import numpy as np
from gpaw.grid_descriptor import GridDescriptor
from gpaw.mpi import world
......@@ -23,13 +22,15 @@ print(dadx.itemsize, dadx.dtype, dadx.shape)
gradx.apply(a, dadx)
# a = [ 0. 1. 2. 3. 4. 5. 6. 7.]
#
# da
# -- = [-2.5 1. 1. 1. 1. 1. 1. -2.5]
# dx
dadx = gd.collect(dadx, broadcast=True)
assert dadx[3, 0, 0] == 1.0 and np.sum(dadx[:, 0, 0]) == 0.0
dAdx = gd.collect(dadx, broadcast=True)
assert not (dAdx[:, 0, 0] - [-3, 1, 1, 1, 1, 1, 1, -3]).any()
# Backwards FD operator:
gradx2 = FDOperator([-1, 1], [[-1, 0, 0], [0, 0, 0]], gd)
gradx2.apply(a, dadx)
dAdx = gd.collect(dadx, broadcast=True)
assert not (dAdx[:, 0, 0] - [-7, 1, 1, 1, 1, 1, 1, 1]).any()
gd = GridDescriptor((1, 8, 1), (1.0, 8.0, 1.0), (1, 0, 1), comm=domain_comm)
dady = gd.zeros()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment