diff --git a/tests/+problemtests/validateallencahn.m b/tests/+problemtests/validateallencahn.m new file mode 100644 index 00000000..eeed6c29 --- /dev/null +++ b/tests/+problemtests/validateallencahn.m @@ -0,0 +1,72 @@ +function validateqg + +fprintf(' Testing Allen-Cahn Equations\n'); + +forcings = {8, @(t, ~, ~) cos(t)}; + +for fi = 1:numel(forcings) + forcing = forcings{fi}; + + model = otp.allencahn.presets.Canonical('Size', 16, 'Forcing', forcing); + model.Parameters.LinearizationPoint = model.Y0; + + %% Linear + [~] = otp.utils.Solver.Nonstiff(model.RHSLinear.F, model.TimeSpan, model.Y0); + fprintf(' Alternate RHSLinear RHS passed\n'); + + %% Non-Linear + [~] = otp.utils.Solver.Nonstiff(model.RHSLinear.F, model.TimeSpan, model.Y0); + fprintf(' Alternate RHSNonlinear RHS passed\n'); + + %% Reaction + [~] = otp.utils.Solver.Nonstiff(model.RHSReaction.F, model.TimeSpan, model.Y0); + fprintf(' Alternate RHSReaction RHS passed\n'); + + tc = model.TimeSpan(1); + y0 = model.Y0; + + f = @(t, y) model.RHSReaction.F(t, y); + japprox = model.RHSReaction.Jacobian(tc, y0); + jtrue = otp.utils.derivatives.jacobian(f, tc, y0); + + normj = norm(jtrue); + + if normj < eps + err = norm(jtrue - japprox); + else + err = norm(jtrue - japprox)/normj; + end + + tol = 1e-6; + assert(err < tol); + + fprintf(' Reaction Jacobian passed\n'); + + %% Diffusion + [~] = otp.utils.Solver.Nonstiff(model.RHSDiffusion.F, model.TimeSpan, model.Y0); + fprintf(' Alternate RHSDiffusion RHS passed\n'); + + tc = model.TimeSpan(1); + y0 = model.Y0; + + f = @(t, y) model.RHSDiffusion.F(t, y); + japprox = model.RHSDiffusion.Jacobian(tc, y0); + jtrue = otp.utils.derivatives.jacobian(f, tc, y0); + + normj = norm(jtrue); + + if normj < eps + err = norm(jtrue - japprox); + else + err = norm(jtrue - japprox)/normj; + end + + tol = 1e-6; + assert(err < tol); + + fprintf(' Diffusion Jacobian passed\n'); + +end + +end + diff --git a/tests/runalltests.m b/tests/runalltests.m index b309e0e3..b7f4a5f4 100644 --- a/tests/runalltests.m +++ b/tests/runalltests.m @@ -24,5 +24,6 @@ function runalltests(varargin) %% Problem specific tests problemtests.validateqg; +problemtests.validateallencahn; end diff --git a/toolbox/+otp/+allencahn/+presets/Canonical.m b/toolbox/+otp/+allencahn/+presets/Canonical.m index 8592f0ab..8550c0aa 100644 --- a/toolbox/+otp/+allencahn/+presets/Canonical.m +++ b/toolbox/+otp/+allencahn/+presets/Canonical.m @@ -3,7 +3,12 @@ methods function obj = Canonical(varargin) - params = otp.allencahn.AllenCahnParameters('Size', 64, 'Alpha', 0.1, 'Beta', 1, 'Forcing', 0, varargin{:}); + params = otp.allencahn.AllenCahnParameters('Size', 64, ... + 'Alpha', 0.1, ... + 'Beta', 1, ... + 'Forcing', 0, ... + 'LinearizationPoint', 0, ... + varargin{:}); x = linspace(0, 1, params.Size); [xs, ys] = meshgrid(x, x); diff --git a/toolbox/+otp/+allencahn/AllenCahnParameters.m b/toolbox/+otp/+allencahn/AllenCahnParameters.m index a8135659..421f345a 100644 --- a/toolbox/+otp/+allencahn/AllenCahnParameters.m +++ b/toolbox/+otp/+allencahn/AllenCahnParameters.m @@ -13,6 +13,9 @@ %Forcing is a forcing function or constant Forcing %MATLAB ONLY: {mustBeA(Forcing, {'numeric', 'function_handle'})} + + %The linearization point for the linear-non-linear splitting + LinearizationPoint %MATLAB ONLY: {otp.utils.validation.mustBeNumerical} end methods diff --git a/toolbox/+otp/+allencahn/AllenCahnProblem.m b/toolbox/+otp/+allencahn/AllenCahnProblem.m index 31bf036a..ccebe77f 100644 --- a/toolbox/+otp/+allencahn/AllenCahnProblem.m +++ b/toolbox/+otp/+allencahn/AllenCahnProblem.m @@ -1,5 +1,12 @@ classdef AllenCahnProblem < otp.Problem %ALLENCAHNPROBLEM + + properties (SetAccess = private) + RHSLinear + RHSNonlinear + RHSReaction + RHSDiffusion + end methods function obj = AllenCahnProblem(timeSpan, y0, parameters) @@ -9,10 +16,11 @@ methods (Access = protected) function onSettingsChanged(obj) - n = obj.Parameters.Size; - alpha = obj.Parameters.Alpha; - beta = obj.Parameters.Beta; + n = obj.Parameters.Size; + alpha = obj.Parameters.Alpha; + beta = obj.Parameters.Beta; forcing = obj.Parameters.Forcing; + uL = obj.Parameters.LinearizationPoint; if obj.NumVars ~= n^2 warning('OTP:inconsistentNumVars', ... @@ -25,16 +33,44 @@ function onSettingsChanged(obj) if ~isa(forcing, 'function_handle') f = @(t, y) otp.allencahn.fConstForce(t, y, L, alpha, beta, forcing); + + ft = @(~) forcing; + + flinear = @(t, y) otp.allencahn.fLinear(t, y, L, alpha, beta, ft, uL); + fnonlinear = @(t, y) otp.allencahn.fNonlinear(t, y, L, alpha, beta, ft, uL); + + freaction = @(t, y) otp.allencahn.fReaction(t, y, L, alpha, beta, ft); + fdiffusion = @(t, y) otp.allencahn.fDiffusion(t, y, L, alpha, beta, ft); + jreaction = @(t, y) otp.allencahn.jacobianReaction(t, y, L, alpha, beta, ft); + jdiffusion = @(t, y) otp.allencahn.jacobianDiffusion(t, y, L, alpha, beta, ft); else [x, y] = meshgrid(linspace(0, 1, n), linspace(0, 1, n)); x = x(:); y = y(:); ft = @(t) forcing(t, x, y); f = @(t, y) otp.allencahn.f(t, y, L, alpha, beta, ft); + + flinear = @(t, y) otp.allencahn.fLinear(t, y, L, alpha, beta, ft, uL); + fnonlinear = @(t, y) otp.allencahn.fNonLinear(t, y, L, alpha, beta, ft, uL); + + freaction = @(t, y) otp.allencahn.fReaction(t, y, L, alpha, beta, ft); + fdiffusion = @(t, y) otp.allencahn.fDiffusion(t, y, L, alpha, beta, ft); + jreaction = @(t, y) otp.allencahn.jacobianReaction(t, y, L, alpha, beta, ft); + jdiffusion = @(t, y) otp.allencahn.jacobianDiffusion(t, y, L, alpha, beta, ft); end obj.RHS = otp.RHS(f, ... 'Jacobian', @(t, u) otp.allencahn.jacobian(t, u, L, alpha, beta, forcing)); + + obj.RHSLinear = otp.RHS(flinear); + + obj.RHSNonlinear = otp.RHS(fnonlinear); + + obj.RHSReaction = otp.RHS(freaction, ... + 'Jacobian', jreaction); + + obj.RHSDiffusion = otp.RHS(fdiffusion, ... + 'Jacobian', jdiffusion); end end diff --git a/toolbox/+otp/+allencahn/fDiffusion.m b/toolbox/+otp/+allencahn/fDiffusion.m new file mode 100644 index 00000000..edce9c95 --- /dev/null +++ b/toolbox/+otp/+allencahn/fDiffusion.m @@ -0,0 +1,9 @@ +% Process splitting: +% f = f_diffusion + f_reaction +% Jacobian = Jacobian_diffusion + Jacobian_reaction +% +function du = fDiffusion(~, u, L, alpha, ~, ~) + +du = alpha*L*u; + +end diff --git a/toolbox/+otp/+allencahn/fLinear.m b/toolbox/+otp/+allencahn/fLinear.m new file mode 100644 index 00000000..326c9a89 --- /dev/null +++ b/toolbox/+otp/+allencahn/fLinear.m @@ -0,0 +1,11 @@ +% Linear-nonlinear splitting: f = f_linear + f_nonlinear +% f_linear(u) = Jac(uL)*u, f_nonlinear(u) = f(u)-Jac(uL)*u +% uL = linearization point, typically the solution at the beginning of the +% time step + + +function du = fLinear(~, u, L, alpha, beta, ~, uL) + +du = alpha*L*u + beta*(u - 3*(uL.^2).*u); + +end diff --git a/toolbox/+otp/+allencahn/fNonlinear.m b/toolbox/+otp/+allencahn/fNonlinear.m new file mode 100644 index 00000000..ec5adbdd --- /dev/null +++ b/toolbox/+otp/+allencahn/fNonlinear.m @@ -0,0 +1,10 @@ +% Linear-nonlinear splitting: f = f_linear + f_nonlinear +% f_linear(u) = Jac*u, f_nonlinear(u) = f(u)-Jac*u +% uL = linearization point, typically the solution at the beginning of the +% time step + +function du = fNonlinear(t, u, ~, ~, beta, forcing, uL) + +du = beta*( - u.^3 + 3*(uL.^2).*u ) + forcing(t); + +end diff --git a/toolbox/+otp/+allencahn/fReaction.m b/toolbox/+otp/+allencahn/fReaction.m new file mode 100644 index 00000000..bbd18fe2 --- /dev/null +++ b/toolbox/+otp/+allencahn/fReaction.m @@ -0,0 +1,9 @@ +% Process splitting: +% f = f_diffusion + f_reaction +% Jacobian = Jacobian_diffusion + Jacobian_reaction +% +function du = fReaction(t, u, ~, ~, beta, forcing) + +du = beta*(u - u.^3) + forcing(t); + +end diff --git a/toolbox/+otp/+allencahn/jacobianDiffusion.m b/toolbox/+otp/+allencahn/jacobianDiffusion.m new file mode 100644 index 00000000..87de41ef --- /dev/null +++ b/toolbox/+otp/+allencahn/jacobianDiffusion.m @@ -0,0 +1,9 @@ +% Process splitting: +% f = f_diffusion + f_reaction +% Jacobian = Jacobian_diffusion + Jacobian_reaction +% +function j = jacobianDiffusion(~, ~, L, alpha, ~, ~) + +j = alpha*L; + +end diff --git a/toolbox/+otp/+allencahn/jacobianReaction.m b/toolbox/+otp/+allencahn/jacobianReaction.m new file mode 100644 index 00000000..3a54e6bd --- /dev/null +++ b/toolbox/+otp/+allencahn/jacobianReaction.m @@ -0,0 +1,9 @@ +% Process splitting: +% f = f_diffusion + f_reaction +% Jacobian = Jacobian_diffusion + Jacobian_reaction +% +function j = jacobianReaction(~, u, L, ~, beta, ~) + +j = spdiags(beta*(1 - 3*u.^2), 0, size(L, 1), size(L, 2)); + +end diff --git a/toolboxPackaging.prj b/toolboxPackaging.prj index c9ddefa4..b4dd6978 100644 --- a/toolboxPackaging.prj +++ b/toolboxPackaging.prj @@ -1,5 +1,5 @@ - + ODE Test Problems Steven Roberts, Andrey A. Popov, Arash Sarshar, Adrian Sandu @@ -14,7 +14,7 @@ 0970ed36-a788-484f-beca-105eecd2562a - + true @@ -85,5 +85,59 @@ + + /Users/sandu/GIT/GitHub/ODE-Test-Problems/ODE Test Problems.mltbx + + + + /Applications/MATLAB_R2024b.app + + + + + + + + + + true + + + + + true + + + + + true + + + + + true + + + + + true + + + + + true + true + false + false + false + false + false + false + 15.5 + true + false + maci64 + true + \ No newline at end of file