From 265b2038c9293b56ddc453c82d797ff9261d0fd8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 14 Feb 2022 15:01:31 -0800 Subject: [PATCH 1/8] Adding flux relaxation decorators Single and double timescale relaxation decorators. These can wrap other decorators to introduce a timescale. --- examples/shestakov_time_dependent.py | 167 +++++++++++++++++++++++++++ tango/extras/fluxrelaxation.py | 96 +++++++++++++++ 2 files changed, 263 insertions(+) create mode 100644 examples/shestakov_time_dependent.py create mode 100644 tango/extras/fluxrelaxation.py diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py new file mode 100644 index 0000000..086fe56 --- /dev/null +++ b/examples/shestakov_time_dependent.py @@ -0,0 +1,167 @@ +"""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 import shestakov_nonlinear_diffusion +from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation +import tango + + +def initialize_shestakov_problem(): + # Problem Setup + L = 1 # size of domain + N = 500 # number of spatial grid points + dx = L / (N - 1) # spatial grid size + x = np.arange(N) * dx # location corresponding to grid points j=0, ..., N-1 + nL = 1e-2 # right boundary condition + nInitialCondition = 1 - 0.5 * x + return (L, N, dx, x, nL, nInitialCondition) + + +def initialize_parameters(): + maxIterations = 5000 + thetaParams = {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10} + EWMAParamTurbFlux = 0.3 + EWMAParamProfile = 1 + lmParams = { + "EWMAParamTurbFlux": EWMAParamTurbFlux, + "EWMAParamProfile": EWMAParamProfile, + "thetaParams": thetaParams, + } + tol = 1e-11 # tol for convergence... reached when a certain error < tol + return (maxIterations, lmParams, tol) + + +class ComputeAllH(object): + def __init__(self): + pass + + def __call__(self, t, x, profiles, HCoeffsTurb): + # n = profiles['default'] + # Define the contributions to the H coefficients for the Shestakov Problem + H1 = np.ones_like(x) + H7 = shestakov_nonlinear_diffusion.H7contrib_Source(x) + + HCoeffs = tango.multifield.HCoefficients(H1=H1, H7=H7) + HCoeffs = HCoeffs + HCoeffsTurb + + return HCoeffs + + +# ============================================================================== +# MAIN STARTS HERE +# ============================================================================== +tlog.setup() + +plt.figure() + +timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] + +for damping_multiplier in [1.0, 2.0, 5.0]: + iterations = [] + for turb_timescale in timescales: + + tlog.info("Initializing...") + L, N, dx, x, nL, n = initialize_shestakov_problem() + maxIterations, lmParams, tol = initialize_parameters() + fluxModel = FluxDoubleRelaxation( + shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), + turb_timescale, + damping_multiplier * turb_timescale, + ) + + label = "n" + turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) + + compute_all_H_density = ComputeAllH() + lodestroMethod = tango.lodestro_method.lm( + lmParams["EWMAParamTurbFlux"], + lmParams["EWMAParamProfile"], + lmParams["thetaParams"], + ) + field0 = tango.multifield.Field( + label=label, + rightBC=nL, + profile_mminus1=n, + compute_all_H=compute_all_H_density, + 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( + L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields + ) + + tlog.info("Initialization complete.") + tlog.info("Entering main time loop...") + + while solver.ok: + # Implicit time advance: iterate to solve the nonlinear equation! + solver.take_timestep() + + n = solver.profiles[label] # finished solution + nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) + + solutionResidual = (n - nss) / np.max(np.abs(nss)) + solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) + + if solver.reachedEnd == True: + print("The solution has been reached successfully.") + print( + "Error compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + else: + print("The solver failed for some reason.") + print( + "Error at end compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + + iterations.append(len(solver.errHistoryFinal)) + + # plt.plot(solver.errHistoryFinal, label = str(turb_timescale)) + + # plt.legend() + # plt.yscale('log') + # plt.xlabel('iteration number') + # plt.ylabel('rms error') + # plt.savefig('residual_history.pdf') + # plt.savefig('residual_history.png') + # plt.figure() + + plt.plot( + timescales, + iterations, + "-o", + label=r"$\tau_{{damp}} / \tau_{{turb}} = {}$".format(damping_multiplier), + ) + plt.axvline(1.0 / damping_multiplier, linestyle="--", color="k") + +plt.xscale("log") +plt.yscale("log") +plt.xlabel("Turbulence relaxation time") +plt.ylabel("Iterations required") +plt.legend() + +plt.savefig("iteration_count.pdf") +plt.savefig("iteration_count.png") + +plt.show() diff --git a/tango/extras/fluxrelaxation.py b/tango/extras/fluxrelaxation.py new file mode 100644 index 0000000..56b1878 --- /dev/null +++ b/tango/extras/fluxrelaxation.py @@ -0,0 +1,96 @@ +import numpy as np + + +class FluxRelaxation(object): + """Decorator that adds time dependence to fluxes + Simple relaxation on a fixed timescale + """ + + def __init__(self, fluxModel, timescale): + """ + 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, fluxModel, turb_timescale, damp_timescale): + """ + 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 From d2cf4baca25e8625f07d3dec90ae67f05d6f98a5 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 14 Feb 2022 17:22:42 -0800 Subject: [PATCH 2/8] Add noise to time-dependent Shestakov case Use AR noise decorator --- examples/shestakov_time_dependent.py | 207 ++++++++++++++++----------- 1 file changed, 121 insertions(+), 86 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 086fe56..0523a40 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -14,6 +14,7 @@ import tango.tango_logging as tlog from tango.extras import shestakov_nonlinear_diffusion from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation +from tango.extras.noisyflux import NoisyFlux import tango @@ -38,7 +39,7 @@ def initialize_parameters(): "EWMAParamProfile": EWMAParamProfile, "thetaParams": thetaParams, } - tol = 1e-11 # tol for convergence... reached when a certain error < tol + tol = 1e-2 # tol for convergence... reached when a certain error < tol return (maxIterations, lmParams, tol) @@ -58,110 +59,144 @@ def __call__(self, t, x, profiles, HCoeffsTurb): return HCoeffs +# ============================================================================== +# Settings +# ============================================================================== + +timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] +damping_multipliers = [1.0, 2.0, 5.0] + +noise_amplitude = 1e-3 # Amplitude of noise to add +noise_scalelength = 10.0 # Spatial scale length (number of cells) + +num_samples = 20 # Number of repeated solves, to gather statistics + # ============================================================================== # MAIN STARTS HERE # ============================================================================== tlog.setup() -plt.figure() - -timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] +fig1, ax1 = plt.subplots() +fig2, ax2 = plt.subplots() -for damping_multiplier in [1.0, 2.0, 5.0]: - iterations = [] +for damping_multiplier in damping_multipliers: + mean_iterations = [] + max_iterations = [] + min_iterations = [] for turb_timescale in timescales: + its = [] # Iterations for these settings + for sample in range(num_samples): + tlog.info("Initializing...") + L, N, dx, x, nL, n = initialize_shestakov_problem() + maxIterations, lmParams, tol = initialize_parameters() + fluxModel = NoisyFlux( + FluxDoubleRelaxation( + shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), + turb_timescale, + damping_multiplier * turb_timescale, + ), + noise_amplitude, + noise_scalelength, + 1.0, + ) + + label = "n" + turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) - tlog.info("Initializing...") - L, N, dx, x, nL, n = initialize_shestakov_problem() - maxIterations, lmParams, tol = initialize_parameters() - fluxModel = FluxDoubleRelaxation( - shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), - turb_timescale, - damping_multiplier * turb_timescale, - ) - - label = "n" - turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) - - compute_all_H_density = ComputeAllH() - lodestroMethod = tango.lodestro_method.lm( - lmParams["EWMAParamTurbFlux"], - lmParams["EWMAParamProfile"], - lmParams["thetaParams"], - ) - field0 = tango.multifield.Field( - label=label, - rightBC=nL, - profile_mminus1=n, - compute_all_H=compute_all_H_density, - 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( - L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields - ) - - tlog.info("Initialization complete.") - tlog.info("Entering main time loop...") - - while solver.ok: - # Implicit time advance: iterate to solve the nonlinear equation! - solver.take_timestep() - - n = solver.profiles[label] # finished solution - nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) - - solutionResidual = (n - nss) / np.max(np.abs(nss)) - solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) - - if solver.reachedEnd == True: - print("The solution has been reached successfully.") - print( - "Error compared to analytic steady state solution is %f" - % (solutionRmsError) + compute_all_H_density = ComputeAllH() + lodestroMethod = tango.lodestro_method.lm( + lmParams["EWMAParamTurbFlux"], + lmParams["EWMAParamProfile"], + lmParams["thetaParams"], ) - else: - print("The solver failed for some reason.") - print( - "Error at end compared to analytic steady state solution is %f" - % (solutionRmsError) + field0 = tango.multifield.Field( + label=label, + rightBC=nL, + profile_mminus1=n, + compute_all_H=compute_all_H_density, + lodestroMethod=lodestroMethod, ) + fields = [field0] + tango.multifield.check_fields_initialize(fields) - iterations.append(len(solver.errHistoryFinal)) + compute_all_H_all_fields = tango.multifield.ComputeAllHAllFields( + fields, turbHandler + ) - # plt.plot(solver.errHistoryFinal, label = str(turb_timescale)) + tArray = np.array([0, 1e4]) # specify the timesteps to be used. - # plt.legend() - # plt.yscale('log') - # plt.xlabel('iteration number') - # plt.ylabel('rms error') - # plt.savefig('residual_history.pdf') - # plt.savefig('residual_history.png') - # plt.figure() + solver = tango.solver.Solver( + L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields + ) - plt.plot( + tlog.info("Initialization complete.") + tlog.info("Entering main time loop...") + + while solver.ok: + # Implicit time advance: iterate to solve the nonlinear equation! + solver.take_timestep() + + n = solver.profiles[label] # finished solution + nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) + + solutionResidual = (n - nss) / np.max(np.abs(nss)) + solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) + + if solver.reachedEnd == True: + print("The solution has been reached successfully.") + print( + "Error compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + else: + print("The solver failed for some reason.") + print( + "Error at end compared to analytic steady state solution is %f" + % (solutionRmsError) + ) + + its.append(len(solver.errHistoryFinal)) + + # Plot error history on figure 1 (for the last sample) + ax1.plot(solver.errHistoryFinal, label=str(turb_timescale)) + + its = np.array(its) + mean_iterations.append(np.mean(its)) + max_iterations.append(np.amax(its)) + min_iterations.append(np.amin(its)) + print(its) + + # Plot iterations vs timescale on figure 2 + ax2.plot( timescales, - iterations, + mean_iterations, "-o", label=r"$\tau_{{damp}} / \tau_{{turb}} = {}$".format(damping_multiplier), ) - plt.axvline(1.0 / damping_multiplier, linestyle="--", color="k") -plt.xscale("log") -plt.yscale("log") -plt.xlabel("Turbulence relaxation time") -plt.ylabel("Iterations required") -plt.legend() + if num_samples > 1: + print(mean_iterations, max_iterations, min_iterations) + ax2.fill_between(timescales, min_iterations, max_iterations, alpha=0.5) + + ax2.axvline(1.0 / damping_multiplier, linestyle="--", color="k") + +# 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("Turbulence relaxation time") +ax2.set_ylabel("Iterations required") +ax2.legend() -plt.savefig("iteration_count.pdf") -plt.savefig("iteration_count.png") +fig2.savefig("iteration_count.pdf") +fig2.savefig("iteration_count.png") plt.show() From e493b67c6740eb4e085e34f4421306f89d5b16ac Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 11 Jul 2022 11:31:54 -0700 Subject: [PATCH 3/8] Noisy shestakov flux July 2022 Results don't agree with experience with GENE-Tango --- examples/shestakov_time_dependent.py | 397 +++++++++++++++++++-------- tango/extras/noisyflux.py | 63 ++++- tango/lodestro_method.py | 4 +- tango/solver.py | 12 +- 4 files changed, 351 insertions(+), 125 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 0523a40..8c5cec3 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -14,10 +14,9 @@ import tango.tango_logging as tlog from tango.extras import shestakov_nonlinear_diffusion from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation -from tango.extras.noisyflux import NoisyFlux +from tango.extras.noisyflux import NoisyFlux, NoisyFluxSpaceTime import tango - def initialize_shestakov_problem(): # Problem Setup L = 1 # size of domain @@ -28,21 +27,6 @@ def initialize_shestakov_problem(): nInitialCondition = 1 - 0.5 * x return (L, N, dx, x, nL, nInitialCondition) - -def initialize_parameters(): - maxIterations = 5000 - thetaParams = {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10} - EWMAParamTurbFlux = 0.3 - EWMAParamProfile = 1 - lmParams = { - "EWMAParamTurbFlux": EWMAParamTurbFlux, - "EWMAParamProfile": EWMAParamProfile, - "thetaParams": thetaParams, - } - tol = 1e-2 # tol for convergence... reached when a certain error < tol - return (maxIterations, lmParams, tol) - - class ComputeAllH(object): def __init__(self): pass @@ -58,127 +42,244 @@ def __call__(self, t, x, profiles, HCoeffsTurb): return HCoeffs - # ============================================================================== -# Settings +# Solve function # ============================================================================== -timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] -damping_multipliers = [1.0, 2.0, 5.0] - -noise_amplitude = 1e-3 # Amplitude of noise to add -noise_scalelength = 10.0 # Spatial scale length (number of cells) - -num_samples = 20 # Number of repeated solves, to gather statistics +def solve_system(noise_timescale = 1.0, # AR time of random noise + turb_timescale = 1.0, # Flux relaxation time + damping_timescale = 1.0, # Flux damping timescale + noise_amplitude = 0.0, # Amplitude of AR(1) noise + noise_scalelength = 10, # Spatial scale length (cells) + EWMAParamTurbFlux = 0.01, + EWMAParamProfile = 1.0, + thetaParams = {"Dmin": 1e-5, "Dmax": 1e13, "dpdxThreshold": 10}, + tol = 1e-2, + maxIterations = 1000, + plot_convergence = False, + check_tol = 0.1 + ): + tlog.info("Initializing...") + L, N, dx, x, nL, n = initialize_shestakov_problem() + fluxModel = NoisyFluxSpaceTime( + FluxDoubleRelaxation( + shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), + turb_timescale, + damping_timescale, + ), + noise_amplitude, + noise_scalelength, + 1.0, # dx + noise_timescale, # noise timescale + 1.0 # dt + ) -# ============================================================================== -# MAIN STARTS HERE -# ============================================================================== -tlog.setup() + label = "n" + turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) -fig1, ax1 = plt.subplots() -fig2, ax2 = plt.subplots() + compute_all_H_density = ComputeAllH() + lodestroMethod = tango.lodestro_method.lm( + EWMAParamTurbFlux, + EWMAParamProfile, + thetaParams, + ) + field0 = tango.multifield.Field( + label=label, + rightBC=nL, + profile_mminus1=n, + compute_all_H=compute_all_H_density, + lodestroMethod=lodestroMethod, + ) + fields = [field0] + tango.multifield.check_fields_initialize(fields) -for damping_multiplier in damping_multipliers: - mean_iterations = [] - max_iterations = [] - min_iterations = [] - for turb_timescale in timescales: - its = [] # Iterations for these settings - for sample in range(num_samples): - tlog.info("Initializing...") - L, N, dx, x, nL, n = initialize_shestakov_problem() - maxIterations, lmParams, tol = initialize_parameters() - fluxModel = NoisyFlux( - FluxDoubleRelaxation( - shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), - turb_timescale, - damping_multiplier * turb_timescale, - ), - noise_amplitude, - noise_scalelength, - 1.0, - ) + compute_all_H_all_fields = tango.multifield.ComputeAllHAllFields( + fields, turbHandler + ) - label = "n" - turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) + tArray = np.array([0, 1e4]) # specify the timesteps to be used. - compute_all_H_density = ComputeAllH() - lodestroMethod = tango.lodestro_method.lm( - lmParams["EWMAParamTurbFlux"], - lmParams["EWMAParamProfile"], - lmParams["thetaParams"], - ) - field0 = tango.multifield.Field( - label=label, - rightBC=nL, - profile_mminus1=n, - compute_all_H=compute_all_H_density, - lodestroMethod=lodestroMethod, - ) - fields = [field0] - tango.multifield.check_fields_initialize(fields) + solver = tango.solver.Solver( + L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields, + saveFluxesInMemory = plot_convergence + ) - compute_all_H_all_fields = tango.multifield.ComputeAllHAllFields( - fields, turbHandler + tlog.info("Initialization complete.") + tlog.info("Entering main time loop...") + + err_history = [] + profile_history = [] + flux_history = [] + num_iterations = 0 + restarts = [] + while num_iterations < maxIterations: + while solver.ok: + # Implicit time advance: iterate to solve the nonlinear equation! + solver.take_timestep() + + n = solver.profiles[label] # finished solution + nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) + + 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) ) - tArray = np.array([0, 1e4]) # specify the timesteps to be used. + if solutionRmsError > check_tol: + print("False convergence!") + solver.m = 0 # Reset timestep index + solver.t = 0.0 - solver = tango.solver.Solver( - L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields + 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) ) + break # Give up + + 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 +# ============================================================================== - tlog.info("Initialization complete.") - tlog.info("Entering main time loop...") +timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] +damping_multipliers = [] #[1.0, 2.0, 5.0, 10] +noise_multipliers = [1e-2, 1e-1, 1.0] - while solver.ok: - # Implicit time advance: iterate to solve the nonlinear equation! - solver.take_timestep() +tolerance = 1e-2 +alpha_profile = 1 +alpha_flux = 0.05 - n = solver.profiles[label] # finished solution - nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) +noise_amplitude = 0.01 # Amplitude of noise to add +noise_scalelength = 10.0 # Spatial scale length (number of cells) - solutionResidual = (n - nss) / np.max(np.abs(nss)) - solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) +num_samples = 20 # Number of repeated solves, to gather statistics - if solver.reachedEnd == True: - print("The solution has been reached successfully.") - print( - "Error compared to analytic steady state solution is %f" - % (solutionRmsError) - ) - else: - print("The solver failed for some reason.") - print( - "Error at end compared to analytic steady state solution is %f" - % (solutionRmsError) - ) - - its.append(len(solver.errHistoryFinal)) - - # Plot error history on figure 1 (for the last sample) - ax1.plot(solver.errHistoryFinal, label=str(turb_timescale)) - - its = np.array(its) - mean_iterations.append(np.mean(its)) - max_iterations.append(np.amax(its)) - min_iterations.append(np.amin(its)) - print(its) - - # Plot iterations vs timescale on figure 2 - ax2.plot( - timescales, - mean_iterations, - "-o", - label=r"$\tau_{{damp}} / \tau_{{turb}} = {}$".format(damping_multiplier), - ) +# ============================================================================== +# MAIN STARTS HERE +# ============================================================================== +tlog.setup() - if num_samples > 1: - print(mean_iterations, max_iterations, min_iterations) - ax2.fill_between(timescales, min_iterations, max_iterations, alpha=0.5) +fig1, ax1 = plt.subplots() +fig2, ax2 = plt.subplots() +fig3, ax3 = plt.subplots() - ax2.axvline(1.0 / damping_multiplier, linestyle="--", color="k") +for damping_multiplier in damping_multipliers: + for noise_multiplier in noise_multipliers: + mean_iterations = [] + max_iterations = [] + min_iterations = [] + for turb_timescale in timescales: + result = collect_statistics({"noise_timescale": turb_timescale * noise_multiplier, # AR time of random noise + "turb_timescale": turb_timescale, # Flux relaxation time + "damping_timescale": turb_timescale * damping_multiplier, # Flux damping timescale + "noise_amplitude": noise_amplitude, # Amplitude of AR(1) noise + "noise_scalelength": noise_scalelength, # Spatial scale length (cells) + "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=str(turb_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_{{damp}} / \tau_{{turb}} = {}, \tau_{{noise}} / \tau_{{turb}} = {}$".format(damping_multiplier, noise_multiplier), + ) + + if num_samples > 1: + print(mean_iterations, max_iterations, min_iterations) + ax2.fill_between(1./timescales, min_iterations, max_iterations, alpha=0.5) + + #ax2.axvline(1.0 / damping_multiplier, linestyle="--", color="k") + #ax2.axvline(1.0 / noise_multiplier, linestyle=":", color='k') + + # Plot total turbulence run time on figure 3 + ax3.plot( + 1. / timescales, + np.array(mean_iterations) / timescales, + "-o", + label=r"$\tau_{{damp}} / \tau_{{turb}} = {}, \tau_{{noise}} / \tau_{{turb}} = {}$".format(damping_multiplier, noise_multiplier), + ) + + 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") @@ -192,11 +293,69 @@ def __call__(self, t, x, profiles, HCoeffsTurb): # Plot of iteration count against timescale ax2.set_xscale("log") ax2.set_yscale("log") -ax2.set_xlabel("Turbulence relaxation time") +ax2.set_xlabel(r"Simulation run time / $\tau_{{turb}}$") 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_{{turb}}$") +ax3.set_ylabel("Total run time required") +ax3.legend() + +fig3.savefig("run_time.pdf") +fig3.savefig("run_time.png") + +plt.show() + + +# ============================================================================== +# Settings +# ============================================================================== + +turb_timescale = 10.0 +damping_multiplier = 1.0 +noise_multiplier = 1.0 + +tolerance = 1e-2 +alpha_ps = np.logspace(-1, 0, num=10) +alpha_ds = np.logspace(-2.5, -1.5, num=10) + +noise_amplitude = 0.01 # 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": turb_timescale * noise_multiplier, # AR time of random noise + "turb_timescale": turb_timescale, # Flux relaxation time + "damping_timescale": turb_timescale * damping_multiplier, # Flux damping timescale + "noise_amplitude": noise_amplitude, # Amplitude of AR(1) noise + "noise_scalelength": noise_scalelength, # Spatial scale length (cells) + "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/noisyflux.py b/tango/extras/noisyflux.py index 124567f..9bf7814 100644 --- a/tango/extras/noisyflux.py +++ b/tango/extras/noisyflux.py @@ -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, fluxModel, amplitude, tac_x, dx, tac_t, dt): + """ + 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/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 From 07ac28b53d8f37f6bb48cbd168989953ad0583f9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 11 Jul 2022 14:18:55 -0700 Subject: [PATCH 4/8] Separate AnalyticFluxModel and ShestakovTestProblem AnalyticFluxModel contains the flux calculation, for general p and q values. ShestakovTestProblem fixes p=3, q=-2, defines a source and has an analytic solution. --- examples/shestakov_time_dependent.py | 17 +- tango/extras/shestakov_nonlinear_diffusion.py | 181 +++++++++--------- 2 files changed, 99 insertions(+), 99 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 8c5cec3..86e56b7 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -28,14 +28,14 @@ def initialize_shestakov_problem(): return (L, N, dx, x, nL, nInitialCondition) class ComputeAllH(object): - def __init__(self): - pass + def __init__(self, test_problem): + self.test_problem = test_problem def __call__(self, t, x, profiles, HCoeffsTurb): # n = profiles['default'] # Define the contributions to the H coefficients for the Shestakov Problem H1 = np.ones_like(x) - H7 = shestakov_nonlinear_diffusion.H7contrib_Source(x) + H7 = self.test_problem.H7contrib_Source(x) HCoeffs = tango.multifield.HCoefficients(H1=H1, H7=H7) HCoeffs = HCoeffs + HCoeffsTurb @@ -61,9 +61,11 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise ): tlog.info("Initializing...") L, N, dx, x, nL, n = initialize_shestakov_problem() + test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx) + fluxModel = NoisyFluxSpaceTime( FluxDoubleRelaxation( - shestakov_nonlinear_diffusion.AnalyticFluxModel(dx), + test_problem, turb_timescale, damping_timescale, ), @@ -77,7 +79,6 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise label = "n" turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) - compute_all_H_density = ComputeAllH() lodestroMethod = tango.lodestro_method.lm( EWMAParamTurbFlux, EWMAParamProfile, @@ -87,7 +88,7 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise label=label, rightBC=nL, profile_mminus1=n, - compute_all_H=compute_all_H_density, + compute_all_H=ComputeAllH(test_problem), lodestroMethod=lodestroMethod, ) fields = [field0] @@ -118,7 +119,9 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise solver.take_timestep() n = solver.profiles[label] # finished solution - nss = shestakov_nonlinear_diffusion.steady_state_solution(x, nL) + + source = test_problem.GetSource(x) + nss = test_problem.steady_state_solution(x, nL) solutionResidual = (n - nss) / np.max(np.abs(nss)) solutionRmsError = np.sqrt(1 / len(n) * np.sum(solutionResidual ** 2)) diff --git a/tango/extras/shestakov_nonlinear_diffusion.py b/tango/extras/shestakov_nonlinear_diffusion.py index a57b183..f8c9f3f 100644 --- a/tango/extras/shestakov_nonlinear_diffusion.py +++ b/tango/extras/shestakov_nonlinear_diffusion.py @@ -3,104 +3,101 @@ 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 _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) -def get_flux(n, dx): + 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)^3 / n^2 + """ + def __init__(self, dx, 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=3, q=-2, 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 exast steady state solution for the Shestakov test problem + + Inputs: + x Spatial coordinate grid (array) + nL boundary condition n(L) (scalar) + Outputs: + """ + nright = ( nL**(1/3) + 1/3 * (self.S0 * self.delta)**(1/3) *(1-x) )**3 + nleft = ( nL**(1/3) + 1/3 * (self.S0 * self.delta)**(1/3) * (1 - self.delta + (3/4) * (self.delta - x**(4/3) / self.delta**(1/3))))**3 + 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) From afe93a5d6ba620baa830cb464be8142bd31506b3 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 13 Jul 2022 16:44:27 -0700 Subject: [PATCH 5/8] Generalise Shestakov problem to arbitrary p,q Input is now p and q, defaulting to 3 and -2 for previous case. Analytic solution updated to work for arbitrary p,q --- examples/shestakov_time_dependent.py | 46 +++++++++++++------ tango/extras/shestakov_nonlinear_diffusion.py | 19 +++++--- 2 files changed, 44 insertions(+), 21 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 86e56b7..953f4dc 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -61,7 +61,7 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise ): tlog.info("Initializing...") L, N, dx, x, nL, n = initialize_shestakov_problem() - test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx) + test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx, p=15, q=-14) fluxModel = NoisyFluxSpaceTime( FluxDoubleRelaxation( @@ -114,9 +114,18 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise num_iterations = 0 restarts = [] while num_iterations < maxIterations: - while solver.ok: - # Implicit time advance: iterate to solve the nonlinear equation! - solver.take_timestep() + 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 n = solver.profiles[label] # finished solution @@ -152,11 +161,18 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise "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 - err_history = np.concatenate(err_history) - profiles = np.concatenate(profile_history) - fluxes = np.concatenate(flux_history) + 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): @@ -206,18 +222,18 @@ def collect_statistics(inputs, num_samples=20): # Settings # ============================================================================== -timescales = [1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] +timescales = [1.0] #1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] damping_multipliers = [] #[1.0, 2.0, 5.0, 10] -noise_multipliers = [1e-2, 1e-1, 1.0] +noise_multipliers = [1e-2]#, 1e-1, 1.0] tolerance = 1e-2 -alpha_profile = 1 -alpha_flux = 0.05 +alpha_profile = 0.0 +alpha_flux = 0.0 noise_amplitude = 0.01 # Amplitude of noise to add noise_scalelength = 10.0 # Spatial scale length (number of cells) -num_samples = 20 # Number of repeated solves, to gather statistics +num_samples = 1 # Number of repeated solves, to gather statistics # ============================================================================== # MAIN STARTS HERE @@ -319,13 +335,13 @@ def collect_statistics(inputs, num_samples=20): # Settings # ============================================================================== -turb_timescale = 10.0 +turb_timescale = 1.0 damping_multiplier = 1.0 noise_multiplier = 1.0 tolerance = 1e-2 -alpha_ps = np.logspace(-1, 0, num=10) -alpha_ds = np.logspace(-2.5, -1.5, num=10) +alpha_ps = np.logspace(-2, 0, num=10) +alpha_ds = np.logspace(-2, 0, num=10) noise_amplitude = 0.01 # Amplitude of noise to add noise_scalelength = 10.0 # Spatial scale length (number of cells) diff --git a/tango/extras/shestakov_nonlinear_diffusion.py b/tango/extras/shestakov_nonlinear_diffusion.py index f8c9f3f..521a586 100644 --- a/tango/extras/shestakov_nonlinear_diffusion.py +++ b/tango/extras/shestakov_nonlinear_diffusion.py @@ -61,16 +61,16 @@ def _dxCenteredDifference(self, u, dx): 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[n] = -(dn/dx)^p * n^q """ - def __init__(self, dx, S0 = 1, delta = 0.1): + 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=3, q=-2, field='n') + super(ShestakovTestProblem, self).__init__(dx, p=p, q=q, field='n') self.S0 = S0 self.delta = delta @@ -85,15 +85,22 @@ def GetSource(self, x): return S def steady_state_solution(self, x, nL): - """Return the exast steady state solution for the Shestakov test problem + """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: """ - nright = ( nL**(1/3) + 1/3 * (self.S0 * self.delta)**(1/3) *(1-x) )**3 - nleft = ( nL**(1/3) + 1/3 * (self.S0 * self.delta)**(1/3) * (1 - self.delta + (3/4) * (self.delta - x**(4/3) / self.delta**(1/3))))**3 + + 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 From 3db33ecf37987ee80d915dd97f96d79b156ee958 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 28 Jul 2022 15:52:27 -0700 Subject: [PATCH 6/8] FluxRelaxationOscillation and FluxAverage decorators Modify the Shestakov model flux to more closely model the characteristics of GENE turbulence simulations. --- tango/extras/fluxrelaxation.py | 88 +++++++++++++++++++++++++++++++++- tango/extras/noisyflux.py | 6 +-- 2 files changed, 89 insertions(+), 5 deletions(-) diff --git a/tango/extras/fluxrelaxation.py b/tango/extras/fluxrelaxation.py index 56b1878..873ed09 100644 --- a/tango/extras/fluxrelaxation.py +++ b/tango/extras/fluxrelaxation.py @@ -6,7 +6,7 @@ class FluxRelaxation(object): Simple relaxation on a fixed timescale """ - def __init__(self, fluxModel, timescale): + def __init__(self, timescale, fluxModel): """ Inputs: fluxModel fluxmodel to be decorated. @@ -46,7 +46,7 @@ class FluxDoubleRelaxation(object): one timescale for the drive, and one for the damping. """ - def __init__(self, fluxModel, turb_timescale, damp_timescale): + def __init__(self, turb_timescale, damp_timescale, fluxModel): """ Inputs: fluxModel @@ -94,3 +94,87 @@ def get_flux(self, profiles): newFluxes[key] = self.lastTurb[key] ** 2 / self.lastDrive[key] return newFluxes + +class FluxRelaxationOscillation(object): + """Decorator that adds time dependence to fluxes with oscillation + Relaxation on a fixed timescale, with oscillation of specified relative magnitude + """ + + def __init__(self, timescale, amplitude, 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 + + amplitude Relative oscillation amplitude e.g. 0.2 is 20% + """ + 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 + self.amplitude = amplitude # Oscillation relative amplitude + self.timescale = timescale + self.time = 0 # Keeps track of time for the oscillation phase + + 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.copy() # Damping based on flux without oscillation + + # Add a relative oscillation + for key in newFluxes: + newFluxes[key] *= (1. + self.amplitude * np.sin(3. * self.time / self.timescale)) + self.time += 1 + return newFluxes + +class FluxAverage(object): + """Decorator that averages the flux over a given number of iterations + """ + def __init__(self, nsteps, 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 + nsteps Number of steps + """ + assert hasattr(fluxModel, "get_flux") and callable( + getattr(fluxModel, "get_flux") + ) + assert nsteps > 0 + + self.fluxModel = fluxModel + self.nsteps = nsteps + + def get_flux(self, profiles): + # Call the flux model to get the new flux + fluxes = self.fluxModel.get_flux(profiles) + + # Sum fluxes over nsteps + for i in range(self.nsteps - 1): + next_fluxes = self.fluxModel.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] /= self.nsteps + + return fluxes diff --git a/tango/extras/noisyflux.py b/tango/extras/noisyflux.py index 9bf7814..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 @@ -110,7 +110,7 @@ class NoisyFluxSpaceTime(object): fluxModel. The correlation length in space is controlled by the AR parameters """ - def __init__(self, fluxModel, amplitude, tac_x, dx, tac_t, dt): + def __init__(self, amplitude, tac_x, dx, tac_t, dt, fluxModel): """ Inputs: fluxModel fluxmodel to be decorated From cf2e4490865ee3adf05fd4eaaff0f0bc2bf5f61f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 28 Jul 2022 15:55:33 -0700 Subject: [PATCH 7/8] Exploring Shestakov fluxes with different p and q values Modified the shestakov_time_dependent.py example --- examples/shestakov_time_dependent.py | 123 ++++++++++++++------------- 1 file changed, 66 insertions(+), 57 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 953f4dc..2b16ae2 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -13,7 +13,7 @@ import tango.tango_logging as tlog from tango.extras import shestakov_nonlinear_diffusion -from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation +from tango.extras.fluxrelaxation import FluxRelaxation, FluxDoubleRelaxation, FluxRelaxationOscillation, FluxAverage from tango.extras.noisyflux import NoisyFlux, NoisyFluxSpaceTime import tango @@ -47,34 +47,37 @@ def __call__(self, t, x, profiles, HCoeffsTurb): # ============================================================================== def solve_system(noise_timescale = 1.0, # AR time of random noise - turb_timescale = 1.0, # Flux relaxation time - damping_timescale = 1.0, # Flux damping timescale + 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 = False, + plot_convergence = True, check_tol = 0.1 ): tlog.info("Initializing...") L, N, dx, x, nL, n = initialize_shestakov_problem() - test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx, p=15, q=-14) - - fluxModel = NoisyFluxSpaceTime( - FluxDoubleRelaxation( - test_problem, - turb_timescale, - damping_timescale, - ), - noise_amplitude, - noise_scalelength, - 1.0, # dx - noise_timescale, # noise timescale - 1.0 # dt - ) + test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx, p=p, q=q) + + nsteps = 200 + 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))) label = "n" turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) @@ -222,15 +225,18 @@ def collect_statistics(inputs, num_samples=20): # Settings # ============================================================================== -timescales = [1.0] #1e-3, 0.1, 0.2, 0.5, 1, 2, 5, 10] -damping_multipliers = [] #[1.0, 2.0, 5.0, 10] -noise_multipliers = [1e-2]#, 1e-1, 1.0] +timescales = np.logspace(np.log10(0.01), np.log10(10), 7) +oscillation_amplitudes = []#, 1e-2, 1e-1] +noise_multipliers = [0.2] + +tolerance = 1e-1 +alpha_profile = 0.03 +alpha_flux = 0.1 -tolerance = 1e-2 -alpha_profile = 0.0 -alpha_flux = 0.0 +shestakov_p = 15 +shestakov_q = -10 -noise_amplitude = 0.01 # Amplitude of noise to add +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 @@ -244,17 +250,19 @@ def collect_statistics(inputs, num_samples=20): fig2, ax2 = plt.subplots() fig3, ax3 = plt.subplots() -for damping_multiplier in damping_multipliers: +for oscillation_amplitude in oscillation_amplitudes: for noise_multiplier in noise_multipliers: mean_iterations = [] max_iterations = [] min_iterations = [] - for turb_timescale in timescales: - result = collect_statistics({"noise_timescale": turb_timescale * noise_multiplier, # AR time of random noise - "turb_timescale": turb_timescale, # Flux relaxation time - "damping_timescale": turb_timescale * damping_multiplier, # Flux damping timescale + 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) + "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}, @@ -264,7 +272,7 @@ def collect_statistics(inputs, num_samples=20): num_samples = num_samples) # Plot error history on figure 1 (for the last sample) - ax1.plot(result["err_history"], label=str(turb_timescale)) + 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"]) @@ -277,22 +285,19 @@ def collect_statistics(inputs, num_samples=20): 1./timescales, mean_iterations, "-o", - label=r"$\tau_{{damp}} / \tau_{{turb}} = {}, \tau_{{noise}} / \tau_{{turb}} = {}$".format(damping_multiplier, noise_multiplier), + 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) - #ax2.axvline(1.0 / damping_multiplier, linestyle="--", color="k") - #ax2.axvline(1.0 / noise_multiplier, linestyle=":", color='k') - # Plot total turbulence run time on figure 3 ax3.plot( 1. / timescales, np.array(mean_iterations) / timescales, "-o", - label=r"$\tau_{{damp}} / \tau_{{turb}} = {}, \tau_{{noise}} / \tau_{{turb}} = {}$".format(damping_multiplier, noise_multiplier), + label=r"$\tau_{{noise}} / \tau_{{flux}} = {}, A_{{flux}} = {}$".format(noise_multiplier, oscillation_amplitude), ) if num_samples > 1: @@ -312,7 +317,7 @@ def collect_statistics(inputs, num_samples=20): # Plot of iteration count against timescale ax2.set_xscale("log") ax2.set_yscale("log") -ax2.set_xlabel(r"Simulation run time / $\tau_{{turb}}$") +ax2.set_xlabel(r"Simulation run time / $\tau_{{flux}}$") ax2.set_ylabel("Iterations required") ax2.legend() @@ -321,7 +326,7 @@ def collect_statistics(inputs, num_samples=20): ax3.set_xscale("log") ax3.set_yscale("log") -ax3.set_xlabel(r"Simulation run time / $\tau_{{turb}}$") +ax3.set_xlabel(r"Simulation run time / $\tau_{{flux}}$") ax3.set_ylabel("Total run time required") ax3.legend() @@ -330,20 +335,22 @@ def collect_statistics(inputs, num_samples=20): plt.show() - # ============================================================================== # Settings # ============================================================================== -turb_timescale = 1.0 -damping_multiplier = 1.0 -noise_multiplier = 1.0 +shestakov_p = 9 +shestakov_q = -6 + +flux_timescale = 0.2 +noise_multiplier = 0.2 +oscillation_amplitude = 0.1 -tolerance = 1e-2 +tolerance = 1e-1 alpha_ps = np.logspace(-2, 0, num=10) alpha_ds = np.logspace(-2, 0, num=10) -noise_amplitude = 0.01 # Amplitude of noise to add +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 @@ -357,18 +364,20 @@ def collect_statistics(inputs, num_samples=20): for i, alpha_p in enumerate(alpha_ps): for j, alpha_d in enumerate(alpha_ds): - result = collect_statistics({"noise_timescale": turb_timescale * noise_multiplier, # AR time of random noise - "turb_timescale": turb_timescale, # Flux relaxation time - "damping_timescale": turb_timescale * damping_multiplier, # Flux damping timescale - "noise_amplitude": noise_amplitude, # Amplitude of AR(1) noise - "noise_scalelength": noise_scalelength, # Spatial scale length (cells) - "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) + 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) From b2abc903efef2d71cb1a372f824047ac385de3be Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 29 Jul 2022 17:16:59 -0700 Subject: [PATCH 8/8] Restructuring flux model decorators Rather than being classes which implement `get_flux` methods, now functions which replace the `get_flux` method of given object. The motivation is to combine the flux calculation method with the test problem initial and boundary conditions, and exact solution calculation. The goal is to make testing of different iterations and boundary condition calculations simpler. --- examples/shestakov_time_dependent.py | 139 +++++++++++------- tango/extras/fluxrelaxation.py | 123 +++++++++------- tango/extras/shestakov_nonlinear_diffusion.py | 3 + 3 files changed, 158 insertions(+), 107 deletions(-) diff --git a/examples/shestakov_time_dependent.py b/examples/shestakov_time_dependent.py index 2b16ae2..b47e5ce 100644 --- a/examples/shestakov_time_dependent.py +++ b/examples/shestakov_time_dependent.py @@ -12,35 +12,48 @@ import matplotlib.pyplot as plt import tango.tango_logging as tlog -from tango.extras import shestakov_nonlinear_diffusion +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 -def initialize_shestakov_problem(): - # Problem Setup - L = 1 # size of domain - N = 500 # number of spatial grid points - dx = L / (N - 1) # spatial grid size - x = np.arange(N) * dx # location corresponding to grid points j=0, ..., N-1 - nL = 1e-2 # right boundary condition - nInitialCondition = 1 - 0.5 * x - return (L, N, dx, x, nL, nInitialCondition) - -class ComputeAllH(object): - def __init__(self, test_problem): - self.test_problem = test_problem - - def __call__(self, t, x, profiles, HCoeffsTurb): - # n = profiles['default'] - # Define the contributions to the H coefficients for the Shestakov Problem +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.test_problem.H7contrib_Source(x) - - HCoeffs = tango.multifield.HCoefficients(H1=H1, H7=H7) - HCoeffs = HCoeffs + HCoeffsTurb - - return HCoeffs + H7 = self.H7contrib_Source(x) + return HCoeffsTurb + tango.multifield.HCoefficients(H1=H1, H7=H7) # ============================================================================== # Solve function @@ -62,25 +75,31 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise check_tol = 0.1 ): tlog.info("Initializing...") - L, N, dx, x, nL, n = initialize_shestakov_problem() - test_problem = shestakov_nonlinear_diffusion.ShestakovTestProblem(dx, p=p, q=q) - - nsteps = 200 - 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))) - - label = "n" - turbHandler = tango.lodestro_method.TurbulenceHandler(dx, x, fluxModel) + + 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, @@ -88,11 +107,11 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise thetaParams, ) field0 = tango.multifield.Field( - label=label, - rightBC=nL, - profile_mminus1=n, - compute_all_H=ComputeAllH(test_problem), - lodestroMethod=lodestroMethod, + 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) @@ -104,7 +123,13 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise tArray = np.array([0, 1e4]) # specify the timesteps to be used. solver = tango.solver.Solver( - L, x, tArray, maxIterations, tol, compute_all_H_all_fields, fields, + test_problem.L, + test_problem.x, + tArray, + maxIterations, + tol, + compute_all_H_all_fields, + fields, saveFluxesInMemory = plot_convergence ) @@ -130,10 +155,12 @@ def solve_system(noise_timescale = 1.0, # AR time of random noise num_iterations = maxIterations break + label = test_problem.label() + n = solver.profiles[label] # finished solution - source = test_problem.GetSource(x) - nss = test_problem.steady_state_solution(x, nL) + 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)) @@ -226,15 +253,15 @@ def collect_statistics(inputs, num_samples=20): # ============================================================================== timescales = np.logspace(np.log10(0.01), np.log10(10), 7) -oscillation_amplitudes = []#, 1e-2, 1e-1] +oscillation_amplitudes = [0.0]#, 1e-2, 1e-1] noise_multipliers = [0.2] tolerance = 1e-1 -alpha_profile = 0.03 +alpha_profile = 1 alpha_flux = 0.1 -shestakov_p = 15 -shestakov_q = -10 +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) @@ -335,6 +362,8 @@ def collect_statistics(inputs, num_samples=20): plt.show() +sys.exit(0) + # ============================================================================== # Settings # ============================================================================== diff --git a/tango/extras/fluxrelaxation.py b/tango/extras/fluxrelaxation.py index 873ed09..b519333 100644 --- a/tango/extras/fluxrelaxation.py +++ b/tango/extras/fluxrelaxation.py @@ -95,86 +95,105 @@ def get_flux(self, profiles): return newFluxes -class FluxRelaxationOscillation(object): +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 - """ - def __init__(self, timescale, amplitude, fluxModel): - """ - Inputs: + # 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 - timescale Ratio of flux relaxation timescale to coupling period - amplitude Relative oscillation amplitude e.g. 0.2 is 20% - """ - assert hasattr(fluxModel, "get_flux") and callable( - getattr(fluxModel, "get_flux") - ) - assert timescale > 0.0 + Note: The fluxModel object will be modified - 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 - self.amplitude = amplitude # Oscillation relative amplitude - self.timescale = timescale - self.time = 0 # Keeps track of time for the oscillation phase + # Returns - 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 + 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 newFluxes: - newFluxes[key] = ( - self.weight * self.lastFluxes[key] - + (1.0 - self.weight) * newFluxes[key] + for key in new_fluxes: + new_fluxes[key] = ( + weight * last_fluxes[key] + + (1.0 - weight) * new_fluxes[key] ) - self.lastFluxes = newFluxes.copy() # Damping based on flux without oscillation + last_fluxes = new_fluxes.copy() # Damping based on flux without oscillation # Add a relative oscillation - for key in newFluxes: - newFluxes[key] *= (1. + self.amplitude * np.sin(3. * self.time / self.timescale)) - self.time += 1 - return newFluxes + for key in new_fluxes: + new_fluxes[key] *= (1. + amplitude * np.sin(3. * time / timescale)) + time += 1 + return new_fluxes -class FluxAverage(object): + # 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 - """ - def __init__(self, nsteps, fluxModel): - """ - Inputs: + + # 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 - self.fluxModel = fluxModel - self.nsteps = nsteps - - def get_flux(self, profiles): + """ + 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 = self.fluxModel.get_flux(profiles) + fluxes = inner_get_flux(profiles) # Sum fluxes over nsteps - for i in range(self.nsteps - 1): - next_fluxes = self.fluxModel.get_flux(profiles) + 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] /= self.nsteps - + fluxes[key] /= nsteps return fluxes + + # Replace the get_flux function + fluxModel.get_flux = get_flux + return fluxModel diff --git a/tango/extras/shestakov_nonlinear_diffusion.py b/tango/extras/shestakov_nonlinear_diffusion.py index 521a586..508e235 100644 --- a/tango/extras/shestakov_nonlinear_diffusion.py +++ b/tango/extras/shestakov_nonlinear_diffusion.py @@ -41,6 +41,9 @@ def calc_flux_gradient(self, n, dx): flux = self.get_flux(n, dx) return self._dxCenteredDifference(flux, 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.