Code covered by the BSD License

# meshfree exponential integrator

### Stefan Rainer (view profile)

experimetal code for meshfree exponential integrators

[interpolant,points]=main
```function [interpolant,points]=main
%====================================================================
%            Meshfree exponential Integrator 2D
%   Equation: (see function timeEvolution)
%     u_dt = F(u) = epsilon*(u_xx)-alpha*(u_x)+rho*u*(u-1/2)*(1-u)
%   Reference:
%        Meshfree exponential integrators
%        M. Caliari, A. Ostermann and S. Rainer
%        SIAM J. Sci. Comput (2013)
%===================================================================

rbf = @(r,e) max(1-e.*r,0).^6.*(35*(e.*r).^2+18*e.*r+3);
rbf_x = @(r,e,dx) -56*dx.*e.^2.*(5*e.*r+1).*max(1-e.*r,0).^5;
rbf_xx = @(r,e) 56*e.^2.*((35*(e.*r).^2-4*e.*r-1).*max(1-e.*r,0).^4 ...
-(5*e.*r+1).*max(1-e.*r,0).^5);

isInDomain = @(u) (abs(u(:,1))<=2 & abs(u(:,2))<=2);
minDistance = 1e-3;

% shape parameter of RBF function
shapeParameter = 2;

% tolerance in space per time step
tolerance = 1e-3;

% final integration time
finalTime = 0.1;

% set of candidate intepolation points for first residual subsampling
[p1,p2] = meshgrid(linspace(-2,1,10));
candidatePoints = [p1(:),p2(:)];

% initial value
initfun = @(x) exp(-10*((x(:,1)+0.4).^2+(x(:,2)+0.4).^2));

[interpolant,points] = calculateRbfInterpolant(candidatePoints,initfun,tolerance);

[xx,yy] = meshgrid(linspace(-2,2,100));
surf(xx,yy,reshape(interpolant([xx(:),yy(:)]),size(xx)));
drawnow;
pause;

dt=1e-2;
t = 0;
while t<finalTime
t = t+dt;
funTimeEvolution = @(points,tol,eigenValues) timeEvolution(points,dt,...
interpolant,tol,eigenValues);
integrationPoints = [points];

%   integrationPoints = sort(unique_tol(integrationPoints,1e-3));

[newInterpolant,newpoints,newDt]=calculateRbfInterpolant(integrationPoints,...
@(p,eigenValues) funTimeEvolution(p,0.1*tolerance,eigenValues),...
tolerance,dt);
if newDt>0
% The advection constant for the problem is 6
[interpolant,points] = calculateRbfInterpolant(points+6*dt,...
newInterpolant,tolerance);
disp(sprintf('Number of points: %d ',size(points,1)));
surf(xx,yy,reshape(interpolant([xx(:),yy(:)]),size(xx)));
hold on;
plot3(points(:,1),points(:,2),1.1*ones(size(points,1)),'k.');
title(sprintf('Number of points: %d ',size(points,1)));
axis([-2,2,-2,2,0,1.1]);
hold off;
drawnow;
pause(0.001);
else
disp('Timestep rejected');
t = t-dt;
end
dt = abs(newDt);
if t+dt > finalTime
dt = finalTime-t;
end
dt
t
end

%=======================================================================
%        Residual Subsampling procedure
%=======================================================================
function [interpolant,points,dt] = calculateRbfInterpolant(candidatePoints,...
fun,tolerance,varargin)

numberCandidatePoints = length(candidatePoints);
pointsCoarsable = [];

% check for error at check points
[checkPoints,index] = calculateCheckPoints(candidatePoints);
allPoints = [candidatePoints;checkPoints];

if nargout<=2
rhs = fun(allPoints);
else
% time evolution
[rhs,dt,eigenValues] = fun(allPoints,[]);
if dt<0
interpolant =  [];
points = [];
return;
end
end

checkPointsValues1 = rhs(numberCandidatePoints+1:end);
rhs = rhs(1:numberCandidatePoints);

% calculate initial interpolant
DM = distanceMatrix(candidatePoints,candidatePoints);
IM = sparse(rbf(DM,shapeParameter));
coef = IM\rhs;
interpolant = @(dsites) evaluateInterpolant(dsites,candidatePoints,coef);

% evaluate at check points
checkPointsValues2 = interpolant(checkPoints);

% find points to be refined
errorCheckPoints = abs(checkPointsValues2-checkPointsValues1);
pointsToRefine = checkPoints(errorCheckPoints>tolerance,:);

% if no points are refined find coarsable points
if (isempty(pointsToRefine) & nargin<=3 & numberCandidatePoints>2)
testCoarsable = zeros(size(candidatePoints,1),1);
for i=1:size(candidatePoints,1)
testCoarsable(i) = max(abs(errorCheckPoints(index{i})));
end
pointsCoarsable = candidatePoints(testCoarsable<0.1*tolerance,:);
end
if (numberCandidatePoints<=2)
pointsCoarsable =[];
end

while (~isempty(pointsToRefine) | ~isempty(pointsCoarsable))
if ~isempty(pointsToRefine)
if ~isempty(pointsCoarsable)
interpolant = prevInterpolant;
candidatePoints = prevCandidatePoints;
break;
end
candidatePoints = [candidatePoints;pointsToRefine];
else
prevCandidatePoints = candidatePoints;
prevInterpolant = interpolant;

sortedPointsToCoarse = sortrows([candidatePoints,testCoarsable],3);
candidatePoints = setdiff(candidatePoints,sortedPointsToCoarse(1,1:2),'rows');
end

% check for error at check points
[checkPoints,index] = calculateCheckPoints(candidatePoints);
allPoints = [candidatePoints;checkPoints];
numberCandidatePoints = length(candidatePoints);

if nargout<=2
rhs = fun(allPoints);
else
[rhs,dt,eigenValues] = fun(allPoints,eigenValues);
end

checkPointsValues1 = rhs(numberCandidatePoints+1:end);
rhs = rhs(1:numberCandidatePoints);

% calculate inital interpolant
DM = distanceMatrix(candidatePoints,candidatePoints);
IM = rbf(DM,shapeParameter);
coef = IM\rhs;
interpolant = @(dsites) evaluateInterpolant(dsites,candidatePoints,coef);

% evaluate at check points
checkPointsValues2 = interpolant(checkPoints);

% find points to be refined
errorCheckPoints = abs(checkPointsValues2-checkPointsValues1);
pointsToRefine = checkPoints(errorCheckPoints>tolerance,:);

% if no points are refined find coarsable points
if (isempty(pointsToRefine) & nargin<=3 & numberCandidatePoints>2)
testCoarsable = zeros(size(candidatePoints,1),1);
for i=1:size(candidatePoints,1)
testCoarsable(i) = max(abs(errorCheckPoints(index{i})));
end
pointsCoarsable = candidatePoints(testCoarsable<0.1*tolerance,:);
end
if (numberCandidatePoints<=2)
pointsCoarsable =[];
end
end
points = candidatePoints;
end

%=======================================================================
%   Definition of the PDE
%   time march scheme (exponential Rosenbrock method order 3 (exprb32)
%=======================================================================

function [u,dt,eigenValues] = timeEvolution(points,dt,interpolant,tol,eigenValues)
% Equation: u_dt = F(u) = epsilon*(u_xx)-alpha*(u_x)+rho*u*(u-1/2)*(1-u)

epsilon = 5e-2;
alpha = 6;
rho = 70;
%     warning('you are solving a linear problem with
DM = distanceMatrix(points,points);
IM = sparse(rbf(DM,shapeParameter));

[R,dummy,p] = chol(IM,'vector');

if dummy~=0
error('singular IM matrix: some points probably coincide');
end

un = interpolant(points);

% Linearize at un
EM = diffusion(DM,points);
diffusionPart = @(u) evalOperator(u,R,p,EM);
JnSymm = @(u) (epsilon*diffusionPart(u)+rho*(3*un.*(1-un)-1/2).*u);
EM = advection(DM,points);
advectionPart = @(u) evalOperator(u,R,p,EM);
JnSSymm = @(u) -alpha*advectionPart(u);

F = @(u) (epsilon*diffusionPart(u)-alpha*advectionPart(u)+...
rho*u.*(u-1/2).*(1-u));

g = @(u) F(u)-JnSymm(u)-JnSSymm(u);

%time step
Fn = F(un);
[phi1,errest,info,eigenValues] = phileja(dt,JnSymm,JnSSymm,1,dt*Fn);

Un2 = un+phi1;
Dn2 = g(Un2)-g(un);

[phi3,errest,info]=phileja(dt,JnSymm,JnSSymm,3,2*dt*Dn2,[],eigenValues);

err = phi3;

u = Un2+err;

normerror = norm(err,inf);

factor = min(2,max(0.6,abs(0.9*(tol/normerror)^(1/2))));
if rho~=0
dt = factor*dt;
end

if normerror> tol
fprintf('integration error too big %e',normerror);
dt = -dt;
end
end

function EM = advection(DM,points)
dx = differenceMatrix(points(:,1),points(:,1));
EM = rbf_x(DM,shapeParameter,dx);
end

function EM = diffusion(DM,points)
EM = rbf_xx(DM,shapeParameter);
end

function Au = evalOperator(u,R,p,EM)
coef = zeros(size(u));
coef(p) = R\(R'\(u(p)));
Au = EM*coef;
end

%========================================================

function [checkPoints,index] = calculateCheckPoints(points)
tri = delaunay(points(:,1),points(:,2));
checkPoints = circle(tri,points(:,1),points(:,2));
filter = isInDomain(checkPoints);
checkPoints = checkPoints(filter,:);
tri = tri(filter,:);
index = cell(size(points,1),1);
[checkPoints,pointer,idx] = uniqueTol(checkPoints,minDistance);
for i=1:size(points)
index{i}=find(tri(:,1)==i | tri(:,2) == i | tri(:,3) == i);
index{i} = unique(pointer(index{i}));
end
end

function values = evaluateInterpolant(dsites,ctrs,coef)
DM = distanceMatrix(ctrs,dsites);
EM = rbf(DM,shapeParameter);
values = EM*coef;
end

function [result,pointer,idx]=uniqueTol(v,tol)
% result: set with unique points
% pointer: vector that points from the old set to the new set
% idx: indices of the old point contained in the new set
DM = distanceMatrix(v,v);

idx = 1:size(v,1);
i=1;
pointer = zeros(length(v),1);
while ~isempty(idx(i+1:end))
todelete = find(DM(idx(i+1:end),idx(i))<=tol)+i;
pointer(idx(todelete))=i;
idx = idx(setdiff(1:length(idx),todelete));
i = i+1;
end
pointer(pointer==0) = 1:length(idx);
result = v(idx,:);
end

end```

Contact us