Code covered by the BSD License  

Highlights from
meshfree exponential integrator

meshfree exponential integrator

by

 

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