# Tutorial 17 – Iterative Solvers

This tutorial is intended for the full version of the toolbox.

## Introduction

In this tutorial we will present iterative solver in microfuidic application. For more detail about microfluid please refer to tutorial 16.

## Code

```% Clear MATLAB Command Window and Workspace
clc;
clear;

[p,e,t] = importMeshGmsh('microfluidicH3D_coarse.msh');
%[p,e,t] = importMeshGmsh('microfluidicH3D.msh');
displayMesh(p,t);
```

```% Convert mesh to second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);

% Define fluid
nu = 0.1;

% Init solution and convergence criteria
[u, convergence] = initSolution(p,t,[0 0 0],0);
maxres = 1e-3;
maxiter = 25;

% Define inlet velocity profile
vel = [1 0 0];

% Options of the iterative solver
opt.alpha_p = 0.02;
opt.pres_iter = 8;
opt.bicg_iter = 48;
opt.bicg_eps = 1e-10;

% Iterate nonlinear terms
for iter = 1:maxiter
% Compute shear rate and viscosity for non-Newtonian fluids
% gamma = shearRate(p,t,u);
% nu = nuinf + (nu0-nuinf)*(1+(lambda*gamma).^2).^((n-1)/2);

% Assemble matrix and right hand side
[NS, F] = assembleNavierStokesMatrix(p,e,t,nu,u(indices.indu),u(indices.indv),u(indices.indw),'nosupg');
% [NS, F] = assembleNavierStokesMatrix(p,e,t,nu,u(indices.indu),u(indices.indv),u(indices.indw),'supgDoublyAsymptotic');

% Apply boundary conditions
[NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,[77 78],'inlet', [1 0 0]);
[NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,[80 81],'wall');

% Compute and plot residuals
[stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres);
plotResiduals(convergence,2);

% Break if solution converged
if(stop)
break;
end

% Solve equations
tic
%[u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'direct');
%[u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'SIMPLE');
[u, convergence] = iterativeSolver(NS, F, u, indices, convergence, 'SIMPLE',opt);
toc
end

% Plot solution
figure(3);
displaySolution(p,t,u(indices.indu),'x-velocity');
```

```% Solve scalar transport problem
k = 0.001;
%k = 0.1;

% Assemble problem matrix
[D,F] = assembleDiffusionMatrix(p,t,k);
D = D + assembleScalarConvectionMatrix(p,t,k,u(indices.indu),u(indices.indv),u(indices.indw),'supgDoublyAsymptotic');

% Apply boundary conditions
[D,F] = imposeScalarBoundaryCondition(p,e,D,F,77,'value',0);
[D,F] = imposeScalarBoundaryCondition(p,e,D,F,78,'value',1);

% Solve system of equations
phi = D\F;
% Display solution

figure(4);
displaySolution(p,t,phi,'Concentration');
```

```figure(5);
pressure = generatePressureData(u,p,t);
displaySolution(p,t,pressure,'Pressure');
```