Large residual error on the initial value condition

[BOBBY-LUKA]

I tried to modify your code by making it to solve transient partial differential equation. The error on the imposed initial value is extremely high compared to the PDE and the boundary condition. The mapped solution on the circular disc appeared with some patches on it, I have attached one of the solutions here. How can I improve the performance of the PINN to give same solution with the one from FEA?

Find below the code and the Loss function

%create geometry
model = createpde;
r = [1,0,0,0.025]';
g=[r];
g = decsg(g);
geometryFromEdges(model,g);
pdegplot(model,EdgeLabels="on")
axis equal
% apply initial value condition
U= 72201;
setInitialConditions(model,U);
%Specify zero Dirichlet boundary conditions on all edges
applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u', 28.542);
% specify PDE coefficient
pdeCoeffs.m = 0;
pdeCoeffs.d = 0.003;
pdeCoeffs.c= 4.1295E-10;
pdeCoeffs.a = 3.8566E-5;
pdeCoeffs.f = 0.2785;
specifyCoefficients(model,'m',0,'d', 0.003,'c', 4.1295E-10,'a', 3.8566E-5,'f', 0.2785,'face',1);

% Generate the mesh
msh = generateMesh(model,Hmax=0.001,Hgrad=2, ...
Hedge={1:model.Geometry.NumEdges, 0.00002});
figure
pdemesh(model);
axis equal

%Generate Spatial Data for Training PINN
boundaryNodes = findNodes(msh,"region", ...
Edge=1:model.Geometry.NumEdges);
domainNodes = setdiff(1:size(msh.Nodes,2),boundaryNodes);
domainCollocationPoints = msh.Nodes(:,domainNodes)';

% Extract the number of nodes
numNodes = size(msh.Nodes, 2);

% Extract the number of interior nodes

numInteriorNodes = size(domainNodes, 2);

%generate time points for the number of nodes between 0 and 60
numSteps= numInteriorNodes;
startValue=0;
endValue=60;
stepSize=(endValue-startValue)/(numSteps-1);
dataT=startValue:stepSize:endValue;
dataT=transpose(dataT);
%concatenate
domainCollocationPoints =[domainCollocationPoints,dataT];
% Retrieve the node coordinates

nodes = model.Mesh.Nodes;

% Initialize the solution array with the initial condition

IC= 77244;
setInitialConditions(model, IC);
% Training the model requires a data set of collocation points that enforce the boundary conditions, enforce the initial conditions,
BCiNodes = findNodes(model.Mesh,"region","Edge",4);
% initial time for initial value and boundary conditions.
T0 = zeros(size(BCiNodes));

T0IC=zeros(size(domainNodes));

%Define Network Architecture
numLayers = 4;
numNeurons = 30;
%Create the neural network.
layers = featureInputLayer(3);

for i = 1:numLayers-1
layers = [
layers
fullyConnectedLayer(numNeurons)
geluLayer];
end

layers = [
layers
fullyConnectedLayer(1)];

pinn = dlnetwork(layers);
%Specify Training Options
numEpochs =25;
miniBatchSize = 250;
initialLearnRate = 0.001;
learnRateDecay = 0.0005;
%Convert the training data to deep learening array(Deep Learning Toolbox) objects.
ds = arrayDatastore(domainCollocationPoints);
mbq = minibatchqueue(ds,MiniBatchSize=miniBatchSize, ...
MiniBatchFormat="BC");
%Initialize the average gradients and squared average gradients.
averageGrad = [];
averageSqGrad = [];
%Calculate the total number of iterations for the training progress monitor.
numIterations = numEpochs* ...
ceil(size(domainCollocationPoints,1)/miniBatchSize);
%Initialize the (Deep Learning Toolbox) object.
monitor = trainingProgressMonitor(Metrics="Loss", ...
Info="Epoch", ...
XLabel="Iteration")
%Train PINN
iteration = 0;
epoch = 0;
learningRate = initialLearnRate;
while epoch < numEpochs && ~monitor.Stop
epoch = epoch + 1;
reset(mbq);
while hasdata(mbq) && ~monitor.Stop
iteration = iteration + 1;
XYT = next(mbq);
X = XYT(1,:);
Y = XYT(2,:);
T = XYT(3,:);
% Evaluate the model loss and gradients using dlfeval.
[loss,gradients] = dlfeval(@modelLoss,model,pinn,X,Y,T,T0,pdeCoeffs);
% Update the network parameters using the adamupdate function.
[pinn,averageGrad,averageSqGrad] = ...
adamupdate(pinn,gradients,averageGrad, ...
averageSqGrad,iteration,learningRate);
end
% Update learning rate.
learningRate = initialLearnRate / (1+learnRateDecay*iteration);
% Update the training progress monitor.
recordMetrics(monitor,iteration,Loss=loss);
updateInfo(monitor,Epoch=epoch + " of " + numEpochs);
monitor.Progress = 100 * iteration/numIterations;
end

