Comparison PDE solve and \ mldivide

Hello all,

I have a problem and would like you thoughts on it. I guess I miss something in the definition of some matlab functions but can not find out where my error is.

Comparison of solve function with mldivide in a transient heat analysis gives different solutions.

Probably an missuse of the Nullspace transformation but did not find yet where??

I would really like to understand what I got wrong and would be thankfull for your comments.

Phineas Freak

%Following unified modelling approach:
clear all
%Geometry
shape_1=polyshape([-150, 150, 150, -150],[-150, -150, 150,150]);
BlockRaw= triangulation(shape_1);
BlockRaw=fegeometry(BlockRaw);
BlockRaw=extrude(BlockRaw,300);
BlockRaw=scale(BlockRaw, 1/1000); % [mm] to [m]
% Model
model=femodel(AnalysisType="ThermalTransient", Geometry= BlockRaw);
figure
pdegplot(model,FaceLabels="on",EdgeLabels="on")

model.MaterialProperties=materialProperties(ThermalConductivity=46.5, MassDensity=7800, SpecificHeat=.502);
% Initial Temperature
model.CellIC(1) = cellIC(Temperature=273+180);
%BC1
model.FaceLoad(3) = faceLoad(Heat=300);
% BC2
model.FaceBC(1) = faceBC(Temperature=200);
% BC 3
model.FaceLoad(6) = faceLoad(Heat=100);
% Coarse Mesh - exact result not important as comparison between methods
model = generateMesh(model, GeometricOrder='linear', Hmin=.05);
figure
pdemesh(model)

% Time Interval
dt=.25;
t_end=4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Builtin solver
ZylinderRes = solve(model, [0:dt:t_end]);
figure
pdeplot3D(ZylinderRes.Mesh,ColorMapData=ZylinderRes.Temperature(:,end))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mldivide
state = struct('time', ZylinderRes.SolutionTimes(1), 'u', ZylinderRes.Temperature(:,1)');
FEMMatrix=assembleFEMatrices(model,state);
K = FEMMatrix.K;
M=  FEMMatrix.M;
F = FEMMatrix.F;
H=FEMMatrix.H;
R=FEMMatrix.R;
G=FEMMatrix.G;
theta=1;
[B,Or]=pdenullorth(FEMMatrix.H);
ud=Or*((H*Or\R));   % Prüfung - Korrekt

U_=zeros(length(ud),17);
U_(:,1)=ZylinderRes.Temperature(:,1);

% Nullspace Dirichlet Transformation
Kc =B'*(theta*K + M/dt)*B;

dt_sum = 0;
for i=2:t_end/dt+1
    dt_sum=dt_sum+dt;
    Fc=B'*((F+G) - K*ud + ( M/dt - (1-theta)*K )*U_(:,i-1))    ;
    U_(:,i) =  B* (Kc \ Fc ) + ud;
end
figure
pdeplot3D(ZylinderRes.Mesh,ColorMapData=U_(:,end))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Comparison: Missmatch
figure
plot(ZylinderRes.Temperature(1:400,4)-U_(1:400,4),'b-')

[Torsten]

Are you sure that "solve" uses a uniform time stepsize of 0.25 or is the vector [0:dt:t_end] only used to fix the times when a solution should be supplied ?

Where did you find the method to compute U_ ?

[Phineas]

Hello Torsten,

Thanks for the redaction of the question, looks smooth.

1. I can not edit the solve function, its a MathWorks built in function. I took the proceeding for solve from the examples.
2. Calculation of U_ follows different literature sources that are quite coherent (except me making an error). [1] , [2] , [3].

Problem here is, I do not know how M from assembleFEMatrices() is calculated, Matlab manual is somehow ambiguous here. Can be coefficient for d2T/dt2 or dT/dt term: "M is the mass matrix, the integral of the discretized version of the m or d coefficients. " --> Section ALGORITHMS

[Phineas]

In fact, I will try solvepde() might be this changes the solution....

Answer 1

[Phineas]

This seems to work.

clear all

dt=.25;
t_end=4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
model2=createpde();
model2.Geometry=multicuboid(.3,.3,.3);
pdegplot(model2.Geometry,FaceLabels="on",FaceAlpha=0.5)
axis equal
model2.generateMesh(GeometricOrder='linear',Hmin=.05);
setInitialConditions(model2,273+180);
applyBoundaryCondition(model2,"dirichlet", ...
    "Face",1,"u",100)
applyBoundaryCondition(model2,"neumann",Face=[3],q=0,g=200);
specifyCoefficients(model2,m=0,...
    d=1,...
    c=1,...
    a=0,...
    f=1);
Result2 = solvepde(model2, [0:dt:t_end]);
pdeplot3D(model2,"ColorMapData",Result2.NodalSolution(:,5))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Own Solver
FEMMatrix2=assembleFEMatrices(model2)

K = FEMMatrix2.K;
M =  FEMMatrix2.M;
H =  FEMMatrix2.H;
R = FEMMatrix2.R;
G = FEMMatrix2.G;

%theta=1;
[B,Or]=pdenullorth(H);
ud=Or*((H*Or\R));   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
U_=zeros(length(ud),t_end/dt+1);
U_(:,1)=Result2.NodalSolution(:,1);

theta=1;

Kc =B'*(theta*K + M/dt)*B;

dt_sum=0;
for i=2:t_end/dt+1
    dt_sum=dt_sum+dt;
    Fc=B'*( G + ( M/dt -(1-theta)*K )*U_(:,i-1) - K*ud)    ;
    U_(:,i) =  B* (Kc \ Fc ) + ud;
end

for i=1:t_end/dt+1
    mean(Result2.NodalSolution(:,i)-U_(:,i))
end

for i=1:t_end/dt+1
    figure
    plot(Result2.NodalSolution(1:400,10),U_(1:400,10),'b*')
end

Answer 2

[Walter Roberson]

model.MaterialProperties(ThermalConductivity=46.5, MassDensity=7800, SpecificHeat=.502);

That is invalid. It is converted internally into

model.MaterialProperties("ThermalConductivity", 46.5, "MassDensity", 7800, "SpecificHeat", .502);

which attempts to use the string "ThermalConductivity" as an index into the model.MaterialProperties array.

You instead need something like

model.MaterialProperties = materialProperties(ThermalConductivity=46.5, ...
    MassDensity=7800, ...
    SpecificHeat=0.502);

[Phineas]

Hi Walter,

Thanks for your remark, I will adapt the script.

But, changing the parameters had an effect on the simulation. And as these parameters go into mass and stiffness matrices this should not have an effect on the convergence.

Kind Regards

P. Freak

[Phineas]

I changed it but for the moment still the gap. Will switch doe solvepde().

[Walter Roberson]

The call you had was simply invalid. The replacement code uses the properties that you were trying to specify, not different properties.