Задание начального условия для модели с тонкой сеткой с использованием решения, полученного с более грубой сеткой

Задайте начальные условия для модели с тонкой сеткой, используя решение грубой сетки из предыдущего анализа.

Создание модели PDE и включение геометрии встроенной функции squareg.

model = createpde;
geometryFromEdges(model,@squareg);

Задайте коэффициенты, примените граничные условия и задайте начальные условия.

specifyCoefficients(model,'m',0,'d',1,'c',5,'a',0,'f',0.1);
applyBoundaryCondition(model,'dirichlet','Edge',1,'u',1);
setInitialConditions(model,10);

Создайте сравнительно грубую сетку с максимальной длиной кромки целевого элемента 0,1.

generateMesh(model,'Hmax',0.1);

Решите модель за весь промежуток времени от 0 до 0,02 секунды.

tlist = linspace(0,2E-2,20);
Rtotal = solvepde(model,tlist);

Интерполировать решение в начале координат в течение всего периода времени.

singleSpanSol = Rtotal.interpolateSolution(0,0,1:numel(tlist)); 

Теперь решите модель для первой половины периода времени. Это решение будет использоваться в качестве начального условия при решении модели с более тонкой сеткой для второй половины промежутка времени.

tlist1 = linspace(0,1E-2,10);
R1 = solvepde(model,tlist1);

Создайте интерполятор для интерполяции начального условия.

x = model.Mesh.Nodes(1,:)';
y = model.Mesh.Nodes(2,:)';
interpolant = scatteredInterpolant(x,y,R1.NodalSolution(:,end));

Создайте более тонкую сетку, установив максимальную длину кромки целевого элемента 0,05.

generateMesh(model,'Hmax',0.05);

Используйте результаты модели грубой сетки в качестве начального условия для модели с более тонкой сеткой. Для определения icFcn см. раздел Функция начальных условий.

setInitialConditions(model,@(region) icFcn(region,interpolant));

Решите модель для второй половины периода времени.

tlist2 = linspace(1E-2,2E-2,10);
R2 = solvepde(model,tlist2);

Интерполяция решений в начале координат для первой и второй половин интервала времени.

multispanSol1 = R1.interpolateSolution(0,0,1:numel(tlist1));
multispanSol2 = R2.interpolateSolution(0,0,1:numel(tlist2));

Постройте график всех трех решений в начале координат.

figure
plot(tlist,singleSpanSol)
hold on
plot(tlist1, multispanSol1,'r*')
plot(tlist2, multispanSol2,'ko')
legend('Overall solution','Coarse mesh solution', 'Fine mesh solution')

function u0 = icFcn(region,interpolant)
u0 = interpolant(region.x',region.y');
end