CODEWITHSUNDEEP
matlabtoolboxpde
E. Puz
E. Puz

Reputation: 1

Matlab PDE Toolbox Periodic Boundary Conditions

How can we modify the matlab PDE toolbox code to implement periodic boundary conditions?

Upvotes: 0

Views: 1910

Answers (1)

E. Puz
E. Puz

Reputation: 1

I was able to implement periodic boundary conditions in 2016a for use within the pde toolbox. It took quite a bit of time because so many of the functions and data structures are not documented. Please note that my domain is a three dimensional cube, but the approach here will work for 2D too. (I have tested it for both parabolic and hyperbolic problems). It seems that periodic boundary conditions should be possible by altering the PDE toolbox's "H" matrix, but I haven't been able to get such an approach to work for ellipic, parabolic, or hyperbolic problems.

Here is the approach for implementing periodic boundary conditions for parabolic and hyperbolic problems (using the method of lines of the pde toolbox):

  1. Create a mesh using a delaunay triangulation of your points (make sure nodes touching the boundary of the domain have nodes which "match" on the opposite face). Note that the nodes of meshes produced by the pde toolbox for a 3D cube cannot be used, since they do not have the appropriate periodic structure.

  2. Treat the periodic boundary condition as a time dependent dirichlet boundary condition. You may have a bit of trouble getting matlab to treate your boundary conditions at time dependent ones (I have to omit the details here). However, here is the function that I used in the applyBoundaryCondition() function is as follows:

function [ bcMatrix ] = tdependentdiri( region, state )

% nan is used to "trick" matlab into treating the boundary conditions as time-dependent dirichlet

if( isnan( state.time ) )

bcMatrix = NaN;

else

bcMatrix = state.u;

end

end

  1. Update the fas data (from genmeshinternal() within generateMesh()). Unfortunately, we can't see the code for genmeshinternal()). However fas is cell with six entries, each of which is a 2xN matrix. The second row of each entry. The first row designates a tetrahedra whose face is on the boundary of the domain. The second row specifies which node (of the tetrahedra) is located on the interior of the domain. Write an algorithm to update fas for your mesh. In fact, you need to update each of the variables (using your mesh data) on the following line (within generateMesh()):

[nodes, fas, tet, Hmax, Hmin] = genmeshinternal(self.Geometry,Hmax,Hmin,geomOrder);

Note that nodes and tet are taken care of by the delaunay triangulation.

  1. The "periodic mapping" of nodes and the boundary values (which are changing with time) should be done within callValueFuncOnFace(). Basically you need to change the value of the current solution on the boundary to the value at its periodic partner's location (here indicated by pbc_indx):

state.u = self.uN(:,pbc_indx); % self.uN holds the current time's solution

I put this right after the line which reads: appRegion = self.applicationRegion(xyzPts(:,i), faceNormals(:,i));

There could be an easier way, but at least this one eventually works.

Upvotes: 0

Related Questions