Sandu split Allen-Cahn

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
+

tests/runalltests.m

@@ -24,5 +24,6 @@ function runalltests(varargin)
 
 %% Problem specific tests
 problemtests.validateqg;
+problemtests.validateallencahn;
 
 end

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);

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

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

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

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

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

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

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

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

toolboxPackaging.prj

@@ -1,5 +1,5 @@
 <deployment-project plugin="plugin.toolbox" plugin-version="1.0">
-  <configuration name="ODE Test Problems" target="target.toolbox" target-name="Package Toolbox">
+  <configuration file="/Users/sandu/GIT/GitHub/ODE-Test-Problems/toolboxPackaging.prj" location="/Users/sandu/GIT/GitHub/ODE-Test-Problems" name="toolboxPackaging" target="target.toolbox" target-name="Package Toolbox">
     <param.appname>ODE Test Problems</param.appname>
     <param.authnamewatermark>Steven Roberts, Andrey A. Popov, Arash Sarshar, Adrian Sandu</param.authnamewatermark>
     <param.email />
@@ -14,7 +14,7 @@
     <param.products.version />
     <param.platforms />
     <param.guid>0970ed36-a788-484f-beca-105eecd2562a</param.guid>
-    <param.exclude.filters></param.exclude.filters>
+    <param.exclude.filters />
     <param.exclude.pcodedmfiles>true</param.exclude.pcodedmfiles>
     <param.examples />
     <param.demosxml />
@@ -85,5 +85,59 @@
     <fileset.depfun.included />
     <fileset.depfun.excluded />
     <fileset.package />
+    <build-deliverables>
+      <file location="${PROJECT_ROOT}" name="ODE Test Problems.mltbx" optional="false">/Users/sandu/GIT/GitHub/ODE-Test-Problems/ODE Test Problems.mltbx</file>
+    </build-deliverables>
+    <workflow />
+    <matlab>
+      <root>/Applications/MATLAB_R2024b.app</root>
+      <toolboxes>
+        <toolbox name="matlabcoder" />
+        <toolbox name="embeddedcoder" />
+        <toolbox name="fixedpoint" />
+        <toolbox name="matlabhdlcoder" />
+        <toolbox name="neuralnetwork" />
+      </toolboxes>
+      <toolbox>
+        <matlabcoder>
+          <enabled>true</enabled>
+        </matlabcoder>
+      </toolbox>
+      <toolbox>
+        <embeddedcoder>
+          <enabled>true</enabled>
+        </embeddedcoder>
+      </toolbox>
+      <toolbox>
+        <fixedpoint>
+          <enabled>true</enabled>
+        </fixedpoint>
+      </toolbox>
+      <toolbox>
+        <matlabhdlcoder>
+          <enabled>true</enabled>
+        </matlabhdlcoder>
+      </toolbox>
+      <toolbox>
+        <neuralnetwork>
+          <enabled>true</enabled>
+        </neuralnetwork>
+      </toolbox>
+    </matlab>
+    <platform>
+      <unix>true</unix>
+      <mac>true</mac>
+      <windows>false</windows>
+      <win2k>false</win2k>
+      <winxp>false</winxp>
+      <vista>false</vista>
+      <linux>false</linux>
+      <solaris>false</solaris>
+      <osver>15.5</osver>
+      <os32>true</os32>
+      <os64>false</os64>
+      <arch>maci64</arch>
+      <matlab>true</matlab>
+    </platform>
   </configuration>
 </deployment-project>