How to implement periodic boundary conditions for 2D PDE

Hi all,

I'm trying to solve the diffusion equation in a 2D space but I need to set the left and right boundaries to periodic. Unfortunately, I don't think matlab has this functionality built in. Does anyone know of any way that this might be possible?

Any help would be greatly appreciated!

Here's my example code:

N0 = 1e19;
D0 = 1.49e-6;
Q = 1.17; 
kBoltz = 8.617e-5;
T=1000+273;
Deff=D0*exp(-Q/(kBoltz*T));
endTime = 120*60;

model = createpde(1);
width = 30e-4;
height = 30e-4;
dopedwidth = 1e-4;

g=diff_geom(width,height,dopedwidth);
geometryFromEdges(model,g);
hmax = .00005;
msh = generateMesh(model,'Hmax',hmax);
d = 1;
%c = @(~,state) Deff*(state.u/N0).^2;
c = Deff;
m = 0;
a = 0;
f = 0;
specifyCoefficients(model,'m', m , 'd' ,d, 'c',c, 'a',a, 'f',f);
applyBoundaryCondition(model,'dirichlet','Edge',2,'u',N0);
setInitialConditions(model,0);
tlist = 0:50:endTime;
numNodes = size(p,2);
u0(1:numNodes) = 0;
setInitialConditions(model,N0,'edge',2);
model.SolverOptions.RelativeTolerance = 1.0e-3;
model.SolverOptions.AbsoluteTolerance = 1.0e-4;
R = solvepde(model,tlist);
u = R.NodalSolution;

figure;
pdeplot(model,'XYData',u(:,end),'Contour','off','ColorMap','jet');
title(sprintf('Diffusion In Sample, Transient Solution( %d seconds)\n', ...
  tlist(1,end)));
caxis([0 1e19]);
xlabel 'X-coordinate (cm)'
ylabel 'Y-coordinate (cm)'
axis ([1e-3 2e-3 0 1e-3]);

Answer 1

[Ravi Kumar]

Hi Robert,

PDE Toolbox does not have an interface to specify periodic BCs. However, it is easy to modify the system equations to enforce periodicity if your geometry is simple and your mesh has identical number of nodes on the periodic boundary pair.

I can attempt a workaround if you provide me the complete reproduction code. Currently, the function diff_geom is unknown, could you provide it.

Regards, Ravi

[Robert Smith]

Hi Ravi, thanks so much for your help. Here's the function, sorry for not including it:

function [geomat] = diff_geom(regxwidth,regywidth,dopedwidth)
seg1End = regxwidth/2 - dopedwidth/2;
seg3Start = regxwidth/2 + dopedwidth/2;
geomat = [ lineMat(0,seg1End,0,0) ; lineMat(seg1End,seg3Start,0,0) ; lineMat(seg3Start,regxwidth,0,0) ; lineMat(regxwidth,regxwidth,0,regywidth) ; lineMat(regxwidth,0,regywidth,regywidth) ; lineMat(0,0,regywidth,0)]';

[Ravi Kumar]

Hi Robert,

I still need 'lineMat' which looks like something you might be loading from a MAT-file before executing your script. In the absence of the correct geometry I tried solving your setup on a simple square, but that did not give any meaningful solution. It would help if you provide a complete reproduction code.

Regards Ravi

[Robert Smith]

Hi Ravi,

Thanks again for looking into this. I apologize again for forgetting this part of the code:

function [lineMat] = lineMat(xStart,xEnd,yStart,yEnd)
lineMat = [2 , xStart , xEnd , yStart , yEnd , 1 , 0 , 0 , 0 , 0 , 0 , 0];

[Ravi Kumar]

Thanks. Your original code works now.

I did not realize before that you have a transient problem. The workaround I have is suitable for steady-state setup. Is there a possibility that you can convert the problem (by changing the coefficients) such that you can solve an equivalent steady-state problem?

Main issue is since periodic BC is not supported as type of BC, you need to get the system matrices and alter them to mimic periodicity. This involves some coding, which I can provide you, but in transient analysis this process need to be repeated at ever time-step. So essentially implementing periodic BC is equivalent to solving the time-integration part on your own using the system matrices, which involves significant additional coding.

[Robert Smith]

Hmm, that's unfortunate. Do you think there is a way to use the nonconstatn boundary conditions syntax to force periodicity (documented here: specify boundary conditions and Solve PDEs with Nonconstant Boundary Conditions)? I was wondering if there was a way to set u (the solution) at the left boundary equal to the right by using the state.u or state.time functions.

I guess I could work around this by expanding my simulation to include 3 or more periods but I'm worried that this will drastically increase the simulation time...

[Christian Digel]

Hi Ravi,

I just found you're answer here. I would be interested in the steady state implementation of the pereodic boundary condition. Could you give a small example how to manipulate the matrices ? Thank you in advance!

Kind Regads

Chris