Source code for gusto.rexi.rexi_coefficients

import numpy
from gusto import Configuration


[docs] class RexiParameters(Configuration): """ mu and a coefficients from "A high-order time-parallel scheme for solving wave propagation problems via the direct construction of an approximate time-evolution operator", Haut et.al. """ h = 0.2 M = 64 reduce_to_half = False mu = -4.315321510875024 + 1j*0 # The variables below this line should not be changed. L = 11 a = [ -1.0845749544592896e-7 + 1j*2.77075431662228e-8, 1.858753344202957e-8 + 1j*-9.105375434750162e-7, 3.6743713227243024e-6 + 1j*7.073284346322969e-7, -2.7990058083347696e-6 + 1j*0.0000112564827639346, 0.000014918577548849352 + 1j*-0.0000316278486761932, -0.0010751767283285608 + 1j*-0.00047282220513073084, 0.003816465653840016 + 1j*0.017839810396560574, 0.12124105653274578 + 1j*-0.12327042473830248, -0.9774980792734348 + 1j*-0.1877130220537587, 1.3432866123333178 + 1j*3.2034715228495942, 4.072408546157305 + 1j*-6.123755543580666, -9.442699917778205 + 1j*0., 4.072408620272648 + 1j*6.123755841848161, 1.3432860877712938 + 1j*-3.2034712658530275, -0.9774985292598916 + 1j*0.18771238018072134, 0.1212417070363373 + 1j*0.12326987628935386, 0.0038169724770333343 + 1j*-0.017839242222443888, -0.0010756025812659208 + 1j*0.0004731874917343858, 0.000014713754789095218 + 1j*0.000031358475831136815, -2.659323898804944e-6 + 1j*-0.000011341571201752273, 3.6970377676364553e-6 + 1j*-6.517457477594937e-7, 3.883933649142257e-9 + 1j*9.128496023863376e-7, -1.0816457995911385e-7 + 1j*-2.954309729192276e-8 ]
[docs] def b_coefficients(h, M): """ Compute the b coefficients where b_m = h^2 exp(imh) """ m = numpy.arange(-M, M+1) return numpy.exp(h*h)*numpy.exp(-1j*m*h)
[docs] def RexiCoefficients(rexi_parameters): """ Compute the A_n and B_n coefficients in the REXI sum exp(ix) \approx sum_{n=0}^{P} B_n (ix + A_n) where P = 4N+1 if rexi_parameters.reduce_to_half is False else 2N+1 if True. Returns 3 numpy arrays: alpha contains the A_n coefficients beta contains the B_n coefficients beta2 contains zeros if rexi_parameters.reduce_to_half is False else it contains the coefficients that multiply the conjugate terms :arg rexi_parameters: class containing the parameters (h, M and reduce_to_half) necessary to compute the coefficients """ h = float(rexi_parameters.h) M = int(float(rexi_parameters.M)) # get L, mu and the a coefficients L = rexi_parameters.L mu = rexi_parameters.mu a = rexi_parameters.a # calculate the b coefficients b = b_coefficients(h, M) # allocate arrays for alpha, beta_re and beta_im N = M + L alpha = numpy.zeros((2*N+1,), dtype=numpy.complex128) beta_re = numpy.zeros((2*N+1,), dtype=numpy.complex128) beta_im = numpy.zeros((2*N+1,), dtype=numpy.complex128) # compute alpha, beta_re and beta_im for l in range(-L, L+1): for m in range(-M, M+1): n = l+m alpha[n+N] = h*(mu + 1j*n) beta_re[n+N] += b[m+M].real*h*a[l+L] beta_im[n+N] += b[m+M].imag*h*a[l+L] # calculate conj(beta_re) and conj(beta_im), used to define A_n and B_n beta_conj_re = numpy.conjugate(beta_re) beta_conj_im = numpy.conjugate(beta_im) if rexi_parameters.reduce_to_half: # If reducing the number of solvers to (nearly) half, as # described in notes.pdf, we only need alpha_n for n \in [N, # 2N]. We need to retain all the betas (due to the lack of # symmetry coming from the numerical calculation of the a # coefficients) but in order not to special case the Nth # solver (where the solution is real hence equal to its # conjugate) we must divide the Nth value of the betas by 2. alpha = alpha[N:] beta_re[N] /= 2. beta_im[N] /= 2. beta_conj_re[N] /= 2. beta_conj_im[N] /= 2. beta = numpy.concatenate( (beta_re[N:] + 1j*beta_im[N:], -beta_conj_re[N::-1] - 1j*beta_conj_im[N::-1]) )/2 # beta2 is the coefficient that multiplies the conjugate of # the solution beta2 = numpy.concatenate( (beta_re[N::-1] + 1j*beta_im[N::-1], -beta_conj_re[N:] - 1j*beta_conj_im[N:]) )/2 else: beta = numpy.concatenate( (beta_re + 1j*beta_im, -beta_conj_re[::-1] - 1j*beta_conj_im[::-1]) )/2 # when not reducing the number of solvers we return zero here # in order not to special case this option (i.e. we still # calculate the conjugate of the solution but then multiply it # by zero - you're doing (nearly) twice the amount of work necessary # anyway!) beta2 = numpy.zeros(len(beta)) alpha = numpy.concatenate((alpha, -alpha)) return alpha, beta, beta2