diff --git a/cuda/abs.cu b/cuda/abs.cu new file mode 100644 index 000000000..428d1822c --- /dev/null +++ b/cuda/abs.cu @@ -0,0 +1,10 @@ +// dst[i] = abs(a[i]) + +extern "C" __global__ void +unary_abs(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = fabsf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/acos.cu b/cuda/acos.cu new file mode 100644 index 000000000..d299a0048 --- /dev/null +++ b/cuda/acos.cu @@ -0,0 +1,16 @@ +// dst[i] = acos(a[i]), returns 0 for a[i] outside of domain [-1, 1] + +extern "C" __global__ void +unary_acos(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + float x = a[i]; + if (x >= -1.0f && x <= 1.0f){ + dst[i] = acosf(x); + } + else{ + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/acosh.cu b/cuda/acosh.cu new file mode 100644 index 000000000..db1ee9644 --- /dev/null +++ b/cuda/acosh.cu @@ -0,0 +1,15 @@ +//dst[i] = acosh(a[i]), returns 0 for a[i] outside of domain [1, inf) + +extern "C" __global__ void +unary_acosh(float* __restrict__ dst, float* __restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + if (x >= 1.0f) { + dst[i] = acoshf(x); + } else { + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/asin.cu b/cuda/asin.cu new file mode 100644 index 000000000..dcaefb53d --- /dev/null +++ b/cuda/asin.cu @@ -0,0 +1,15 @@ +// dst[i] = asin(a[i]), returns 0 for a[i] outside of domain [-1, 1] + +extern "C" __global__ void +unary_asin(float* __restrict__ dst, float* __restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + if (x >= -1.0f && x <= 1.0f) { + dst[i] = asinf(x); + } else { + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/asinh.cu b/cuda/asinh.cu new file mode 100644 index 000000000..9c05a9e1b --- /dev/null +++ b/cuda/asinh.cu @@ -0,0 +1,10 @@ +//dst[i] = asinh(a[i]) + +extern "C" __global__ void +unary_asinh(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = asinhf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/atan.cu b/cuda/atan.cu new file mode 100644 index 000000000..7f1e99e96 --- /dev/null +++ b/cuda/atan.cu @@ -0,0 +1,10 @@ +//dst[i] = atan(a[i]) + +extern "C" __global__ void +unary_atan(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = atanf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/atan2.cu b/cuda/atan2.cu new file mode 100644 index 000000000..6bbc54c7c --- /dev/null +++ b/cuda/atan2.cu @@ -0,0 +1,10 @@ +//dst[i] = atan2(a[i], b[i]) + +extern "C" __global__ void +pw_atan2(float *__restrict__ dst, float *__restrict__ a, float *__restrict__ b, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = atan2f(a[i], b[i]); + } +} \ No newline at end of file diff --git a/cuda/atanh.cu b/cuda/atanh.cu new file mode 100644 index 000000000..bd2893c20 --- /dev/null +++ b/cuda/atanh.cu @@ -0,0 +1,16 @@ +// dst[i] = atanh(a[i]), returns 0 for a[i] outside of domain (-1, 1) + +extern "C" __global__ void +unary_atanh(float* __restrict__ dst, float* __restrict__ a, int N) { + + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + if (x > -1.0f && x < 1.0f) { + dst[i] = atanhf(x); + } else { + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/cos.cu b/cuda/cos.cu new file mode 100644 index 000000000..24530354d --- /dev/null +++ b/cuda/cos.cu @@ -0,0 +1,9 @@ +//dst[i] = cos(a[i]) +extern "C" __global__ void +unary_cos(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = cosf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/cosh.cu b/cuda/cosh.cu new file mode 100644 index 000000000..6f4f22603 --- /dev/null +++ b/cuda/cosh.cu @@ -0,0 +1,9 @@ +//dst[i] = cosh(a[i]) +extern "C" __global__ void +unary_cosh(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = coshf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/erf.cu b/cuda/erf.cu new file mode 100644 index 000000000..a7bcc1a00 --- /dev/null +++ b/cuda/erf.cu @@ -0,0 +1,9 @@ +//dst[i] = erf(a[i]) +extern "C" __global__ void +unary_erf(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = erff(a[i]); + } +} \ No newline at end of file diff --git a/cuda/erfc.cu b/cuda/erfc.cu new file mode 100644 index 000000000..d992a8689 --- /dev/null +++ b/cuda/erfc.cu @@ -0,0 +1,9 @@ +//dst[i] = erfc(a[i]) +extern "C" __global__ void +unary_erfc(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = erfcf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/exp.cu b/cuda/exp.cu new file mode 100644 index 000000000..3deae42c3 --- /dev/null +++ b/cuda/exp.cu @@ -0,0 +1,9 @@ +//dst[i] = exp(a[i]) +extern "C" __global__ void +unary_exp(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = expf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/gamma.cu b/cuda/gamma.cu new file mode 100644 index 000000000..4be7a3764 --- /dev/null +++ b/cuda/gamma.cu @@ -0,0 +1,15 @@ +// dst[i] = gamma(a[i]), returns 0 for x <= 0, returns 0 for a[i] outside of domain (0, inf) with poles at non-positive integers +extern "C" __global__ void +unary_gamma(float* __restrict__ dst, float* __restrict__ a, int N) { + + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + if (x > 0.0f) { + dst[i] = tgammaf(x); + } else { + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/heaviside.cu b/cuda/heaviside.cu new file mode 100644 index 000000000..df157aea6 --- /dev/null +++ b/cuda/heaviside.cu @@ -0,0 +1,19 @@ +//dst[i] = heaviside(a[i]) +// returns 0 for x < 0, 0.5 for x == 0, 1 for x > 0 +extern "C" __global__ void +unary_heaviside(float* __restrict__ dst, float* __restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + + if (x > 0.0f) { + dst[i] = 1.0f; + } else if (x < 0.0f) { + dst[i] = 0.0f; + } else { + // x == 0 + dst[i] = 0.5f; + } + } +} \ No newline at end of file diff --git a/cuda/log.cu b/cuda/log.cu new file mode 100644 index 000000000..d84ce7cb4 --- /dev/null +++ b/cuda/log.cu @@ -0,0 +1,13 @@ +// dst[i] = log(a[i]), returns 0 for non-positive input +extern "C" __global__ void +unary_log(float* __restrict__ dst, float* __restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + if (a[i] > 0.0f) { + dst[i] = logf(a[i]); + } else { + dst[i] = 0.0f; + } + } +} \ No newline at end of file diff --git a/cuda/mod.cu b/cuda/mod.cu new file mode 100644 index 000000000..cf887d6ec --- /dev/null +++ b/cuda/mod.cu @@ -0,0 +1,10 @@ +//dst[i] = mod(a[i], b[i]) +extern "C" __global__ void +pw_mod(float *__restrict__ dst, float *__restrict__ a, float *__restrict__ b, int N) +{ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + if (i < N) + { + dst[i] = fmodf(a[i], b[i]); + } +} \ No newline at end of file diff --git a/cuda/pow.cu b/cuda/pow.cu new file mode 100644 index 000000000..4b164dd62 --- /dev/null +++ b/cuda/pow.cu @@ -0,0 +1,28 @@ +// dst[i] = pow(a[i], b[i]), pow(a,b) for negative a and b returns -pow(-a,b) for fractional exponents, and for 0^0 returns 1 +extern "C" __global__ void +pw_pow(float* __restrict__ dst, float* __restrict__ a, float* __restrict__ b, int N) { +int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + float x = a[i]; + float y = b[i]; + + // Guard: 0^negative → return 0 + if (x == 0.0f && y < 0.0f) { + dst[i] = 0.0f; + return; + } + + // Check if exponent is effectively integer + float y_floor = floorf(y); + bool y_is_int = fabsf(y - y_floor) < 1e-6f; + + if (x < 0.0f && !y_is_int) { + // Fractional exponent → move minus sign outside + dst[i] = -1.0f * powf(fabsf(x), y); + } else { + // Integer exponent or positive base + dst[i] = powf(x, y); + } + } +} \ No newline at end of file diff --git a/cuda/sin.cu b/cuda/sin.cu new file mode 100644 index 000000000..27ac4a552 --- /dev/null +++ b/cuda/sin.cu @@ -0,0 +1,10 @@ +//dst[i] = sin(a[i]) +extern "C" __global__ void +unary_sin(float *__restrict__ dst, float *__restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) + { + dst[i] = sinf(a[i]); + } +} diff --git a/cuda/sinc.cu b/cuda/sinc.cu new file mode 100644 index 000000000..e9201e384 --- /dev/null +++ b/cuda/sinc.cu @@ -0,0 +1,11 @@ +// dst[i] = sinc(a[i]) +extern "C" __global__ void +unary_sinc(float *__restrict__ dst, float *__restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) + { + float x = a[i]; + dst[i] = (x == 0.0f) ? 1.0f : sinf(x) / x; + } +} \ No newline at end of file diff --git a/cuda/sinh.cu b/cuda/sinh.cu new file mode 100644 index 000000000..494f60adf --- /dev/null +++ b/cuda/sinh.cu @@ -0,0 +1,10 @@ +//dst[i] = sinh(a[i]) +extern "C" __global__ void +unary_sinh(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) + { + dst[i] = sinhf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/tan.cu b/cuda/tan.cu new file mode 100644 index 000000000..16f29751f --- /dev/null +++ b/cuda/tan.cu @@ -0,0 +1,9 @@ +//dst[i] = tan(a[i]), for poles at π/2 + nπ, returns 0 +extern "C" __global__ void +unary_tan(float* __restrict__ dst, float* __restrict__ a, int N) { + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N) { + dst[i] = tanf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/tanh.cu b/cuda/tanh.cu new file mode 100644 index 000000000..e0ef43517 --- /dev/null +++ b/cuda/tanh.cu @@ -0,0 +1,9 @@ +//dst[i] = tanh(a[i]) +extern "C" __global__ void +unary_tanh(float *__restrict__ dst, float *__restrict__ a, int N){ + int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x; + + if (i < N){ + dst[i] = tanhf(a[i]); + } +} \ No newline at end of file diff --git a/cuda/unary_math.go b/cuda/unary_math.go new file mode 100644 index 000000000..2289731bd --- /dev/null +++ b/cuda/unary_math.go @@ -0,0 +1,123 @@ +package cuda + +import ( + "unsafe" + + "github.com/mumax/3/data" + "github.com/mumax/3/util" +) + +// applyUnary is a helper that applies a pointwise async kernel +// to every component of dst and a. dst and a must have the same shape. +func applyUnary(dst, a *data.Slice, kernel func(dst, a unsafe.Pointer, N int, cfg *config)) { + N := dst.Len() + cfg := make1DConf(N) + for c := 0; c < dst.NComp(); c++ { + kernel(dst.DevPtr(c), a.DevPtr(c), N, cfg) + } +} + +// for inputs outside of a function's domain, return 0 +// + +// dst[i] = sin(a[i]) +func QSin(dst, a *data.Slice) { applyUnary(dst, a, k_unary_sin_async) } + +// dst[i] = cos(a[i]) +func QCos(dst, a *data.Slice) { applyUnary(dst, a, k_unary_cos_async) } + +// dst[i] = tan(a[i]) +// poles at π/2 + nπ, returns 0 +func QTan(dst, a *data.Slice) { applyUnary(dst, a, k_unary_tan_async) } + +// dst[i] = exp(a[i]) +func QExp(dst, a *data.Slice) { applyUnary(dst, a, k_unary_exp_async) } + +// dst[i] = log(a[i]) +func QLog(dst, a *data.Slice) { applyUnary(dst, a, k_unary_log_async) } + +// dst[i] = |a[i]| +func QAbs(dst, a *data.Slice) { applyUnary(dst, a, k_unary_abs_async) } + +// dst[i] = acos(a[i]) +// returns 0 for a[i] outside of domain [-1, 1] +func QAcos(dst, a *data.Slice) { applyUnary(dst, a, k_unary_acos_async) } + +// dst[i] = acosh(a[i]) +// returns 0 for a[i] outside of domain [1, inf) +func QAcosh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_acosh_async) } + +// dst[i] = asin(a[i]) +// returns 0 for a[i] outside of domain [-1, 1] +func QAsin(dst, a *data.Slice) { applyUnary(dst, a, k_unary_asin_async) } + +// dst[i] = asinh(a[i]) +func QAsinh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_asinh_async) } + +// dst[i] = atan(a[i]) +func QAtan(dst, a *data.Slice) { applyUnary(dst, a, k_unary_atan_async) } + +// dst[i] = atanh(a[i]) +// returns 0 for a[i] outside of domain (-1, 1) +func QAtanh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_atanh_async) } + +// dst[i] = cosh(a[i]) +func QCosh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_cosh_async) } + +// dst[i] = sinh(a[i]) +func QSinh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_sinh_async) } + +// dst[i] = tanh(a[i]) +func QTanh(dst, a *data.Slice) { applyUnary(dst, a, k_unary_tanh_async) } + +// dst[i] = erf(a[i]) +func QErf(dst, a *data.Slice) { applyUnary(dst, a, k_unary_erf_async) } + +// dst[i] = erfc(a[i]) +func QErfc(dst, a *data.Slice) { applyUnary(dst, a, k_unary_erfc_async) } + +// dst[i] = tgamma(a[i]) +// returns 0 for a[i] outside of domain (0, inf) with poles at non-positive integers +func QGamma(dst, a *data.Slice) { applyUnary(dst, a, k_unary_gamma_async) } + +// dst[i] = 0 if a[i] < 0, 0.5 if a[i] == 0, 1 if a[i] > 0. +func QHeaviside(dst, a *data.Slice) { applyUnary(dst, a, k_unary_heaviside_async) } + +// dst[i] = sin(a[i])/a[i], 1 for a[i] == 0. +func QSinc(dst, a *data.Slice) { applyUnary(dst, a, k_unary_sinc_async) } + +// dst[i] = a[i]^b[i] +// pow(a,b) for negative a and b returns -pow(-a,b) for fractional exponents, and for 0^0 returns 1 +func QPow(dst, src1, src2 *data.Slice) { + N := dst.Len() + nComp := dst.NComp() + util.Assert(src1.Len() == N && src1.NComp() == nComp && src2.Len() == N && src2.NComp() == nComp) + + cfg := make1DConf(N) + for c := 0; c < nComp; c++ { + k_pw_pow_async(dst.DevPtr(c), src1.DevPtr(c), src2.DevPtr(c), N, cfg) + } +} + +// dst[i] = a[i] mod b[i] +func QMod(dst, src1, src2 *data.Slice) { + N := dst.Len() + nComp := dst.NComp() + util.Assert(src1.Len() == N && src1.NComp() == nComp && src2.Len() == N && src2.NComp() == nComp) + cfg := make1DConf(N) + for c := 0; c < nComp; c++ { + k_pw_mod_async(dst.DevPtr(c), src1.DevPtr(c), src2.DevPtr(c), N, cfg) + } +} + +// dst[i] = atan2(a[i], b[i]) +// returns 0 for atan2(0,0) +func QAtan2(dst, a, b *data.Slice) { + N := dst.Len() + nComp := dst.NComp() + util.Assert(a.Len() == N && a.NComp() == nComp && b.Len() == N && b.NComp() == nComp) + cfg := make1DConf(N) + for c := 0; c < nComp; c++ { + k_pw_atan2_async(dst.DevPtr(c), a.DevPtr(c), b.DevPtr(c), N, cfg) + } +} diff --git a/engine/unary_math.go b/engine/unary_math.go new file mode 100644 index 000000000..81c55fe7a --- /dev/null +++ b/engine/unary_math.go @@ -0,0 +1,137 @@ +package engine + +import ( + "fmt" + + "github.com/mumax/3/cuda" + "github.com/mumax/3/data" +) + +func init() { + DeclFunc("Qsin", QSin, "Element-wise sin of a Quantity (scalar or vector)") + DeclFunc("Qcos", QCos, "Element-wise cos of a Quantity (scalar or vector)") + DeclFunc("Qtan", QTan, "Element-wise tan of a Quantity (scalar or vector)") + DeclFunc("Qexp", QExp, "Element-wise exp of a Quantity (scalar or vector)") + DeclFunc("Qlog", QLog, "Element-wise log of a Quantity (scalar or vector)") + DeclFunc("Qabs", QAbs, "Element-wise absolute value of a Quantity (scalar or vector)") + DeclFunc("Qacos", QAcos, "Element-wise arccos of a Quantity (scalar or vector)") + DeclFunc("Qacosh", QAcosh, "Element-wise arccosh of a Quantity (scalar or vector)") + DeclFunc("Qasin", QAsin, "Element-wise arcsin of a Quantity (scalar or vector)") + DeclFunc("Qasinh", QAsinh, "Element-wise arcsinh of a Quantity (scalar or vector)") + DeclFunc("Qatan", QAtan, "Element-wise arctan of a Quantity (scalar or vector)") + DeclFunc("Qatan2", QAtan2, "Element-wise atan2(y, x) of two Quantities (scalar or vector)") + DeclFunc("Qatanh", QAtanh, "Element-wise arctanh of a Quantity (scalar or vector)") + DeclFunc("Qcosh", QCosh, "Element-wise cosh of a Quantity (scalar or vector)") + DeclFunc("Qsinh", QSinh, "Element-wise sinh of a Quantity (scalar or vector)") + DeclFunc("Qtanh", QTanh, "Element-wise tanh of a Quantity (scalar or vector)") + DeclFunc("Qerf", QErf, "Element-wise error function of a Quantity (scalar or vector)") + DeclFunc("Qerfc", QErfc, "Element-wise complementary error function of a Quantity (scalar or vector)") + DeclFunc("Qgamma", QGamma, "Element-wise gamma function of a Quantity (scalar or vector)") + DeclFunc("Qheaviside", QHeaviside, "Element-wise Heaviside step function of a Quantity (scalar or vector). Returns 0 for negative inputs, 1 for positive inputs, and 0.5 for zero.") + DeclFunc("Qsinc", QSinc, "Element-wise sinc function of a Quantity (scalar or vector)") + DeclFunc("Qpow", QPow, "Element-wise power: Qpow(base, exp). Supports fractional and negative powers, with the convention that for negative base and fractional exponent, the result is -pow(-base, exp), and for 0^0 returns 1.") + DeclFunc("Qmod", QMod, "Element-wise modulo: Qmod(a, b)") +} + +type unaryFuncQuantity struct { + a Quantity + f func(dst, a *data.Slice) + name string +} + +func (q *unaryFuncQuantity) NComp() int { return q.a.NComp() } + +func (q *unaryFuncQuantity) EvalTo(dst *data.Slice) { + a := ValueOf(q.a) + defer cuda.Recycle(a) + cuda.Zero(dst) + q.f(dst, a) +} + +func newUnaryQ(a Quantity, f func(dst, a *data.Slice), name string) Quantity { + return &unaryFuncQuantity{a, f, name} +} + +//pow, mod and atan2 are binary functions, so they need a separate struct and constructor + +type binaryFuncQuantity struct { + a, b Quantity + nComp int + f func(dst, a, b *data.Slice) + name string +} + +func (q *binaryFuncQuantity) NComp() int { return q.nComp } + +func (q *binaryFuncQuantity) EvalTo(dst *data.Slice) { + a := ValueOf(q.a) + defer cuda.Recycle(a) + b := ValueOf(q.b) + defer cuda.Recycle(b) + cuda.Zero(dst) + + switch { + //vector-vector, scalar-scalar, or same number of components + case a.NComp() == b.NComp(): + q.f(dst, a, b) + //vector-scalar + case a.NComp() == 1: + // broadcast a across each component + for c := 0; c < b.NComp(); c++ { + q.f(dst.Comp(c), a, b.Comp(c)) + } + //vector-scalar + case b.NComp() == 1: + // broadcast b across each component + for c := 0; c < a.NComp(); c++ { + q.f(dst.Comp(c), a.Comp(c), b) + } + } +} + +func newBinaryQ(a, b Quantity, f func(dst, a, b *data.Slice), name string) Quantity { + nComp := -1 + switch { + case a.NComp() == b.NComp(): + nComp = a.NComp() + case a.NComp() == 1: + nComp = b.NComp() + case b.NComp() == 1: + nComp = a.NComp() + default: + panic(fmt.Sprintf("Cannot apply %v to %v and %v components", name, a.NComp(), b.NComp())) + } + return &binaryFuncQuantity{a, b, nComp, f, name} +} + +// each unary function returns a Quantity representing the function applied pointwise +// for inputs outside of the function's domain, returns 0 +// limited domains for: tan poles at π/2 + nπ, acos: [-1, 1], acosh: [1, inf), asin: [-1, 1], atanh: (-1, 1), log: (0, inf), gamma: (0, inf) with poles at non-positive integers +func QSin(a Quantity) Quantity { return newUnaryQ(a, cuda.QSin, "Qsin") } +func QCos(a Quantity) Quantity { return newUnaryQ(a, cuda.QCos, "Qcos") } +func QTan(a Quantity) Quantity { return newUnaryQ(a, cuda.QTan, "Qtan") } +func QExp(a Quantity) Quantity { return newUnaryQ(a, cuda.QExp, "Qexp") } +func QLog(a Quantity) Quantity { return newUnaryQ(a, cuda.QLog, "Qlog") } +func QAbs(a Quantity) Quantity { return newUnaryQ(a, cuda.QAbs, "Qabs") } +func QAcos(a Quantity) Quantity { return newUnaryQ(a, cuda.QAcos, "Qacos") } +func QAcosh(a Quantity) Quantity { return newUnaryQ(a, cuda.QAcosh, "Qacosh") } +func QAsin(a Quantity) Quantity { return newUnaryQ(a, cuda.QAsin, "Qasin") } +func QAsinh(a Quantity) Quantity { return newUnaryQ(a, cuda.QAsinh, "Qasinh") } +func QAtan(a Quantity) Quantity { return newUnaryQ(a, cuda.QAtan, "Qatan") } +func QAtanh(a Quantity) Quantity { return newUnaryQ(a, cuda.QAtanh, "Qatanh") } +func QCosh(a Quantity) Quantity { return newUnaryQ(a, cuda.QCosh, "Qcosh") } +func QSinh(a Quantity) Quantity { return newUnaryQ(a, cuda.QSinh, "Qsinh") } +func QTanh(a Quantity) Quantity { return newUnaryQ(a, cuda.QTanh, "Qtanh") } +func QErf(a Quantity) Quantity { return newUnaryQ(a, cuda.QErf, "Qerf") } +func QErfc(a Quantity) Quantity { return newUnaryQ(a, cuda.QErfc, "Qerfc") } +func QGamma(a Quantity) Quantity { return newUnaryQ(a, cuda.QGamma, "Qgamma") } +func QHeaviside(a Quantity) Quantity { return newUnaryQ(a, cuda.QHeaviside, "Qheaviside") } +func QSinc(a Quantity) Quantity { return newUnaryQ(a, cuda.QSinc, "Qsinc") } + +// each binary function returns a Quantity representing the function applied pointwise +// for inputs outside of the function's domain, returns 0 +// pow(a,b), mod(a,b) and atan2(y,x) +// pow(a,b) for negative a and b returns -pow(-a,b) for fractional exponents, and for 0^0 returns 1 +func QPow(a, b Quantity) Quantity { return newBinaryQ(a, b, cuda.QPow, "Qpow") } +func QMod(a, b Quantity) Quantity { return newBinaryQ(a, b, cuda.QMod, "Qmod") } +func QAtan2(y, x Quantity) Quantity { return newBinaryQ(y, x, cuda.QAtan2, "Qatan2") } diff --git a/test/unarymath.mx3 b/test/unarymath.mx3 new file mode 100644 index 000000000..a2e9bda8c --- /dev/null +++ b/test/unarymath.mx3 @@ -0,0 +1,801 @@ +Nx := 32 +Ny :=32 +Nz := 2 +c := 5e-9 +setgridsize(Nx, Ny, Nz) +setcellsize(c, c, c) +Aex = 10e-12 +Msat = 600e3 +m = uniform(1, 0, 0) + +tol := 1e-6 + +// one mask per function, sized to its domain +vmask_sin := newVectorMask(Nx, Ny, Nz) +vmask_cos := newVectorMask(Nx, Ny, Nz) +vmask_tan := newVectorMask(Nx, Ny, Nz) +vmask_exp := newVectorMask(Nx, Ny, Nz) +vmask_log := newVectorMask(Nx, Ny, Nz) +vmask_abs := newVectorMask(Nx, Ny, Nz) +vmask_acos := newVectorMask(Nx, Ny, Nz) +vmask_acosh := newVectorMask(Nx, Ny, Nz) +vmask_asin := newVectorMask(Nx, Ny, Nz) +vmask_asinh := newVectorMask(Nx, Ny, Nz) +vmask_atan := newVectorMask(Nx, Ny, Nz) +vmask_atanh := newVectorMask(Nx, Ny, Nz) +vmask_cosh := newVectorMask(Nx, Ny, Nz) +vmask_sinh := newVectorMask(Nx, Ny, Nz) +vmask_tanh := newVectorMask(Nx, Ny, Nz) +vmask_erf := newVectorMask(Nx, Ny, Nz) +vmask_erfc := newVectorMask(Nx, Ny, Nz) +vmask_gamma := newVectorMask(Nx, Ny, Nz) +vmask_heaviside := newVectorMask(Nx, Ny, Nz) +vmask_sinc := newVectorMask(Nx, Ny, Nz) +vmask_pow_base := newVectorMask(Nx, Ny, Nz) +vmask_pow_exp := newVectorMask(Nx, Ny, Nz) +vmask_mod_a := newVectorMask(Nx, Ny, Nz) +vmask_mod_b := newVectorMask(Nx, Ny, Nz) +vmask_atan2_y := newVectorMask(Nx, Ny, Nz) +vmask_atan2_x := newVectorMask(Nx, Ny, Nz) + +for i := 0; i < Nx; i++ { + for j := 0; j < Ny; j++ { + for k := 0; k < Nz; k++ { + // sin/cos: (-1.57, 1.57) + vmask_sin.setVector(i,j,k, vector((rand()-0.5)*3.14, (rand()-0.5)*3.14, (rand()-0.5)*3.14)) + vmask_cos.setVector(i,j,k, vector((rand()-0.5)*3.14, (rand()-0.5)*3.14, (rand()-0.5)*3.14)) + // tan: (-1.5, 1.5) — away from pi/2 singularities + vmask_tan.setVector(i,j,k, vector((rand()-0.5)*3.0, (rand()-0.5)*3.0, (rand()-0.5)*3.0)) + // exp: (-1, 1) + vmask_exp.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // log: (0.1, 2.1) — strictly positive + vmask_log.setVector(i,j,k, vector(rand()*2.0+0.1, rand()*2.0+0.1, rand()*2.0+0.1)) + // abs: (-2, 2) + vmask_abs.setVector(i,j,k, vector((rand()-0.5)*4.0, (rand()-0.5)*4.0, (rand()-0.5)*4.0)) + // acos/asin: (-0.99, 0.99) + vmask_acos.setVector(i,j,k, vector((rand()-0.5)*1.98, (rand()-0.5)*1.98, (rand()-0.5)*1.98)) + vmask_asin.setVector(i,j,k, vector((rand()-0.5)*1.98, (rand()-0.5)*1.98, (rand()-0.5)*1.98)) + // acosh: (1.01, 3.01) + vmask_acosh.setVector(i,j,k, vector(rand()*2.0+1.01, rand()*2.0+1.01, rand()*2.0+1.01)) + // asinh: (-1, 1) + vmask_asinh.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // atan: (-1, 1) + vmask_atan.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // atanh: (-0.9, 0.9) + vmask_atanh.setVector(i,j,k, vector((rand()-0.5)*1.8, (rand()-0.5)*1.8, (rand()-0.5)*1.8)) + // cosh: (-1, 1) + vmask_cosh.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // sinh: (-1, 1) + vmask_sinh.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // tanh: (-2, 2) + vmask_tanh.setVector(i,j,k, vector((rand()-0.5)*4.0, (rand()-0.5)*4.0, (rand()-0.5)*4.0)) + // erf: (-1, 1) + vmask_erf.setVector(i,j,k, vector((rand()-0.5)*2.0, (rand()-0.5)*2.0, (rand()-0.5)*2.0)) + // erfc: (0, 1.5) + vmask_erfc.setVector(i,j,k, vector(rand()*1.5, rand()*1.5, rand()*1.5)) + // gamma: (0.5, 2.5) — avoids poles at 0 and negative integers + vmask_gamma.setVector(i,j,k, vector(rand()*2.0+0.5, rand()*2.0+0.5, rand()*2.0+0.5)) + // heaviside: (-2, 2) + vmask_heaviside.setVector(i,j,k, vector((rand()-0.5)*4.0, (rand()-0.5)*4.0, (rand()-0.5)*4.0)) + // sinc: (-pi, pi) + vmask_sinc.setVector(i,j,k, vector((rand()-0.5)*6.28, (rand()-0.5)*6.28, (rand()-0.5)*6.28)) + // pow: base (0.1, 2.1), exp (0.5, 2.5) + vmask_pow_base.setVector(i,j,k, vector(rand()*2.0+0.1, rand()*2.0+0.1, rand()*2.0+0.1)) + vmask_pow_exp.setVector(i,j,k, vector(rand()*2.0+0.5, rand()*2.0+0.5, rand()*2.0+0.5)) + // mod: dividend (-5, 5), divisor (0.5, 2.5) + vmask_mod_a.setVector(i,j,k, vector((rand()-0.5)*10.0, (rand()-0.5)*10.0, (rand()-0.5)*10.0)) + vmask_mod_b.setVector(i,j,k, vector(rand()*2.0+0.5, rand()*2.0+0.5, rand()*2.0+0.5)) + // atan2: (-2, 2) for both + vmask_atan2_y.setVector(i,j,k, vector((rand()-0.5)*4.0, (rand()-0.5)*4.0, (rand()-0.5)*4.0)) + vmask_atan2_x.setVector(i,j,k, vector((rand()-0.5)*4.0, (rand()-0.5)*4.0, (rand()-0.5)*4.0)) + } + } +} + +// precompute reference masks +ref_sin := newVectorMask(Nx, Ny, Nz) +ref_cos := newVectorMask(Nx, Ny, Nz) +ref_tan := newVectorMask(Nx, Ny, Nz) +ref_exp := newVectorMask(Nx, Ny, Nz) +ref_log := newVectorMask(Nx, Ny, Nz) +ref_abs := newVectorMask(Nx, Ny, Nz) +ref_acos := newVectorMask(Nx, Ny, Nz) +ref_acosh := newVectorMask(Nx, Ny, Nz) +ref_asin := newVectorMask(Nx, Ny, Nz) +ref_asinh := newVectorMask(Nx, Ny, Nz) +ref_atan := newVectorMask(Nx, Ny, Nz) +ref_atanh := newVectorMask(Nx, Ny, Nz) +ref_cosh := newVectorMask(Nx, Ny, Nz) +ref_sinh := newVectorMask(Nx, Ny, Nz) +ref_tanh := newVectorMask(Nx, Ny, Nz) +ref_erf := newVectorMask(Nx, Ny, Nz) +ref_erfc := newVectorMask(Nx, Ny, Nz) +ref_gamma := newVectorMask(Nx, Ny, Nz) +ref_heaviside := newVectorMask(Nx, Ny, Nz) +ref_sinc := newVectorMask(Nx, Ny, Nz) +ref_pow := newVectorMask(Nx, Ny, Nz) +ref_mod := newVectorMask(Nx, Ny, Nz) +ref_atan2 := newVectorMask(Nx, Ny, Nz) + +for i := 0; i < Nx; i++ { + for j := 0; j < Ny; j++ { + for k := 0; k < Nz; k++ { + ref_sin.setVector(i,j,k, vector(sin(vmask_sin.Get(0,i,j,k)), sin(vmask_sin.Get(1,i,j,k)), sin(vmask_sin.Get(2,i,j,k)))) + ref_cos.setVector(i,j,k, vector(cos(vmask_cos.Get(0,i,j,k)), cos(vmask_cos.Get(1,i,j,k)), cos(vmask_cos.Get(2,i,j,k)))) + ref_tan.setVector(i,j,k, vector(tan(vmask_tan.Get(0,i,j,k)), tan(vmask_tan.Get(1,i,j,k)), tan(vmask_tan.Get(2,i,j,k)))) + ref_exp.setVector(i,j,k, vector(exp(vmask_exp.Get(0,i,j,k)), exp(vmask_exp.Get(1,i,j,k)), exp(vmask_exp.Get(2,i,j,k)))) + ref_log.setVector(i,j,k, vector(log(vmask_log.Get(0,i,j,k)), log(vmask_log.Get(1,i,j,k)), log(vmask_log.Get(2,i,j,k)))) + ref_abs.setVector(i,j,k, vector(abs(vmask_abs.Get(0,i,j,k)), abs(vmask_abs.Get(1,i,j,k)), abs(vmask_abs.Get(2,i,j,k)))) + ref_acos.setVector(i,j,k, vector(acos(vmask_acos.Get(0,i,j,k)), acos(vmask_acos.Get(1,i,j,k)), acos(vmask_acos.Get(2,i,j,k)))) + ref_acosh.setVector(i,j,k, vector(acosh(vmask_acosh.Get(0,i,j,k)), acosh(vmask_acosh.Get(1,i,j,k)), acosh(vmask_acosh.Get(2,i,j,k)))) + ref_asin.setVector(i,j,k, vector(asin(vmask_asin.Get(0,i,j,k)), asin(vmask_asin.Get(1,i,j,k)), asin(vmask_asin.Get(2,i,j,k)))) + ref_asinh.setVector(i,j,k, vector(asinh(vmask_asinh.Get(0,i,j,k)), asinh(vmask_asinh.Get(1,i,j,k)), asinh(vmask_asinh.Get(2,i,j,k)))) + ref_atan.setVector(i,j,k, vector(atan(vmask_atan.Get(0,i,j,k)), atan(vmask_atan.Get(1,i,j,k)), atan(vmask_atan.Get(2,i,j,k)))) + ref_atanh.setVector(i,j,k, vector(atanh(vmask_atanh.Get(0,i,j,k)), atanh(vmask_atanh.Get(1,i,j,k)), atanh(vmask_atanh.Get(2,i,j,k)))) + ref_cosh.setVector(i,j,k, vector(cosh(vmask_cosh.Get(0,i,j,k)), cosh(vmask_cosh.Get(1,i,j,k)), cosh(vmask_cosh.Get(2,i,j,k)))) + ref_sinh.setVector(i,j,k, vector(sinh(vmask_sinh.Get(0,i,j,k)), sinh(vmask_sinh.Get(1,i,j,k)), sinh(vmask_sinh.Get(2,i,j,k)))) + ref_tanh.setVector(i,j,k, vector(tanh(vmask_tanh.Get(0,i,j,k)), tanh(vmask_tanh.Get(1,i,j,k)), tanh(vmask_tanh.Get(2,i,j,k)))) + ref_erf.setVector(i,j,k, vector(erf(vmask_erf.Get(0,i,j,k)), erf(vmask_erf.Get(1,i,j,k)), erf(vmask_erf.Get(2,i,j,k)))) + ref_erfc.setVector(i,j,k, vector(erfc(vmask_erfc.Get(0,i,j,k)), erfc(vmask_erfc.Get(1,i,j,k)), erfc(vmask_erfc.Get(2,i,j,k)))) + ref_gamma.setVector(i,j,k, vector(gamma(vmask_gamma.Get(0,i,j,k)), gamma(vmask_gamma.Get(1,i,j,k)), gamma(vmask_gamma.Get(2,i,j,k)))) + ref_heaviside.setVector(i,j,k,vector(heaviside(vmask_heaviside.Get(0,i,j,k)), heaviside(vmask_heaviside.Get(1,i,j,k)), heaviside(vmask_heaviside.Get(2,i,j,k)))) + ref_sinc.setVector(i,j,k, vector(sinc(vmask_sinc.Get(0,i,j,k)), sinc(vmask_sinc.Get(1,i,j,k)), sinc(vmask_sinc.Get(2,i,j,k)))) + ref_pow.setVector(i,j,k, vector(pow(vmask_pow_base.Get(0,i,j,k), vmask_pow_exp.Get(0,i,j,k)), pow(vmask_pow_base.Get(1,i,j,k), vmask_pow_exp.Get(1,i,j,k)), pow(vmask_pow_base.Get(2,i,j,k), vmask_pow_exp.Get(2,i,j,k)))) + ref_mod.setVector(i,j,k, vector(mod(vmask_mod_a.Get(0,i,j,k), vmask_mod_b.Get(0,i,j,k)), mod(vmask_mod_a.Get(1,i,j,k), vmask_mod_b.Get(1,i,j,k)), mod(vmask_mod_a.Get(2,i,j,k), vmask_mod_b.Get(2,i,j,k)))) + ref_atan2.setVector(i,j,k, vector(atan2(vmask_atan2_y.Get(0,i,j,k), vmask_atan2_x.Get(0,i,j,k)), atan2(vmask_atan2_y.Get(1,i,j,k), vmask_atan2_x.Get(1,i,j,k)), atan2(vmask_atan2_y.Get(2,i,j,k), vmask_atan2_x.Get(2,i,j,k)))) + } + } +} + +// test points +p1x := 0; p1y := 15; p1z := 1 +p2x := 10; p2y := 2; p2z := 0 +p3x := 25; p3y := 12; p3z := 1 +p4x := 4; p4y := 4; p4z := 0 +p5x := 31; p5y := 20; p5z := 0 + +OutputFormat = OVF2_TEXT + +// --- sin --- +B_ext.add(ref_sin, 1) +AddFieldTerm(Qsin(CustomQuantity(vmask_sin))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref := LoadFile("testmath2.out/B_ext000000.ovf") +test := LoadFile("testmath2.out/B_custom000000.ovf") +Expect("sin (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("sin (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("sin (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("sin (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("sin (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("sin (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("sin (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("sin (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("sin (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("sin (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("sin (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("sin (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("sin (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("sin (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("sin (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- cos --- +B_ext.add(ref_cos, 1) +AddFieldTerm(Qcos(CustomQuantity(vmask_cos))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000001.ovf") +test = LoadFile("testmath2.out/B_custom000001.ovf") +Expect("cos (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("cos (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("cos (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("cos (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("cos (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("cos (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("cos (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("cos (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("cos (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("cos (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("cos (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("cos (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("cos (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("cos (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("cos (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- tan --- +B_ext.add(ref_tan, 1) +AddFieldTerm(Qtan(CustomQuantity(vmask_tan))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000002.ovf") +test = LoadFile("testmath2.out/B_custom000002.ovf") +Expect("tan (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("tan (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("tan (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("tan (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("tan (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("tan (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("tan (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("tan (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("tan (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("tan (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("tan (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("tan (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("tan (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("tan (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("tan (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- exp --- +B_ext.add(ref_exp, 1) +AddFieldTerm(Qexp(CustomQuantity(vmask_exp))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000003.ovf") +test = LoadFile("testmath2.out/B_custom000003.ovf") +Expect("exp (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("exp (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("exp (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("exp (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("exp (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("exp (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("exp (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("exp (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("exp (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("exp (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("exp (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("exp (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("exp (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("exp (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("exp (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- log --- +B_ext.add(ref_log, 1) +AddFieldTerm(Qlog(CustomQuantity(vmask_log))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000004.ovf") +test = LoadFile("testmath2.out/B_custom000004.ovf") +Expect("log (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("log (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("log (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("log (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("log (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("log (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("log (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("log (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("log (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("log (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("log (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("log (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("log (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("log (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("log (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- abs --- +B_ext.add(ref_abs, 1) +AddFieldTerm(Qabs(CustomQuantity(vmask_abs))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000005.ovf") +test = LoadFile("testmath2.out/B_custom000005.ovf") +Expect("abs (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("abs (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("abs (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("abs (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("abs (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("abs (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("abs (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("abs (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("abs (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("abs (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("abs (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("abs (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("abs (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("abs (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("abs (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- acos --- +B_ext.add(ref_acos, 1) +AddFieldTerm(Qacos(CustomQuantity(vmask_acos))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000006.ovf") +test = LoadFile("testmath2.out/B_custom000006.ovf") +Expect("acos (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("acos (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("acos (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("acos (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("acos (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("acos (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("acos (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("acos (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("acos (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("acos (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("acos (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("acos (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("acos (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("acos (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("acos (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- acosh --- +B_ext.add(ref_acosh, 1) +AddFieldTerm(Qacosh(CustomQuantity(vmask_acosh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000007.ovf") +test = LoadFile("testmath2.out/B_custom000007.ovf") +Expect("acosh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("acosh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("acosh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("acosh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("acosh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("acosh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("acosh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("acosh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("acosh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("acosh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("acosh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("acosh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("acosh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("acosh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("acosh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- asin --- +B_ext.add(ref_asin, 1) +AddFieldTerm(Qasin(CustomQuantity(vmask_asin))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000008.ovf") +test = LoadFile("testmath2.out/B_custom000008.ovf") +Expect("asin (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("asin (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("asin (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("asin (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("asin (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("asin (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("asin (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("asin (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("asin (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("asin (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("asin (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("asin (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("asin (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("asin (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("asin (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- asinh --- +B_ext.add(ref_asinh, 1) +AddFieldTerm(Qasinh(CustomQuantity(vmask_asinh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000009.ovf") +test = LoadFile("testmath2.out/B_custom000009.ovf") +Expect("asinh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("asinh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("asinh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("asinh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("asinh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("asinh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("asinh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("asinh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("asinh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("asinh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("asinh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("asinh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("asinh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("asinh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("asinh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- atan --- +B_ext.add(ref_atan, 1) +AddFieldTerm(Qatan(CustomQuantity(vmask_atan))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000010.ovf") +test = LoadFile("testmath2.out/B_custom000010.ovf") +Expect("atan (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("atan (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("atan (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("atan (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("atan (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("atan (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("atan (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("atan (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("atan (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("atan (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("atan (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("atan (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("atan (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("atan (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("atan (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- atanh --- +B_ext.add(ref_atanh, 1) +AddFieldTerm(Qatanh(CustomQuantity(vmask_atanh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000011.ovf") +test = LoadFile("testmath2.out/B_custom000011.ovf") +Expect("atanh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("atanh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("atanh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("atanh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("atanh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("atanh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("atanh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("atanh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("atanh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("atanh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("atanh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("atanh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("atanh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("atanh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("atanh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- cosh --- +B_ext.add(ref_cosh, 1) +AddFieldTerm(Qcosh(CustomQuantity(vmask_cosh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000012.ovf") +test = LoadFile("testmath2.out/B_custom000012.ovf") +Expect("cosh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("cosh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("cosh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("cosh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("cosh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("cosh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("cosh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("cosh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("cosh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("cosh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("cosh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("cosh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("cosh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("cosh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("cosh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- sinh --- +B_ext.add(ref_sinh, 1) +AddFieldTerm(Qsinh(CustomQuantity(vmask_sinh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000013.ovf") +test = LoadFile("testmath2.out/B_custom000013.ovf") +Expect("sinh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("sinh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("sinh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("sinh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("sinh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("sinh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("sinh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("sinh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("sinh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("sinh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("sinh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("sinh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("sinh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("sinh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("sinh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- tanh --- +B_ext.add(ref_tanh, 1) +AddFieldTerm(Qtanh(CustomQuantity(vmask_tanh))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000014.ovf") +test = LoadFile("testmath2.out/B_custom000014.ovf") +Expect("tanh (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("tanh (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("tanh (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("tanh (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("tanh (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("tanh (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("tanh (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("tanh (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("tanh (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("tanh (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("tanh (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("tanh (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("tanh (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("tanh (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("tanh (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- erf --- +B_ext.add(ref_erf, 1) +AddFieldTerm(Qerf(CustomQuantity(vmask_erf))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000015.ovf") +test = LoadFile("testmath2.out/B_custom000015.ovf") +Expect("erf (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("erf (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("erf (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("erf (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("erf (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("erf (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("erf (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("erf (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("erf (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("erf (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("erf (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("erf (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("erf (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("erf (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("erf (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- erfc --- +B_ext.add(ref_erfc, 1) +AddFieldTerm(Qerfc(CustomQuantity(vmask_erfc))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000016.ovf") +test = LoadFile("testmath2.out/B_custom000016.ovf") +Expect("erfc (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("erfc (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("erfc (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("erfc (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("erfc (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("erfc (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("erfc (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("erfc (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("erfc (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("erfc (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("erfc (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("erfc (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("erfc (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("erfc (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("erfc (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- gamma --- +B_ext.add(ref_gamma, 1) +AddFieldTerm(Qgamma(CustomQuantity(vmask_gamma))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000017.ovf") +test = LoadFile("testmath2.out/B_custom000017.ovf") +Expect("gamma (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("gamma (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("gamma (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("gamma (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("gamma (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("gamma (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("gamma (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("gamma (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("gamma (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("gamma (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("gamma (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("gamma (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("gamma (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("gamma (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("gamma (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- heaviside --- +B_ext.add(ref_heaviside, 1) +AddFieldTerm(Qheaviside(CustomQuantity(vmask_heaviside))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000018.ovf") +test = LoadFile("testmath2.out/B_custom000018.ovf") +Expect("heaviside (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("heaviside (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("heaviside (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("heaviside (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("heaviside (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("heaviside (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("heaviside (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("heaviside (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("heaviside (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("heaviside (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("heaviside (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("heaviside (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("heaviside (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("heaviside (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("heaviside (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- sinc --- +B_ext.add(ref_sinc, 1) +AddFieldTerm(Qsinc(CustomQuantity(vmask_sinc))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000019.ovf") +test = LoadFile("testmath2.out/B_custom000019.ovf") +Expect("sinc (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("sinc (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("sinc (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("sinc (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("sinc (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("sinc (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("sinc (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("sinc (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("sinc (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("sinc (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("sinc (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("sinc (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("sinc (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("sinc (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("sinc (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- pow --- +B_ext.add(ref_pow, 1) +AddFieldTerm(Qpow(CustomQuantity(vmask_pow_base), CustomQuantity(vmask_pow_exp))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000020.ovf") +test = LoadFile("testmath2.out/B_custom000020.ovf") +Expect("pow (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("pow (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("pow (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("pow (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("pow (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("pow (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("pow (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("pow (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("pow (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("pow (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("pow (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("pow (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("pow (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("pow (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("pow (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- mod --- +B_ext.add(ref_mod, 1) +AddFieldTerm(Qmod(CustomQuantity(vmask_mod_a), CustomQuantity(vmask_mod_b))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000021.ovf") +test = LoadFile("testmath2.out/B_custom000021.ovf") +Expect("mod (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("mod (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("mod (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("mod (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("mod (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("mod (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("mod (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("mod (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("mod (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("mod (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("mod (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("mod (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("mod (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("mod (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("mod (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() + +// --- atan2 --- +B_ext.add(ref_atan2, 1) +AddFieldTerm(Qatan2(CustomQuantity(vmask_atan2_y), CustomQuantity(vmask_atan2_x))) +m = uniform(1, 0, 0) + +save(B_ext) +save(B_custom) +ref = LoadFile("testmath2.out/B_ext000022.ovf") +test = LoadFile("testmath2.out/B_custom000022.ovf") +Expect("atan2 (0,15,1) x", test.Get(0,p1x,p1y,p1z), ref.Get(0,p1x,p1y,p1z), tol) +Expect("atan2 (0,15,1) y", test.Get(1,p1x,p1y,p1z), ref.Get(1,p1x,p1y,p1z), tol) +Expect("atan2 (0,15,1) z", test.Get(2,p1x,p1y,p1z), ref.Get(2,p1x,p1y,p1z), tol) +Expect("atan2 (64,32,0) x", test.Get(0,p2x,p2y,p2z), ref.Get(0,p2x,p2y,p2z), tol) +Expect("atan2 (64,32,0) y", test.Get(1,p2x,p2y,p2z), ref.Get(1,p2x,p2y,p2z), tol) +Expect("atan2 (64,32,0) z", test.Get(2,p2x,p2y,p2z), ref.Get(2,p2x,p2y,p2z), tol) +Expect("atan2 (77,10,1) x", test.Get(0,p3x,p3y,p3z), ref.Get(0,p3x,p3y,p3z), tol) +Expect("atan2 (77,10,1) y", test.Get(1,p3x,p3y,p3z), ref.Get(1,p3x,p3y,p3z), tol) +Expect("atan2 (77,10,1) z", test.Get(2,p3x,p3y,p3z), ref.Get(2,p3x,p3y,p3z), tol) +Expect("atan2 (33,50,1) x", test.Get(0,p4x,p4y,p4z), ref.Get(0,p4x,p4y,p4z), tol) +Expect("atan2 (33,50,1) y", test.Get(1,p4x,p4y,p4z), ref.Get(1,p4x,p4y,p4z), tol) +Expect("atan2 (33,50,1) z", test.Get(2,p4x,p4y,p4z), ref.Get(2,p4x,p4y,p4z), tol) +Expect("atan2 (100,20,0) x", test.Get(0,p5x,p5y,p5z), ref.Get(0,p5x,p5y,p5z), tol) +Expect("atan2 (100,20,0) y", test.Get(1,p5x,p5y,p5z), ref.Get(1,p5x,p5y,p5z), tol) +Expect("atan2 (100,20,0) z", test.Get(2,p5x,p5y,p5z), ref.Get(2,p5x,p5y,p5z), tol) +RemoveCustomFields() +RemoveCustomEnergies() +B_ext.RemoveExtraTerms() \ No newline at end of file