diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py new file mode 100644 index 0000000..b47e5ce --- /dev/null +++ b/examples/shestakov_time_dependent.py @@ -0,0 +1,418 @@ +"""Example for how to use tango to solve a turbulence and transport problem. + +Using a solver class, saving to a file, and using the analysis package to load the data and save a plot + +Here, the "turbulent flux" is specified analytically, using the example in the Shestakov et al. (2003) paper. +This example is a nonlinear diffusion equation with specified diffusion coefficient and source. There is a +closed form answer for the steady state solution which can be compared with the numerically found solution. +""" + +from __future__ import division, absolute_import +import numpy as np +import matplotlib.pyplot as plt + +import tango.tango_logging as tlog +from tango.extras.shestakov_nonlinear_diffusion import ShestakovTestProblem +from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation, FluxRelaxationOscillation, FluxAverage +from tango.extras.noisyflux import NoisyFlux, NoisyFluxSpaceTime +import tango + +class ShestakovProblem(ShestakovTestProblem): + def __init__(self, + p = 3, + q = -2, + L = 1, # size of domain + N = 500, # number of spatial grid points + rightBC = 1e-2, # right boundary condition + initial = lambda x: 1 - 0.5*x): + self.L = L # size of domain + self.N = N # number of spatial grid points + self.dx = L / (N - 1) # spatial grid size + self.x = np.arange(N) * self.dx # location corresponding to grid points j=0, ..., N-1 + self._rightBC = rightBC # right boundary condition + self._initial = initial + + super(ShestakovProblem, self).__init__(self.dx, p, q) + + def InitialCondition(self): + return self._initial(self.x) + + def GetSource(self, x = None): + if x is None: + x = self.x + return super(ShestakovProblem, self).GetSource(x) + + def steady_state_solution(self): + return super(ShestakovProblem, self).steady_state_solution( + self.x, + self.RightBC()) + + def RightBC(self): + return self._rightBC + + def ComputeAllH(self, t, x, profiles, HCoeffsTurb): + H1 = np.ones_like(x) + H7 = self.H7contrib_Source(x) + return HCoeffsTurb + tango.multifield.HCoefficients(H1=H1, H7=H7) + +# ============================================================================== +# Solve function +# ============================================================================== + +def solve_system(noise_timescale = 1.0, # AR time of random noise + flux_timescale = 1.0, # Flux relaxation time + oscillation_amplitude = 0.0, # Flux oscillation relative amplitude + noise_amplitude = 0.0, # Amplitude of AR(1) noise + noise_scalelength = 10, # Spatial scale length (cells) + EWMAParamTurbFlux = 0.01, + EWMAParamProfile = 1.0, + p = 3, + q = -2, + thetaParams = {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10}, + tol = 1e-2, + maxIterations = 1000, + plot_convergence = True, + check_tol = 0.1 + ): + tlog.info("Initializing...") + + nsteps = 1 + + # test_problem = ShestakovProblem(p=p, q=q) + + test_problem = FluxRelaxationOscillation( + flux_timescale * nsteps, + oscillation_amplitude, + ShestakovProblem(p=p, q=q)) + + # fluxModel = FluxAverage( + # nsteps, # Number of steps + # NoisyFluxSpaceTime( + # noise_amplitude, + # noise_scalelength, + # 1.0, # dx + # noise_timescale * nsteps, # noise timescale + # 1.0, # dt + # FluxRelaxationOscillation( + # flux_timescale * nsteps, + # oscillation_amplitude, + # test_problem))) + + turbHandler = tango.lodestro_method.TurbulenceHandler( + test_problem.dx, test_problem.x, test_problem) + + lodestroMethod = tango.lodestro_method.lm( + EWMAParamTurbFlux, + EWMAParamProfile, + thetaParams, + ) + field0 = tango.multifield.Field( + label = test_problem.label(), + rightBC = test_problem.RightBC(), + profile_mminus1 = test_problem.InitialCondition(), + compute_all_H = test_problem.ComputeAllH, + lodestroMethod = lodestroMethod, + ) + fields = [field0] + tango.multifield.check_fields_initialize(fields) + + compute_all_H_all_fields = tango.multifield.ComputeAllHAllFields( + fields, turbHandler + ) + + tArray = np.array([0, 1e4]) # specify the timesteps to be used. + + solver = tango.solver.Solver( + test_problem.L, + test_problem.x, + tArray, + maxIterations, + tol, + compute_all_H_all_fields, + fields, + saveFluxesInMemory = plot_convergence + ) + + tlog.info("Initialization complete.") + tlog.info("Entering main time loop...") + + err_history = [] + profile_history = [] + flux_history = [] + num_iterations = 0 + restarts = [] + while num_iterations < maxIterations: + try: + while solver.ok: + # Implicit time advance: iterate to solve the nonlinear equation! + solver.take_timestep() + except Exception as e: + # Failed + print("Solver failed with error: {}".format(e)) + solver.reachedEnd = False + solutionResidual = 1e10 + solutionRmsError = 1e10 + num_iterations = maxIterations + break + + label = test_problem.label() + + n = solver.profiles[label] # finished solution + + source = test_problem.GetSource() + nss = test_problem.steady_state_solution() + + solutionResidual = (n - nss) / np.max(np.abs(nss)) + solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) + + num_iterations += len(solver.errHistoryFinal) + err_history.append(solver.errHistoryFinal.copy()) + profile_history.append(solver.profilesAllIterations[label].copy()) + flux_history.append(solver.fluxesAllIterations[label].copy()) + + if solver.reachedEnd == True: + print("The solution has been reached successfully.") + print( + "Error compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + + if solutionRmsError > check_tol: + print("False convergence!") + solver.m = 0 # Reset timestep index + solver.t = 0.0 + + restarts.append(num_iterations) + else: + break + else: + print("The solver failed for some reason.") + print( + "Error at end compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + if solutionRmsError < check_tol: + # Succeeded anyway + restarts.append(num_iterations) + solver.reachedEnd = True + else: + num_iterations = maxIterations + break # Give up + + if err_history != []: + err_history = np.concatenate(err_history) + profiles = np.concatenate(profile_history) + fluxes = np.concatenate(flux_history) + + if plot_convergence: + for j in range(num_iterations): + fig, [ax1, ax2, ax3] = plt.subplots(1,3) + + # Profile and flux history in grey + for i in range(num_iterations): + ax1.plot(profiles[i,:], color='0.6') + ax2.plot(fluxes[i,:], color='0.6') + + # Plot the converged profiles and fluxes in black + ax1.plot(profiles[-1,:], color='k') + ax2.plot(fluxes[-1,:], color='k') + + # Current step in blue + ax1.plot(profiles[j,:], color='b') + ax2.plot(fluxes[j,:], color='b') + + ax3.plot(err_history, color='k') + ax3.plot([j], [err_history[j]], 'bo') + ax3.set_yscale('log') + ax3.set_xlabel("Iteration") + ax3.set_ylabel("Residual") + plt.show() + + import sys + sys.exit(0) + + return {"converged": solver.reachedEnd, + "num_iterations": num_iterations, + "residual": solutionResidual, + "rms_error": solutionRmsError, + "err_history": err_history} + +def collect_statistics(inputs, num_samples=20): + iterations = [] + for i in range(num_samples): + result = solve_system(**inputs) + iterations.append(result["num_iterations"]) + + result["max_iterations"] = max(iterations) + result["min_iterations"] = min(iterations) + result["mean_iterations"] = np.mean(np.array(iterations)) + return result + +# ============================================================================== +# Settings +# ============================================================================== + +timescales = np.logspace(np.log10(0.01), np.log10(10), 7) +oscillation_amplitudes = [0.0]#, 1e-2, 1e-1] +noise_multipliers = [0.2] + +tolerance = 1e-1 +alpha_profile = 1 +alpha_flux = 0.1 + +shestakov_p = 3 +shestakov_q = -2 + +noise_amplitude = 0.0 # Amplitude of noise to add +noise_scalelength = 10.0 # Spatial scale length (number of cells) + +num_samples = 1 # Number of repeated solves, to gather statistics + +# ============================================================================== +# MAIN STARTS HERE +# ============================================================================== +tlog.setup() + +fig1, ax1 = plt.subplots() +fig2, ax2 = plt.subplots() +fig3, ax3 = plt.subplots() + +for oscillation_amplitude in oscillation_amplitudes: + for noise_multiplier in noise_multipliers: + mean_iterations = [] + max_iterations = [] + min_iterations = [] + for flux_timescale in timescales: + result = collect_statistics({"noise_timescale": flux_timescale * noise_multiplier, # AR time of random noise + "flux_timescale": flux_timescale, # Flux relaxation time + "oscillation_amplitude": oscillation_amplitude, + "noise_amplitude": noise_amplitude, # Amplitude of AR(1) noise + "noise_scalelength": noise_scalelength, # Spatial scale length (cells), + "p": shestakov_p, + "q": shestakov_q, + "EWMAParamTurbFlux": alpha_flux, + "EWMAParamProfile": alpha_profile, + "thetaParams": {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10}, + "tol": tolerance, + "maxIterations": 1000, + "plot_convergence": False}, + num_samples = num_samples) + + # Plot error history on figure 1 (for the last sample) + ax1.plot(result["err_history"], label=r"Time/$\tau_{{flux}}$ = {:.2f}".format(1./flux_timescale)) + + mean_iterations.append(result["mean_iterations"]) + max_iterations.append(result["max_iterations"]) + min_iterations.append(result["min_iterations"]) + + timescales = np.array(timescales) + + # Plot iterations vs timescale on figure 2 + ax2.plot( + 1./timescales, + mean_iterations, + "-o", + label=r"$\tau_{{noise}} / \tau_{{flux}} = {}, A_{{flux}} = {}$".format(noise_multiplier, oscillation_amplitude), + ) + + if num_samples > 1: + print(mean_iterations, max_iterations, min_iterations) + ax2.fill_between(1./timescales, min_iterations, max_iterations, alpha=0.5) + + # Plot total turbulence run time on figure 3 + ax3.plot( + 1. / timescales, + np.array(mean_iterations) / timescales, + "-o", + label=r"$\tau_{{noise}} / \tau_{{flux}} = {}, A_{{flux}} = {}$".format(noise_multiplier, oscillation_amplitude), + ) + + if num_samples > 1: + ax3.fill_between(1. / timescales, + np.array(min_iterations) / timescales, + np.array(max_iterations) / timescales, alpha=0.5) + +# Plot of error against iteration +ax1.set_yscale("log") +ax1.set_xlabel("iteration number") +ax1.set_ylabel("rms error") +ax1.legend() + +fig1.savefig("residual_history.pdf") +fig1.savefig("residual_history.png") + +# Plot of iteration count against timescale +ax2.set_xscale("log") +ax2.set_yscale("log") +ax2.set_xlabel(r"Simulation run time / $\tau_{{flux}}$") +ax2.set_ylabel("Iterations required") +ax2.legend() + +fig2.savefig("iteration_count.pdf") +fig2.savefig("iteration_count.png") + +ax3.set_xscale("log") +ax3.set_yscale("log") +ax3.set_xlabel(r"Simulation run time / $\tau_{{flux}}$") +ax3.set_ylabel("Total run time required") +ax3.legend() + +fig3.savefig("run_time.pdf") +fig3.savefig("run_time.png") + +plt.show() + +sys.exit(0) + +# ============================================================================== +# Settings +# ============================================================================== + +shestakov_p = 9 +shestakov_q = -6 + +flux_timescale = 0.2 +noise_multiplier = 0.2 +oscillation_amplitude = 0.1 + +tolerance = 1e-1 +alpha_ps = np.logspace(-2, 0, num=10) +alpha_ds = np.logspace(-2, 0, num=10) + +noise_amplitude = 0.1 # Amplitude of noise to add +noise_scalelength = 10.0 # Spatial scale length (number of cells) + +num_samples = 1 # Number of repeated solves, to gather statistics + +# ============================================================================== +# MAIN STARTS HERE +# ============================================================================== +tlog.setup() + +iterations = np.zeros((len(alpha_ps), len(alpha_ds))) + +for i, alpha_p in enumerate(alpha_ps): + for j, alpha_d in enumerate(alpha_ds): + result = collect_statistics({"noise_timescale": flux_timescale * noise_multiplier, # AR time of random noise + "flux_timescale": flux_timescale, # Flux relaxation time + "oscillation_amplitude": oscillation_amplitude, + "noise_amplitude": noise_amplitude, # Amplitude of AR(1) noise + "noise_scalelength": noise_scalelength, # Spatial scale length (cells) + "p": shestakov_p, + "q": shestakov_q, + "EWMAParamTurbFlux": alpha_d, + "EWMAParamProfile": alpha_p, + "thetaParams": {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10}, + "tol": tolerance, + "maxIterations": 1000, + "plot_convergence": False}, + num_samples = num_samples) + iterations[i,j] = result["mean_iterations"] + +plt.contourf(np.log10(alpha_ds), np.log10(alpha_ps), np.log10(iterations), 50) +plt.xlabel(r'Diffusion smoothing $\log_{10} \alpha_D$') +plt.ylabel(r'Profile smoothing $\log_{10} \alpha_p$') +plt.title("$\log_{10}$ iterations") +plt.colorbar() +plt.savefig("iterations_contour.png") +plt.show() diff --git a/tango/extras/fluxrelaxation.py b/tango/extras/fluxrelaxation.py new file mode 100644 index 0000000..b519333 --- /dev/null +++ b/tango/extras/fluxrelaxation.py @@ -0,0 +1,199 @@ +import numpy as np + + +class FluxRelaxation(object): + """Decorator that adds time dependence to fluxes + Simple relaxation on a fixed timescale + """ + + def __init__(self, timescale, fluxModel): + """ + Inputs: + fluxModel fluxmodel to be decorated. + Should have a get_flux(profiles) method which takes + a dictionary input and returns a dictionary of fluxes + timescale Ratio of flux relaxation timescale to coupling period + """ + assert hasattr(fluxModel, "get_flux") and callable( + getattr(fluxModel, "get_flux") + ) + assert timescale > 0.0 + + self.fluxModel = fluxModel + # Weight between 0 and 1 on last fluxes (-> 0 as timescale becomes shorter) + self.weight = np.exp(-1.0 / timescale) + self.lastFluxes = None # No previous flux + + def get_flux(self, profiles): + # Call the flux model to get the new flux + newFluxes = self.fluxModel.get_flux(profiles) + if self.lastFluxes is None: + self.lastFluxes = newFluxes + + # Apply relaxation to each flux channel + for key in newFluxes: + newFluxes[key] = ( + self.weight * self.lastFluxes[key] + + (1.0 - self.weight) * newFluxes[key] + ) + + self.lastFluxes = newFluxes + return newFluxes + + +class FluxDoubleRelaxation(object): + """Decorator that adds time dependence to fluxes, with two timescales: + one timescale for the drive, and one for the damping. + """ + + def __init__(self, turb_timescale, damp_timescale, fluxModel): + """ + Inputs: + fluxModel + fluxmodel to be decorated. + Should have a get_flux(profiles) method which takes + a dictionary input and returns a dictionary of fluxes + turb_timescale + Ratio of flux relaxation timescale to coupling period + damp_timescale + Ratio of damping timescale to coupling period + """ + assert hasattr(fluxModel, "get_flux") and callable( + getattr(fluxModel, "get_flux") + ) + assert turb_timescale > 0.0 + assert damp_timescale > 0.0 + + self.fluxModel = fluxModel + # Weight between 0 and 1 on last fluxes (-> 0 as timescale becomes shorter) + self.turb_weight = np.exp(-1.0 / turb_timescale) + self.damp_weight = np.exp(-1.0 / damp_timescale) + self.lastTurb = None # No previous turbulent drive or damping + self.lastDamp = None + + def get_flux(self, profiles): + # Call the flux model to get the new flux + newFluxes = self.fluxModel.get_flux(profiles) + if self.lastTurb is None: + # Shallow copies of the flux dictionary + self.lastTurb = newFluxes.copy() + self.lastDrive = newFluxes.copy() + + # Apply relaxation to each flux channel + for key in newFluxes: + # Relax both turbulence drive and damping towards new fluxes + self.lastTurb[key] = ( + self.turb_weight * self.lastTurb[key] + + (1.0 - self.turb_weight) * newFluxes[key] + ) + self.lastDrive[key] = ( + self.damp_weight * self.lastDrive[key] + + (1.0 - self.damp_weight) * newFluxes[key] + ) + # Calculate a ratio which goes towards the input flux at long time + newFluxes[key] = self.lastTurb[key] ** 2 / self.lastDrive[key] + + return newFluxes + +def FluxRelaxationOscillation(timescale, amplitude, fluxModel): + """Decorator that adds time dependence to fluxes with oscillation + Relaxation on a fixed timescale, with oscillation of specified relative magnitude + + # Inputs + timescale Ratio of flux relaxation timescale to coupling period + + amplitude Relative oscillation amplitude e.g. 0.2 is 20% + + fluxModel fluxmodel to be decorated. + Should have a get_flux(profiles) method which takes + a dictionary input and returns a dictionary of fluxes + + Note: The fluxModel object will be modified + + # Returns + + The modified fluxModel object, with new get_flux method + + # Example + + """ + assert hasattr(fluxModel, "get_flux") and callable( + getattr(fluxModel, "get_flux") + ) + assert timescale > 0.0 + + # We're going to wrap this function + inner_get_flux = fluxModel.get_flux + + # Weight between 0 and 1 on last fluxes (-> 0 as timescale becomes shorter) + weight = np.exp(-1.0 / timescale) + last_fluxes = None # No previous flux + time = 0 # Keeps track of time for the oscillation phase + + # Replacement flux calculation + def get_flux(profiles): + nonlocal inner_get_flux + nonlocal last_fluxes + nonlocal time + nonlocal weight + # Call the wrapped flux model to get the new flux + new_fluxes = inner_get_flux(profiles) + if last_fluxes is None: + last_fluxes = new_fluxes + + # Apply relaxation to each flux channel + for key in new_fluxes: + new_fluxes[key] = ( + weight * last_fluxes[key] + + (1.0 - weight) * new_fluxes[key] + ) + + last_fluxes = new_fluxes.copy() # Damping based on flux without oscillation + + # Add a relative oscillation + for key in new_fluxes: + new_fluxes[key] *= (1. + amplitude * np.sin(3. * time / timescale)) + time += 1 + return new_fluxes + + # Replace the get_flux method + fluxModel.get_flux = get_flux + return fluxModel + +def FluxAverage(nsteps, fluxModel): + """Decorator that averages the flux over a given number of iterations + + # Inputs: + fluxModel fluxmodel to be decorated. + Should have a get_flux(profiles) method which takes + a dictionary input and returns a dictionary of fluxes + nsteps Number of steps + + """ + assert hasattr(fluxModel, "get_flux") and callable( + getattr(fluxModel, "get_flux") + ) + assert nsteps > 0 + + inner_get_flux = fluxModel.get_flux + + def get_flux(profiles): + nonlocal inner_get_flux + + # Call the flux model to get the new flux + fluxes = inner_get_flux(profiles) + + # Sum fluxes over nsteps + for i in range(nsteps - 1): + next_fluxes = inner_get_flux(profiles) + for key in fluxes: + fluxes[key] += next_fluxes[key] + + # Divide by nsteps to get the average flux + for key in fluxes: + fluxes[key] /= nsteps + return fluxes + + # Replace the get_flux function + fluxModel.get_flux = get_flux + return fluxModel diff --git a/tango/extras/noisyflux.py b/tango/extras/noisyflux.py index 124567f..a21b7ec 100644 --- a/tango/extras/noisyflux.py +++ b/tango/extras/noisyflux.py @@ -15,7 +15,7 @@ class NoisyFlux(object): to the flux returned by an inherent fluxModel. The correlation length in space is controlled by the AR parameter. """ - def __init__(self, fluxModel, amplitude, tac_x, dx): + def __init__(self, amplitude, tac_x, dx, fluxModel): """ Inputs: fluxModel fluxmodel to be decorated @@ -66,7 +66,7 @@ class UniformNoisyFlux(object): """Decorator that adds random noise uniform in space (i.e., *perfectly* correlated) to the flux returned by an inherent fluxModel. """ - def __init__(self, fluxModel, amplitude): + def __init__(self, amplitude, fluxModel): """ Inputs: fluxModel fluxmodel to be decorated @@ -102,7 +102,64 @@ def _add_noise(v, amplitude): h = dampen_sides(h) noisy_v = (1 + h) * v return noisy_v - + + +class NoisyFluxSpaceTime(object): + """Decorator that adds random noise correlated in space and time, + generated by an AR(1) process, to the flux returned by an inherent + fluxModel. The correlation length in space is controlled by the + AR parameters + """ + def __init__(self, amplitude, tac_x, dx, tac_t, dt, fluxModel): + """ + Inputs: + fluxModel fluxmodel to be decorated + amplitude amplitude of noise (scalar) + tac_x autocorrelation time/length measured in units of x + dx grid spacing (scalar) + """ + self.fluxModel = fluxModel + self.amplitude = amplitude # scalar + self.tac_x = tac_x / dx # autocorrelation length measured in discrete units (does not have to be integer) + self.phi_t = np.exp(- dt / tac_t) # AR(1) parameter in time + print("Phi = {}".format(self.phi_t)) + self._last_noise = {} # Stores the noise from the last call, to provide correlation in time + + def get_flux(self, profiles): + fluxes = self.fluxModel.get_flux(profiles) + noisyFluxes = self._add_noise_to_fluxes(fluxes) + return noisyFluxes + + def _get_last_noise(self, label, numSamples): + """Returns the last noise for the given label. + If this is the first call, then generate noise""" + if label in self._last_noise: + return self._last_noise[label] + # No noise so generate + return self.amplitude * ar1noise(numSamples, self.tac_x) + + def _add_noise_to_fluxes(self, fluxes): + """Add noise to each field. The noise added to each field is statistically independent.""" + noisyFluxes = {} + for label in fluxes: + flux = fluxes[label] + numSamples = len(flux) + + noise = self.amplitude * ar1noise(numSamples, self.tac_x) # Correlated in space, uncorrelated in time + last_noise = self._get_last_noise(label, numSamples) + + # AR(1) process in time + # Note: Need to scale the noise so that the variance of the result is unchanged + noise = self.phi_t * last_noise + noise * np.sqrt(1 - self.phi_t**2) + + self._last_noise[label] = noise.copy() + + # Damp noise close to both boundaries + noise = dampen_sides(noise) + + noisyFluxes[label] = (1 + noise) * flux + return noisyFluxes + def ar1noise(numSamples, tac): """Generate an AR(1) process. @@ -123,7 +180,7 @@ def ar1noise(numSamples, tac): Outputs: noise noise with unit variance and autocorrelation time tac. (array of length numSamples) """ - lamb = (tac - 1) / (tac + 1) + lamb = np.exp(-1. / tac) cc = 1000 m = int( -np.log(cc) / (2*np.log(lamb)) ) + 1 numSamplesPlusm = numSamples + m @@ -159,4 +216,4 @@ def dampen_sides(v, numPts=None): v_dampened = v.copy() v_dampened[0:numPts] = v_dampened[0:numPts] * np.linspace(0, 1, numPts) v_dampened[-numPts:] = v_dampened[-numPts:] * np.linspace(1, 0, numPts) - return v_dampened \ No newline at end of file + return v_dampened diff --git a/tango/extras/shestakov_nonlinear_diffusion.py b/tango/extras/shestakov_nonlinear_diffusion.py index a57b183..508e235 100644 --- a/tango/extras/shestakov_nonlinear_diffusion.py +++ b/tango/extras/shestakov_nonlinear_diffusion.py @@ -3,104 +3,111 @@ from __future__ import division import numpy as np -"""Shestakov Test Module""" -class AnalyticFluxModel(object): - """Class-based interface to the flux model in the Shestakov analytic - test problem. - - Slight modification of boundary conditions so that n(1) = nL rather than n(1)=0, - in order to make the computation of the flux at the boundary much easier by - avoiding a divide by zero condition. - """ - def __init__(self, dx): +class AnalyticFluxModel: + """Flux model from Shestakov et al. (2003) + Return the flux Gamma, which depends on the density profile n as follows: + Gamma[n] = -(dn/dx)^p * n^q + """ + def __init__(self, dx, p = 3, q = -2, field='n'): + """ + # Inputs + dx : Grid spacing, used to calculate gradients + p : int + Flux depends on gradient to this power + """ self.dx = dx - + self.p = p + self.q = q + self.field = field + def get_flux(self, profiles): - n = profiles['n'] + """Return flux as a dictionary, using profiles dictionary + """ + n = profiles[self.field] flux = {} - flux['n'] = get_flux(n, self.dx) + flux[self.field] = self.calc_flux(n, self.dx) return flux -#class shestakov_analytic_fluxmodel(object): -# """Class-based interface to the flux model in the Shestakov analytic -# test problem. -# -# Slight modification of boundary conditions so that n(1) = nL rather than n(1)=0, -# in order to make the computation of the flux at the boundary much easier by -# avoiding a divide by zero condition. -# -# Alias to AnalyticFluxModel which conforms to style guidleines. This class is kept around for backwards compatibility reasons. -# """ -# def __init__(self, dx): -# self.dx = dx -# -# def get_flux(self, n): -# return get_flux(n, self.dx) - - - -#============================================================================== -# Functions specifying the Shestakov analytic test problem start here -#============================================================================== -def H7contrib_Source(x): - S = GetSource(x) - H7 = S - return H7 + def calc_flux(self, n, dx): + """Return flux Gamma on the same grid as n + """ + dndx = self._dxCenteredDifference(n, dx) + return - dndx**self.p * n**self.q + def calc_flux_gradient(self, n, dx): + """ + Return divergence of flux + """ + flux = self.get_flux(n, dx) + return self._dxCenteredDifference(flux, dx) -def get_flux(n, dx): + def label(self): + return self.field + + def _dxCenteredDifference(self, u, dx): + """Compute du/dx. + du/dx is computed using centered differences on the same grid as u. For the edge points, one-point differences are used. + + Inputs: + u profile (array) + dx grid spacing (scalar) + + Outputs: + dudx (array, same length as u) + """ + dudx = np.zeros_like(u) + dudx[0] = (u[1] - u[0]) / dx + dudx[1:-1] = (u[2:] - u[:-2]) / (2*dx) + dudx[-1] = (u[-1] - u[-2]) / dx + return dudx + +class ShestakovTestProblem(AnalyticFluxModel): """Test problem from Shestakov et al. (2003) - Return the flux Gamma, which depends on the density profile n as follows: - Gamma[n] = -(dn/dx)^3 / n^2 - """ - Gamma = np.zeros_like(n) - - # Return flux Gamma on the same grid as n - dndx = _dxCenteredDifference(n, dx) - Gamma = - dndx**3 / n**2 - return Gamma - -def GetSource(x, S0=1, delta=0.1): - """Test problem from Shestakov et al. (2003). - Return the source S.""" - S = np.zeros_like(x) - S[x < delta] = S0 - return S - -def _dxCenteredDifference(u, dx): - """Compute du/dx. - du/dx is computed using centered differences on the same grid as u. For the edge points, one-point differences are used. - - Inputs: - u profile (array) - dx grid spacing (scalar) - - Outputs: - dudx (array, same length as u) - """ - dudx = np.zeros_like(u) - dudx[0] = (u[1] - u[0]) / dx - dudx[1:-1] = (u[2:] - u[:-2]) / (2*dx) - dudx[-1] = (u[-1] - u[-2]) / dx - return dudx - -def steady_state_solution(x, nL, S0=1, delta=0.1): - """Return the exast steady state solution for the Shestakov test problem - - Inputs: - x Spatial coordinate grid (array) - nL boundary condition n(L) (scalar) - S0 parameter in source term --- amplitude (scalar) - delta parameter in source term --- location where it turns off (scalar) - Outputs: - """ - nright = ( nL**(1/3) + 1/3 * (S0 * delta)**(1/3) *(1-x) )**3 - # nleft = (L - delta + 0.75*(delta - x**(4/3) / delta**(1/3)))**3- - nleft = ( nL**(1/3) + 1/3 * (S0 * delta)**(1/3) * (1 - delta + (3/4) * (delta - x**(4/3) / delta**(1/3))))**3 - nss = nright - nss[x < delta] = nleft[x < delta] - return nss + Return the flux Gamma, which depends on the density profile n as follows: + Gamma[n] = -(dn/dx)^p * n^q + """ + def __init__(self, dx, p = 3, q = -2, S0 = 1, delta = 0.1): + """ + # Inputs + dx Grid spacing + S0 parameter in source term --- amplitude (scalar) + delta parameter in source term --- location where it turns off (scalar) + """ + super(ShestakovTestProblem, self).__init__(dx, p=p, q=q, field='n') + self.S0 = S0 + self.delta = delta + + def H7contrib_Source(self, x): + return self.GetSource(x) + + def GetSource(self, x): + """Test problem from Shestakov et al. (2003). + Return the source S.""" + S = np.zeros_like(x) + S[x < self.delta] = self.S0 + return S + + def steady_state_solution(self, x, nL): + """Return the exact steady state solution for the Shestakov test problem + Generalised for arbitrary p, q + + Inputs: + x Spatial coordinate grid (array) + nL boundary condition n(L) (scalar) + Outputs: + """ + + coef = self.q / self.p + 1. # This is 1/3 for p=3,q=-2 standard case + L = 1.0 + + # Solution in region delta < x < L + nright = (nL**coef + coef * (self.S0 * self.delta)**(1./self.p) * (L - x))**(1./coef) + nleft = (nL**coef + coef * (self.S0 * self.delta)**(1./self.p) * (L - self.delta + (self.p/(self.p + 1)) * (self.delta - x * (x / self.delta)**(1./self.p))))**(1./coef) + + nss = nright + nss[x < self.delta] = nleft[x < self.delta] + return nss -def GetSteadyStateSolution(x, nL, S0=1, delta=0.1): - """alias to steady_state_solution which conforms to style guidelines, but this function is kept around for backwards compatibility.""" - return steady_state_solution(x, nL, S0=S0, delta=delta) \ No newline at end of file + def GetSteadyStateSolution(self, x, nL): + """alias to steady_state_solution which conforms to style guidelines, but this function is kept around for backwards compatibility.""" + return self.steady_state_solution(x, nL) diff --git a/tango/lodestro_method.py b/tango/lodestro_method.py index 53cefb9..5105304 100644 --- a/tango/lodestro_method.py +++ b/tango/lodestro_method.py @@ -377,7 +377,7 @@ def turbflux_to_Hcoeffs_multifield(self, fields, profiles): extradataAllFields[label] = { 'D': D, 'c': c, 'profileTurbGrid': profilesTurbGrid[label], 'profileEWMATurbGrid': profilesEWMATurbGrid[label], - 'fluxTurbGrid': fluxesTurbGrid[label], 'smoothedFluxTurbGrid': smoothedFluxTurbGrid, + 'fluxTurbGrid': fluxesTurbGrid[label], 'smoothedFluxTurbGrid': smoothedFluxTurbGrid, 'fluxEWMATurbGrid':fluxesTurbGrid[label], 'DTurbGrid': DTurbGrid, 'cTurbGrid': cTurbGrid, 'DEWMATurbGrid': DEWMATurbGrid, 'cEWMATurbGrid': cEWMATurbGrid, 'DHatTurbGrid': DcDataTurbGrid['DHat'], 'cHatTurbGrid': DcDataTurbGrid['cHat'], 'thetaTurbGrid': DcDataTurbGrid['theta']} @@ -742,4 +742,4 @@ def get_x_transport_grid(self): def get_x_turbulence_grid(self): return self.x - \ No newline at end of file + diff --git a/tango/solver.py b/tango/solver.py index 4a69523..0552d0b 100644 --- a/tango/solver.py +++ b/tango/solver.py @@ -116,6 +116,7 @@ def __init__(self, L, x, tArray, maxIterations, tol, compute_all_H_all_fields, f self.reachedEnd = False self.profilesAllIterations = None self.fluxesAllIterations = None + self.fluxesAllIterationsEWMA = None # initialize data containers for storing profiles and fluxes when using multiple timesteps self.profilesAllTimesteps = {} @@ -140,9 +141,11 @@ def take_timestep(self): self.errHistory[:] = 0 self.profilesAllIterations = {} self.fluxesAllIterations = {} + self.fluxesAllIterationsEWMA = {} for field in self.fields: self.profilesAllIterations[field.label] = np.zeros((self.maxIterations, self.N)) self.fluxesAllIterations[field.label] = np.zeros((self.maxIterations, self.N)) + self.fluxesAllIterationsEWMA[field.label] = np.zeros((self.maxIterations, self.N)) tango_logging.info("Timestep m={}: Beginning iteration loop ...".format(self.m)) @@ -165,7 +168,8 @@ def take_timestep(self): for field in self.fields: self.profilesAllIterations[field.label] = self.profilesAllIterations[field.label][0:self.countStoredIterations, :] self.fluxesAllIterations[field.label] = self.fluxesAllIterations[field.label][0:self.countStoredIterations, :] - + self.fluxesAllIterationsEWMA[field.label] = self.fluxesAllIterationsEWMA[field.label][0:self.countStoredIterations, :] + # save the last iteration for the profiles and fluxes to represent the converged values at each timestep for field in self.fields: self.profilesAllTimesteps[field.label][self.m, :] = self.profiles[field.label] @@ -257,6 +261,7 @@ def compute_next_iteration(self): if self.saveFluxesInMemory: if extradataAllFields is not None: self.fluxesAllIterations[field.label][index, :] = extradataAllFields[field.label]['fluxTurbGrid'] + self.fluxesAllIterationsEWMA[field.label][index, :] = extradataAllFields[field.label]['fluxEWMATurbGrid'] self.fileHandlerExecutor.execute_scheduled(datadict, self.iterationNumber) @@ -405,14 +410,19 @@ def is_unacceptable(self, oldProfiles, newProfiles): def ok(self): """True unless a stop condition is reached.""" if self.solutionError == True: + print("solutionError") return False elif self.m >= len(self.tArray) - 1: + print("m") return False elif self.t >= self.tFinal: + print("tFinal") return False elif self.iterationNumber >= self.maxIterations: + print("iterationNumber") return False elif self.countStoredIterations >= self.maxIterationsPerSet: + print("countStoredIterations") return False else: return True