Nonlinear material simulation in pdetool

I tried to solve the similar problem explained here with a nonlinear material. And I wrote this code (with a simpler geometry): (And I defined the c coefficient with an interpolation function with respect to the magnitude of the gradient of the solution|grad(solution)|)

u0=4e-7*pi;
model = createpde(1);
gd = [  3,4,+2.0,-2.0,-2.0,+2.0,1.5,1.5,-1.5,-1.5;...
        3,4,+0.8,-0.8,-0.8,+0.8,0.8,0.8,-0.8,-0.8;...
        3,4,+0.4,-0.4,-0.4,+0.4,0.4,0.4,-0.4,-0.4;...
        3,4,+1.2,+0.8,+0.8,+1.2,0.4,0.4,-0.4,-0.4;...
        3,4,-0.8,-1.2,-1.2,-0.8,0.4,0.4,-0.4,-0.4;]'; % geometry gd
ns = char('bound','c1','c2','c3','c4')';
sf = 'bound+c1+c2+c3+c4';
[ dl, bt] = decsg(gd,sf,ns);
model.SolverOptions.ReportStatistics = 'on';
geometryFromEdges(model,dl); % include geometry
generateMesh(model,'Hmax',0.04,'GeometricOrder','Linear');
applyBoundaryCondition(model,'dirichlet','Edge',[1,2,13,14],'u',0);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0,'Face',1);
cCoef = @(~,state) interp1([0 1-4 2e-4 2],[1/(1000) 1/(500) 1/(100) 1/(10)],sqrt((state.ux).^2 + (state.uy).^2));
specifyCoefficients(model,'m',0,'d',0,'c',cCoef,'a',0,'f',0,'Face',2);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',u0,'Face',3);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0,'Face',4);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0,'Face',5);
results = solvepde(model);
B=sqrt((results.XGradients.^2+results.YGradients.^2)); %B=curl(A);
pdeplot(model.Mesh.Nodes,model.Mesh.Elements,'XYData',B,'Mesh','on');
figure; quiver(results.Mesh.Nodes(1,:),results.Mesh.Nodes(2,:),results.YGradients',-results.XGradients');

I supposed that the solver will handle the problem iteratively, however, I got this instead:

Iteration     Residual     Step size  Jacobian: Full
   0          1.0842e-19

And the solution is exactly the same as the time I replace the c coefficient with 1/1000. I expected to see some sort of saturation in the material.

[Alan Weiss]

I believe what occurs is that the residual is very small, so the solver thinks that it is done, that the initial solution is close enough to the correct one. But this comment is based only on viewing the reported iteration, not understanding the problem.

So my question for you is, did the solver find a good solution?

Alan Weiss

MATLAB mathematical toolbox documentation

[Hamed Eskandari]

Dear Alan

Thank you for your comment. You were right about my miscalculation. The initial solution was to close to the final solution and that's why the iterative solution terminated in the first iteration.

I fixed the problem by increasing the excitation value. And it worked pretty well and accurate. Here is the final code:

u0=4e-7*pi;
model = createpde(1);
gd = [  3,4,+2.0,-2.0,-2.0,+2.0,1.5,1.5,-1.5,-1.5;...
        3,4,+0.8,-0.8,-0.8,+0.8,0.8,0.8,-0.8,-0.8;...
        3,4,+0.4,-0.4,-0.4,+0.4,0.4,0.4,-0.4,-0.4;...
        3,4,+1.2,+0.8,+0.8,+1.2,0.4,0.4,-0.4,-0.4;...
        3,4,-0.8,-1.2,-1.2,-0.8,0.4,0.4,-0.4,-0.4;]'; % geometry gd
ns = char('bound','c1','c2','c3','c4')';
sf = 'bound+c1+c2+c3+c4';
[ dl, bt] = decsg(gd,sf,ns);
model.SolverOptions.ReportStatistics = 'on';
geometryFromEdges(model,dl); % include geometry
generateMesh(model,'Hmax',0.04,'GeometricOrder','Linear');
applyBoundaryCondition(model,'dirichlet','Edge',[1,2,13,14],'u',0);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0,'Face',1);
cCoef = @(~,state) interp1([0 1-4 2e-4 2],[1/(1000) 1/(500) 1/(100) 1/(10)],sqrt((state.ux).^2 + (state.uy).^2));
specifyCoefficients(model,'m',0,'d',0,'c',cCoef,'a',0,'f',0,'Face',2);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',u0*10000,'Face',3);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0*10000,'Face',4);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0*10000,'Face',5);
results = solvepde(model);
B=sqrt((results.XGradients.^2+results.YGradients.^2)); %B=curl(A);
pdeplot(model.Mesh.Nodes,model.Mesh.Elements,'XYData',B,'Mesh','on');
figure; quiver(results.Mesh.Nodes(1,:),results.Mesh.Nodes(2,:),results.YGradients',-results.XGradients');

[Viktor Szuhai]

Nice job. But I believe I found a small typo: an "e" is missing in the cCoef definition. Correct:

cCoef = @(~,state) interp1([0 1e-4 2e-4 2],[1/(1000) 1/(500) 1/(100) 1/(10)],sqrt((state.ux).^2 + (state.uy).^2));