% find and plot the solution predicted by the PINN for the number of nodes.
numSteps_1= numNodes;
startValue=0;
endValue=60;
stepSize=(endValue-startValue)/(numSteps_1-1);
dataT=startValue:stepSize:endValue;
T_1=dlarray(dataT,"CB");
nodesDLarray = dlarray(msh.Nodes,"CB");
BOBBY =[nodesDLarray;T_1];
Upinn = gather(extractdata(predict(pinn,BOBBY)));

figure;
pdegplot(model);
axis equal
hold on

pdeplot(model,XYData=Upinn, ColorMap= "jet")
cmap=colormap;
colormap(cmap);
set(colorbar,'Ticks',[]);
set(colorbar,'TickLabels',[]);
set(colorbar,'Visible','on');
hold off
xlabel('$x$',interpreter='latex')
ylabel('$y$',interpreter='latex')
zlabel('$u(x,y)$',interpreter='latex')
title(sprintf(['PINN Predicted Solution: ' ...
'L2 Error = %0.1e'],norm(Upinn-true)))

LOSS FUNCTION
function [loss,gradients] = modelLoss(model,net,X,Y,T,T0,pdeCoeffs)
dim=1;
XYT=cat(dim,X,Y,T);
[U,~] = forward(net,XYT);

% Compute gradients of U.
gradU = dlgradient(sum(U,"all"),{X,Y,T},EnableHigherDerivatives=true);
Ux = gradU{1};
Uy = gradU{2};
Ut = gradU{3};

% Calculate second-order derivatives with respect to X,Y,T.
for i=1:dim
U2x = dlgradient(sum(Ux,"all"),X,EnableHigherDerivatives=true);
U2y = dlgradient(sum(Uy,"all"),Y,EnableHigherDerivatives=true);
U2xy=U2x+U2y;
end

% Enforce PDE. Calculate lossF.
res = pdeCoeffs.d.*Ut-pdeCoeffs.c.*U2xy+pdeCoeffs.a.*U-pdeCoeffs.f;
zeroTarget = zeros(size(res),"like",res);
mseF = l2loss(res,zeroTarget);

% Enforce boundary conditions. Calculate lossU.
actualBC = []; % Boundary information
BC_XYT0 = []; % Boundary coordinates
% Loop over the boundary edges and find boundary
% coordinates and actual BC assigned to PDE model.
numBoundaries = model.Geometry.NumEdges;
for i=1:numBoundaries
BCi = findBoundaryConditions(model.BoundaryConditions,"Edge",i);
BCiNodes = findNodes(model.Mesh,"region","Edge",i);
BC_XY =model.Mesh.Nodes(:,BCiNodes);
BC_XYT0 = [BC_XY;T0]; % Concatenate T0 with the current boundary coordinates
actualBCi = ones(1,numel(BCiNodes))*BCi.u;
actualBC = actualBCi;
end
BC_XYT0 = dlarray(BC_XYT0,"CB"); % Format the coordinates
[predictedBC,~] = forward(net,BC_XYT0);
mseU = l2loss(predictedBC, actualBC);

%Enforce initial conditions. Calculate lossIC.
actualIC = []; % Boundary information
IC_XYT0IC = []; % Boundary coordinates
% Loop over the boundary edges and find boundary
% coordinates and actual BC assigned to PDE model
% Enforce initial conditions. Calculate lossIC.
% Find boundary nodes
for i=1:model.Geometry.NumEdges
boundaryNodes = findNodes(model.Mesh, "region", "Edge", i);
% Find interior nodes by excluding boundary nodes
ICiNodes = setdiff(1:size(model.Mesh.Nodes, 2), boundaryNodes);
u0=72201 ;
actualICi =ones(1,numel(ICiNodes))*u0;
IC_XY = model.Mesh.Nodes(:, ICiNodes);
T0IC=zeros(size(ICiNodes));
IC_XYT0IC=[IC_XY;T0IC];

actualIC=actualICi;

end
IC_XYT0IC = dlarray(IC_XYT0IC,"CB"); % Format the coordinates
[predictedIC,~] = forward(net,IC_XYT0IC);
mseIC = l2loss(predictedIC, actualIC);

% Combine the losses. Use a weighted linear combination
% of loss (lossF lossU and lossIC) terms.

loss=mseF+mseIC+mseU;

% Calculate gradients with respect to the learnable parameters.
gradients = dlgradient(loss,net.Learnables);
end