Skip to content
Snippets Groups Projects
Commit a4046513 authored by Stefan Possanner's avatar Stefan Possanner
Browse files

Merge branch 'fix.full_3d_eval' into 'master' 

more tests to evaluate  at 3d arrays which are not meshgrid

See merge request spossann/gvec_to_python!9
parents f22d5e56 8a4048c0
No related branches found
No related tags found
1 merge request!9more tests to evaluate at 3d arrays which are not meshgrid
Pipeline #155150 passed
Stage: test
Show whitespace changes
Inline Side-by-side
src/gvec_to_python/geometry/domain.py
+ 87
− 89
View file @ a4046513
......@@ -172,11 +172,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays
'''
return self.hmap(*self.Xmap(s, th, ze))
......@@ -186,11 +186,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation),
A 2d numpy array (single-point evaluation) or a 5d numpy array,
with the last two indices referring to the Jacobian entries at each point.'''
tmp1 = self.dX(s, th, ze) # Jacobian of the Xmap
......@@ -211,7 +211,7 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -227,11 +227,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.df_gvec(s, th, ze))
......@@ -241,11 +241,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
tmp = np.ascontiguousarray(self.df_gvec(s, th, ze)) # Jacobian of f_gvec
......@@ -264,11 +264,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.g_gvec(s, th, ze))
......@@ -354,7 +354,7 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -368,11 +368,11 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation),
A 2d numpy array (single-point evaluation) or a 5d numpy array,
with the last two indices referring to the Jacobian entries at each point.'''
tmp1 = self.dP(s2, th2, ze2) # Jacobian of the Pmap
......@@ -393,7 +393,7 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -409,11 +409,11 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.df_pest(s2, th2, ze2))
......@@ -423,11 +423,11 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
tmp = np.ascontiguousarray(self.df_pest(s2, th2, ze2)) # Jacobian of f_gvec
......@@ -446,11 +446,11 @@ class GVEC_domain:
Parameters
----------
s2, th2, ze2 : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.g_pest(s2, th2, ze2))
......@@ -461,7 +461,7 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -475,11 +475,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation),
A 2d numpy array (single-point evaluation) or a 5d numpy array,
with the last two indices referring to the Jacobian entries at each point.'''
tmp1 = self.dL(s, u, v) # Jacobian of the Lmap
......@@ -500,7 +500,7 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -514,11 +514,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.df_unit(s, u, v))
......@@ -528,11 +528,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
tmp = np.ascontiguousarray(self.df_unit(s, u, v)) # Jacobian of f_unit
......@@ -551,11 +551,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.g_unit(s, u, v))
......@@ -567,7 +567,7 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -581,11 +581,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation),
A 2d numpy array (single-point evaluation) or a 5d numpy array,
with the last two indices referring to the Jacobian entries at each point.'''
tmp1 = self.dL(s, u, v) # Jacobian of the Lmap
......@@ -606,7 +606,7 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
......@@ -620,11 +620,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.df_unit_pest(s, u, v))
......@@ -634,11 +634,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
tmp = np.ascontiguousarray(self.df_unit_pest(s, u, v)) # Jacobian of f_unit_pest
......@@ -657,11 +657,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.g_unit_pest(s, u, v))
......@@ -672,11 +672,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.
'''
return self.Xmap(s, th, ze)
......@@ -686,11 +686,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation),
A 2d numpy array (single-point evaluation) or a 5d numpy array,
with the last two indices referring to the Jacobian entries at each point.
'''
return self.dX(s, th, ze)
......@@ -701,11 +701,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A float (single-point evaluation) or a 3d numpy array (meshgrid evaluation).
A float (single-point evaluation) or a 3d numpy array.
'''
return self.det_dX(s, th, ze)
......@@ -715,11 +715,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.df_wo_hmap(s, th, ze))
......@@ -729,11 +729,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
tmp = np.ascontiguousarray(self.df_wo_hmap(s, th, ze)) # Jacobian of f_wo_hmap
......@@ -752,11 +752,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 2d numpy array (single-point evaluation) or a 5d numpy array (meshgrid evaluation).'''
A 2d numpy array (single-point evaluation) or a 5d numpy array.'''
return np.linalg.inv(self.g_wo_hmap(s, th, ze))
......@@ -772,11 +772,11 @@ class GVEC_domain:
Parameters
----------
s, th, ze : float or array-like
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).'''
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.'''
s, th, ze = self.prepare_args(s, th, ze, sparse=False)
......@@ -896,11 +896,11 @@ class GVEC_domain:
Parameters
----------
q1, q2, ze : float or array-like
Coordinates in R^2 x [0, 2*pi]. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in R^2 x [0, 2*pi]. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).'''
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.'''
x = q1 * np.cos(ze)
y = - q1 * np.sin(ze)
......@@ -958,11 +958,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).'''
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.'''
s, u, v = self.prepare_args(s, u, v, sparse=False)
return s, (2*np.pi)*u, (2*np.pi)*self.tor_fraction*v
......@@ -999,38 +999,32 @@ class GVEC_domain:
return np.ones_like(s) * (2*np.pi)**2 * self.tor_fraction
# pest map
def Pmap(self, s2, th2, ze2):
'''The map (s, th, ze) = P(s2, th2, ze2) giving the gvec coordinates (s, th, ze) in [0, 1] x [0, 2*pi]^2
in terms of straight-field-line (PEST) coordinates (s2, th2, ze2) in [0, 1] x [0, 2*pi]^2.
def Pmap(self, s_in, th2_in, ze_in):
'''The map (s, th, ze) = P(s, th2, ze) giving the gvec coordinates (s, th, ze) in [0, 1] x [0, 2*pi]^2
in terms of straight-field-line (PEST) coordinates (s, th2, ze) in [0, 1] x [0, 2*pi]^2.
Parameters
----------
s2, th2, ze2 : float or array-like
PEST coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d or 3d (from dense meshgrid).
s_in, th2_in, ze_in : float or array-like
PEST coordinates in [0, 1] x [0, 2*pi]^2. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).'''
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.'''
s, th2, ze = self.prepare_args(s2, th2, ze2, sparse=False)
def find_th(s_val,th2_val,ze_val):
''' 1D root finder solving th+lambda(s,th,ze)-th2=0 for th, given s,th2,ze '''
func = lambda th: th + self.LA_base.eval_stz(s_val, th, ze_val, self.LA_coef) - th2_val
fprime = lambda th: 1 + self.LA_base.eval_stz(s_val, th, ze_val, self.LA_coef, der='th')
return newton(func, th2_val, fprime) #th2_val as startvalue
if isinstance(s, float):
assert isinstance(th2, float)
assert isinstance(ze, float)
func = lambda th: th + self.LA_base.eval_stz(s, th, ze, self.LA_coef) - th2
fprime = lambda th: 1 + self.LA_base.eval_stz(s, th, ze, self.LA_coef, der='th')
th = newton(func, th2, fprime)
else:
th = np.empty((s + th2 + ze).shape) # broadcast shape if sparse meshgrid
for i in range(s.shape[0]):
for j in range(th2.shape[1]):
for k in range(ze.shape[2]):
func = lambda th: th + self.LA_base.eval_stz(s[i, 0, 0], th, ze[0, 0, k], self.LA_coef) - th2[i, j, k]
fprime = lambda th: 1 + self.LA_base.eval_stz(s[i, 0, 0], th, ze[0, 0, k], self.LA_coef, der='th')
th[i, j, k] = newton(func, th2[i, j, k], fprime)
s, th2, ze = self.prepare_args(s_in, th2_in, ze_in, sparse=True)
vfind_th=np.vectorize(find_th, otypes=[float])
th=vfind_th(s,th2,ze)
return s + 0*th, th, ze + 0*th
def dP(self, s2, th2, ze2):
'''Jacobian of the Pmap.'''
......@@ -1068,11 +1062,11 @@ class GVEC_domain:
Parameters
----------
s, u, v : float or array-like
Coordinates in [0, 1]^3. If list or np.array, must be either 1d or 3d (from meshgrid).
Coordinates in [0, 1]^3. If list or np.array, must be either 1d (meshgrid evaluation) or 3d numpy arrays.
Returns
-------
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).
A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays.
'''
return self.Pmap(*self.Lmap(s, u, v))
......@@ -1100,7 +1094,7 @@ class GVEC_domain:
return np.linalg.det(self.dPL(s, u, v))
@staticmethod
def prepare_args(s, u, v, sparse=True):
def prepare_args(s, u, v, sparse=False):
'''Checks consistency, prepares arguments and performs meshgrid if necessary.'''
if isinstance(s, (list, tuple)):
......@@ -1119,12 +1113,16 @@ class GVEC_domain:
assert isinstance(
v, np.ndarray), '3rd argument should be of type `np.ndarray`. Got {} instead.'.format(type(v))
if s.ndim == 1:
assert s.ndim == u.ndim, '2nd argument has different dimensions than the 1st. Expected {}, got {} instead.'.format(
s.ndim, u.ndim)
assert s.ndim == v.ndim, '3rd argument has different dimensions than the 1st. Expected {}, got {} instead.'.format(
s.ndim, v.ndim)
if s.ndim == 1:
s, u, v = np.meshgrid(s, u, v, indexing='ij', sparse=sparse)
elif s.ndim == 3:
is_meshgrid=(np.all(np.minimum.reduce(s,axis=(1,2))==np.maximum.reduce(s,axis=(1,2))) and
np.all(np.minimum.reduce(u,axis=(0,2))==np.maximum.reduce(u,axis=(0,2))) and
np.all(np.minimum.reduce(v,axis=(0,1))==np.maximum.reduce(v,axis=(0,1))))
return s, u, v
......
tests/test_gvec_domain.py
+ 15
− 6
View file @ a4046513
......@@ -110,7 +110,7 @@ def test_values(use_pyccel,gvec_file,unit_tor_domain):
s = s0[0:n0]
th = th0[0:n1]
ze = ze0[0:n2]
s, th, ze = gvec.domain.prepare_args(s, th, ze)
s, th, ze = gvec.domain.prepare_args(s, th, ze,sparse=True)
th2 = th + gvec._lambda(s, th, ze)
......@@ -129,11 +129,20 @@ def test_values(use_pyccel,gvec_file,unit_tor_domain):
assert np.allclose(gvec.domain.f_gvec(s, th, ze)[1], gvec.domain.f_unit(s, th/(2*np.pi), ze/(2*np.pi*gvec.domain.tor_fraction))[1])
assert np.allclose(gvec.domain.f_gvec(s, th, ze)[2], gvec.domain.f_unit(s, th/(2*np.pi), ze/(2*np.pi*gvec.domain.tor_fraction))[2])
# check Pmap, computing theta with 3d arrays of non-meshgrid coordinates
s =np.random.rand(2,5,3)
th=np.random.rand(2,5,3)*2.*np.pi
ze=np.random.rand(2,5,3)*np.pi
th2 = th + gvec._lambda(s, th, ze)
assert np.allclose(s, gvec.domain.Pmap(s, th2, ze)[0])
assert np.allclose(th, gvec.domain.Pmap(s, th2, ze)[1])
assert np.allclose(ze, gvec.domain.Pmap(s, th2, ze)[2])
# compute theta with equidistant in PEST, compare f and f_pest again
s2=s0
ze2=ze0
th2=np.linspace(0.,1.,6)*2*np.pi
# compute theta with 3d arrays of non-meshgrid coordinates, compare f and f_pest again
s2 = np.random.rand(2,5,3)
th2= np.random.rand(2,5,3)*2*np.pi
ze2= np.random.rand(2,5,3)*np.pi
s,th,ze=gvec.domain.Pmap(s2, th2, ze2)
assert np.allclose(gvec.domain.f_gvec(s, th, ze)[0], gvec.domain.f_pest(s2, th2, ze2)[0])
......
tests/test_gvec_equil.py
+ 81
− 2
View file @ a4046513
......@@ -135,6 +135,85 @@ def test_mhd_fields(mapping, use_pyccel,gvec_file):
assert np.all([tmp[i][j].shape == shp3 for j in range(3) for i in range(2)])
@pytest.mark.parametrize('use_pyccel,gvec_file', [(True ,"testcases/ellipstell/newBC_E4D6_M6N6/GVEC_ELLIPSTELL_E4D6_M6N6_State_0001_00200000")])
@pytest.mark.parametrize('unit_tor_domain', ["one-fp","full"])
def test_values(use_pyccel,gvec_file,unit_tor_domain):
'''Compare values of b_cart for different mappings provided by GVEC_domain.'''
import os
import numpy as np
from gvec_to_python.reader.gvec_reader import create_GVEC_json
from gvec_to_python import GVEC
import gvec_to_python as _
# give absolute paths to the files
pkg_path = _.__path__[0]
dat_file_in = os.path.join( pkg_path, gvec_file+'.dat')
json_file_out = os.path.join( pkg_path, gvec_file+'.json')
create_GVEC_json(dat_file_in, json_file_out)
# main object
gvec = GVEC(json_file_out, use_pyccel=use_pyccel,unit_tor_domain=unit_tor_domain)
# compare f and f_pest
s0 = np.array([.5,.99])
th0 = np.random.rand(5)*2*np.pi
ze0 = np.random.rand(3)*np.pi/gvec.nfp
# all combinations of 1D arrays
for n0,n1,n2 in [(-1,-1,-1),(1,1,1),(2,1,1),(2,5,3),(2,5,1),(2,1,3),(1,5,3),(1,5,1),(1,1,3)]:
if((n0,n1,n2)==(-1,-1,-1)): # also test a fully 3d array
s =np.random.rand(3,4,2)
th=np.random.rand(3,4,2)*2*np.pi
ze=np.random.rand(3,4,2)*np.pi/gvec.nfp
else:
s = s0[0:n0]
th = th0[0:n1]
ze = ze0[0:n2]
s, th, ze = gvec.domain.prepare_args(s, th, ze,sparse=False)
th2 = th + gvec._lambda(s, th, ze)
gvec.mapping="gvec"
b_gvec, xyz_gvec = gvec.b_cart(s,th,ze)
a_gvec, xyz_gvec = gvec.a_cart(s,th,ze)
# compare f and f_pest
gvec.mapping="pest"
b_pest, xyz_pest = gvec.b_cart(s,th2,ze)
a_pest, xyz_pest = gvec.a_cart(s,th2,ze)
assert np.allclose(b_gvec[0], b_pest[0])
assert np.allclose(b_gvec[1], b_pest[1])
assert np.allclose(b_gvec[2], b_pest[2])
assert np.allclose(a_gvec[0], a_pest[0])
assert np.allclose(a_gvec[1], a_pest[1])
assert np.allclose(a_gvec[2], a_pest[2])
# compare f and f_unit, where the tor_fraction is needed for unit to have the same point positions
gvec.mapping="unit"
b_unit, xyz_unit = gvec.b_cart(s,th/(2*np.pi),ze/(2*np.pi*gvec.domain.tor_fraction))
a_unit, xyz_unit = gvec.a_cart(s,th/(2*np.pi),ze/(2*np.pi*gvec.domain.tor_fraction))
assert np.allclose(b_gvec[0], b_unit[0])
assert np.allclose(b_gvec[1], b_unit[1])
assert np.allclose(b_gvec[2], b_unit[2])
assert np.allclose(a_gvec[0], a_unit[0])
assert np.allclose(a_gvec[1], a_unit[1])
assert np.allclose(a_gvec[2], a_unit[2])
# compare f_unit_pest and f_gvec
gvec.mapping="unit_pest"
b_unit_pest, xyz_unit_pest = gvec.b_cart(s,th2/(2*np.pi),ze/(2*np.pi*gvec.domain.tor_fraction))
a_unit_pest, xyz_unit_pest = gvec.a_cart(s,th2/(2*np.pi),ze/(2*np.pi*gvec.domain.tor_fraction))
assert np.allclose(b_unit_pest[0], b_gvec[0])
assert np.allclose(b_unit_pest[1], b_gvec[1])
assert np.allclose(b_unit_pest[2], b_gvec[2])
assert np.allclose(a_unit_pest[0], a_gvec[0])
assert np.allclose(a_unit_pest[1], a_gvec[1])
assert np.allclose(a_unit_pest[2], a_gvec[2])
if __name__ == '__main__':
test_mhd_fields('gvec', use_pyccel=False)
test_mhd_fields('pest', use_pyccel=True)
test_mhd_fields('gvec', use_pyccel=False,gvec_file='testcases/ellipstell/newBC_E1D6_M6N6/GVEC_ELLIPSTELL_E1D6_M6N6_State_0000_00200000')
test_mhd_fields('pest', use_pyccel=True ,gvec_file='testcases/ellipstell/newBC_E4D6_M6N6/GVEC_ELLIPSTELL_E4D6_M6N6_State_0001_00200000')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment