clear uIntDecompMatrixEvo uNLDecompMatrix timeGrid thetaMeshGridPlot thetaGridAnimate clear rGridAnimate rMeshGridPlot pEigenBasis pNLMatrix pNLZeroDecomp pIntMatrixEvo clear pNLDecompMatrix pNLMatrixSliceTemp vEigenBasis vEigenFncZero vNLDecompMatrix vNLZeroDecompnb clear pEigenFncZer o entropyConstXGrid JScaledValues eigenFreqMatrix timePeriod pEigenFncPertTemp clear vNLMatrix rGrid thetaGrid pConstIntMatrix vBar pIntMatrixIteration clear smallDivisorsBasicNoZero zeroCorrectionBasic nonResCorrectionBasic clear pEigenFncIterationLoop1 pEigenFncIterationLoop2 clear pOrigNLMatrix pOrigNLDecompMatrix pOrigNLZeroDecomp uOrigNLDecompMatrix clear pEigenFncIterationLoop pEigenFncPertTemp2 clear pIntMatrixEvoCorrection % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %Choose the parameters that will determine the number of eigenfunctions % that are used. The set of eigenfunctions is a two index set, these % indicies being "rBasisSize" and "thetaBasisSize" %The eigenfunctions are indexed as % phi_(n,m,1) = cos(n*theta)*J_n(lambda_(n,m)*sigmaBar*r) % phi_(n,m,2) = sin(n*theta)*J_n(lambda_(n,m)*sigmaBar*r) % where thetaBasis = n for n = 0,1,...,N-1 % and rBasisSize = m for m = 1,...,M %Each pair (n,m) corresponds to two different eigenfunctions % (the cos and sin variant) This is excepting the case of n = 0, m = 0, % which corresponds to the constant eigenfunction. % so there will be a total of (N-1)*M + 1 eigenfunctions in the p-basis and % effectively (N-1)*M eigenfunctions in the u-basis (because we will always % have that phi_(0,0,1) = c, we know psi_(0,0,1) = grad(phi_(0,0,1)) = 0) %However, really we lose another M eigenfunctions, as for n=0 the sin % variant phi_(0,m,2) = 0 for all m, but these still have a row in % matricies for simplicity % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %Set whether to run the initialization part of the code that: % 1. sets all the constants (background pressure, entropy, number of grid % points, which eigenfunction we are choosing, number of eigenfunctions, etc.) % 2. Calculates the basis of eigenfunctions for those constants perform_Initialization = true; %Set whether to run an iteration method or to just evolve the starting % guess nonlinearly and return the result perform_Basic_Iteration = true; %Set whether or not to improve the initial guess of pBar + alpha*phi_k % before iteration starts, (this means the O(alpha^2) terms that come % from (1/2)*alpha^2*D^2f(pBar)[phi_k,phi_k]) and if so whether to still % evolve the uncorrected initial guess anyway performInitialGuessCorrection = false; performUncorrectedEvoAsWell = true; number_Iterations = 10; % NOTE: to use the initialization correction function, one must ensure that % if we are perturbing the eigenfunction (k1,k2,1) or (k1,k2,2), then % 2*k1 is within our calculated eigen basis, i.e. that 2*k1 < N %Set whether the main iteration scripts will prompt the user to choose the % theta-index and r-index of the fixed k-mode prompt_User_K_Mode_Indicies = true; %Set whether the main iteration scripts will prompt the user to choose the %size of the eigenfunction basis (thetaBasisSize and rBasisSize) prompt_User_Eigen_Basis_Sizes = true; % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if perform_Initialization == true %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %prompt user to choose between an isentropic disk (constant entropy) % or a disk with a radially piecewise constant entropy profile %This is done by prompting the user to specify the number of entropy levels % input of 0 = isentropic case % input of n = radial profile s(r) with n-1 jumps between n levels prompt = ['Specify the number of distinct levels for a radial piecewise constant entropy profile s(r) (1 corresponds to the isentropic case): ']; radial_Entropy_Levels = input([prompt,newline]) % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % if radial_Entropy_Levels == 1 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % run("TwoD_Disk_TimeEvo_Isentropic_Main_Iteration_v2.m") % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % end % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % if radial_Entropy_Levels > 1 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % run("TwoD_Disk_TimeEvo_NonIsentropic_Iteration_Main.m") % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % end % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
twoD_Disc_TimeEvo_Main.m