I would like to solve an advection-diffusion problem on a torus domain. There are three Dirichlet conditions: One at the inlet (concentration $c=0$), one at the outlet ($c=0.5$) and one at the wall ($c=1$).

The geometry is set up in the following code. Note that RegionCapillary is the full solution domain, RegionCapillary2D is a 2D projection for convenient plotting and the list points goes through the centre of the tube. In the plot below RegionCapillary is orange and RegionCapillary2D is blue.

radiusMajor = 4;
radiusMinor = 1;
RegionCapillary = 
  ImplicitRegion[( radiusMajor - Sqrt[x^2 + y^2])^2 + z^2 <= 
     radiusMinor^2 && x >= 0, {x, y, z}];
RegionCapillary2D = 
  ImplicitRegion[( radiusMajor - Sqrt[x^2 + y^2])^2 <= radiusMinor^2 &&
     x >= 0, {x, y}];
RP = RegionPlot3D[
   DiscretizeRegion /@ {RegionCapillary, RegionCapillary2D} // 
    Evaluate, Axes -> True, AxesLabel -> {x, y, z}, 
   PlotStyle -> {Opacity[0.4], Opacity[1]}];
LP = ListPointPlot3D[
   points = 
    Reverse@Table[{radiusMajor Cos[\[Phi]], radiusMajor Sin[\[Phi]] , 
       0}, {\[Phi], -Pi/2, Pi/2,  Pi/1000.0}], 
   PlotStyle -> PointSize[0.03]];
Show[RP, LP]

enter image description here

I solve the problem very similarly to user21's reply to my question (here). There are two parts to the problem: Stokes flow and then an advection-diffusion in which the previous (Stokes) solution provides the fluid velocity.

The Stokes problem is:

<< NDSolve`FEM`
StokesOp := {
   -Inactive[Laplacian][u[x, y, z], {x, y, z}] + D[p[x, y, z], x],
   -Inactive[Laplacian][v[x, y, z], {x, y, z}] + D[p[x, y, z], y],
   -Inactive[Laplacian][w[x, y, z], {x, y, z}] + D[p[x, y, z], z],
   Div[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}]};

inletCondition[x_, y_, z_] := x == 0 && y >= 0
outletCondition[x_, y_, z_] := x == 0 && y < 0
capillaryCondition[x_, y_, z_] := x > 0

bcs := {
   DirichletCondition[p[x, y, z] == 1, inletCondition[x, y, z]],
   DirichletCondition[p[x, y, z] == 0, outletCondition[x, y, z]],
   DirichletCondition[{u[x, y, z] == 0., v[x, y, z] == 0., 
     w[x, y, z] == 0.}, capillaryCondition[x, y, z]]};

VelMMA[x_, y_, z_] := Sqrt[
 xVel[x, y, z]^2 + yVel[x, y, z]^2 + zVel[x, y, z]^2]

Clear[xVel, yVel, pressure];
AbsoluteTiming[{xVel, yVel, zVel, pressure} = 
   NDSolveValue[{StokesOp == {0, 0, 0, 0}, bcs}, {u, v, w, 
     p}, {x, y, z} \[Element] RegionCapillary, 
    Method -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 5 10^-3}}, 
    "ExtrapolationHandler" -> {0. &, "WarningMessage" -> False}];]

The advection-diffusion problem is:

Pe = 0.7*10^4;

eqn = Pe {xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]}.Grad[
      c[x, y, z], {x, y, z}] - Laplacian[c[x, y, z], {x, y, z}];

bcs = {
   DirichletCondition[c[x, y, z] == 0, inletCondition[x, y, z]],
   DirichletCondition[c[x, y, z] == 0.5, outletCondition[x, y, z]],
   DirichletCondition[c[x, y, z] == 1, capillaryCondition[x, y, z]]

 conMMA = NDSolveValue[{eqn == 0, bcs}, 
    c, {x, y, z} \[Element] RegionCapillary, 
    Method -> {"FiniteElement", "InterpolationOrder" -> {c -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 10^-3}}, 
    "ExtrapolationHandler" -> {0. &, "WarningMessage" -> False}];]

The velocity magnitude and concentration are plotted like this:

 Sqrt[xVel[x, y, 0]^2 + yVel[x, y, 0]^2 + zVel[x, y, 0]^2], {x, 
   y} \[Element] DiscretizeRegion@RegionCapillary2D, PlotRange -> All,
  PlotLabel -> "velocity magnitude", PlotPoints -> 100, 
 ColorFunction -> "TemperatureMap", AspectRatio -> 2, 
 PlotLegends -> Automatic]
 conMMA[x, y, 0], {x, y} \[Element] 
  DiscretizeRegion@RegionCapillary2D, PlotRange -> All, 
 PlotLabel -> "concentration", PlotPoints -> 100, 
 ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, 
 AspectRatio -> 2]

enter image description here

I have compared these results with an identical problem setup in COMSOL, finding that COMSOL's solution looks smoother and more plausible: enter image description here

If you trace the velocity and concentration along the centreline which goes right through the tube, you can see that while both solutions satisfy $c=0.5$ at the outlet (i.e. the right hand side of the ListPlots), Mathematica has a very weird behaviour for the diffusion problem. The Stokes problem is largely the same between the two software packages: enter image description here I am puzzled by this discrepancy and would very much appreciate some help. My suspicion is that I am specifying the boundary conditions in Mathematica incorrectly.


With a finer mesh on both Mathematica and COMSOL, the solutions do converge somewhat, but the behaviour near the outlet of the Mathematica solution is still somewhat puzzling: enter image description here

EDIT 2: Do things get better for lower Péclet number?

Yes. I made a comparison for Pe=10^3(original version is Pe = 0.7*10^4), finding a better convergence and less weirdness near the boundary:

enter image description here enter image description here

The weird effects at the boundary of the DensityPlot are just interpolation artifacts from importing Comsol results from text file.

EDIT 3: A more complex example

OK, if I'm doing this, might as well go all in. Here is a more complex example of the tube, based on a looping segment with quite irregular cross-sections. Here it is plotted with centreline as well as planes which I use in inequalities to set boundary conditions in Mathematica:

enter image description here

A major difference to the above examples is that at the outlet there is a Neumann condition, $\boldsymbol{n}\cdot \nabla c = 0$, which in Mathematica is NeumannValue[0, pred]. Here is a plot of concentration and velocity accross the centreline, just as in the original examples:

enter image description here

You can see that once again the Stokes problem is a superb match, but the advection-diffusion problem has some serious discrepancy between solutions. Still, it is interesting that the curves are so similar and only appear to be offset. The negative concentration produced near the inlet by Mathematica is however unphysical.

Here is a pretty visualisation of with SliceContourPlot, which shoes how similar the cross-section profiles look:

enter image description here

EDIT 4: A conclusion on 3D FEM with imported geometries

I ended up working a lot on similarly complex imported 3D geometries as above, switching completely to Comsol was essential. The boundary conditions were getting so tricky it really helped having a visual interface. The meshing control on boundary layers was also helpful, and their constructive geometry (boolean operations on discrete regions) works ok. I ended up being very careful about ensuring flux conservation in Comsol (The best way to evaluate advective and diffusive flux is via the tds.ncflux_c and tds.ndflux_c). My simulation results to similar problems to edit 3, made in Comsol and checked very carefully, can be seen in papers here and here. In fact, in the first linked paper we even looked fluxes at cross-sections of sub-meshes ($N_j$) of quite complex vasculature, so a visual interface was enormously helpful for that:

enter image description here´

As for Mathematica, I gave up on trying to solve complex advection-diffusion problems on discrete 3D domains, simply because the meshing got too complicated. If you import a discrete BoundaryMesh in MMA, the boundary cannot easily refined: As it does not come from predefined MMA regions, options like MaxCellMeasure and AccuracyGoal cannot be used on the boundary (related question here). Maybe the workflow with an external mesher like gmsh is the way to go here, but I haven't tried.

  • 1
    $\begingroup$ Does the behavior change if you use a smaller Pe? $\endgroup$
    – user21
    Commented Nov 9, 2017 at 10:27
  • $\begingroup$ Yes, see Edit 2. $\endgroup$ Commented Nov 9, 2017 at 11:50
  • $\begingroup$ I added another more complex example, perhaps it is of interest for you @user21. $\endgroup$ Commented Nov 9, 2017 at 12:28
  • 1
    $\begingroup$ Thank you! I ran some more simulations with low Peclet numbers and the agreement is generally very good. I think that the fine boundary layer structure of the high Pe case cannot be resolved by the relatively coarse mesh I chose. But the lower Pe, the less of a problem thin boundary layers become. I am hoping to write a proper answer to my question with the low Pe plots soon. $\endgroup$ Commented Nov 14, 2017 at 12:16
  • 1
    $\begingroup$ If you could briefly outline how you addressed the issue that would be great. It's nearly impossible to answer your question without having access to COMSOL and to see what caused the difference between the models. $\endgroup$
    – user21
    Commented Jan 2, 2018 at 7:16

1 Answer 1


enter image description hereI have a suggestion how to set the correct outlet boundary condition (I think this is corrupting the entire solution... at least in Comsol). The point is that the advection diffusion equation can be written as $$ \nabla \cdot J = 0, $$ with $$ J=v\,c - D\nabla c $$ where $v$ is the velocity field and $D$ is the diffusion coefficient.

The bulk lagrangian density (weak form) is $$ \mathcal{L}_\mathrm{bulk}=-\nabla c^\dagger\cdot J $$ where $c^\dagger$ is the test functtion, and the boundary contribution is $$ \mathcal{L}_\mathrm{boundary}=c^\dagger J\cdot\hat{n} $$ this is the weak form contribution to use in place of the outlet boundary condition. This avoid the ripples for the concentration field at the outlet. Here is the solution with a computed velocity field for $Pe=100$ (almost advection dominated problem)

  • $\begingroup$ In Comsol I am using the transport of dilute species package, which has a predefined type "outlet" which is $n\cdot\nabla c=0$. I think in Mathematica this should simply be taken care of by NeumannValue[0, pred]. $\endgroup$ Commented Nov 14, 2017 at 20:39
  • 2
    $\begingroup$ You can add a "weak contribution" in comsol and obtain the results I posted... isn't this becoming a comsol-related discussion?! :) $\endgroup$
    – Fabio
    Commented Nov 14, 2017 at 20:42
  • $\begingroup$ I'd be quite interested to have a look at your Comsol file, feel free to email me at [email protected] if you like. $\endgroup$ Commented Nov 14, 2017 at 20:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.