I don't know if this solution is what you are going for - it is more of a kludge, and I can't say if it correctly models diffusion from a center to a regions boundaries.
Method 1 - add an interior boundary and second DirichletCondition
But this is a quick method. First you use the "IncludePoints"
option of ToBoundaryMesh
to include a set of points describing a circular region around Paris. These points are then considered part of the boundary for NDSolve
. Now you add a second DirichletCondition
setting the value of u
to some high constant at that interior boundary.
The two adjustable parameters here are the radius around Paris that you consider, and the value you set the function to at that boundary.
ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
parisradius = 0.3;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
GeoGridPosition[
GeoPosition[
CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
"Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1,
"IncludePoints" -> CirclePoints[paris, parisradius, 50]];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}] - 20;
usol = NDSolveValue[{op == 1,
DirichletCondition[u[x, y] == 0, Norm[{x, y} - paris] > .6],
DirichletCondition[u[x, y] == parisdiffusionvalue,
Norm[{x, y} - paris] < .6]}, u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh,
PlotPoints -> 100, ColorFunction -> "Temperature"],
bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PlotPoints -> 50]
Thanks to kglr for showing how to convert city coordinates to the Mercator projection.
Method 2 - Modify the differential equations
This seems to give a better result, but again it isn't based on any transport theory. All I'm doing is setting the Laplacian equal to a Gaussian function centered around Paris. The only important adjustable parameter is the width of the Gaussian,
ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
parisradius = 2;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
GeoGridPosition[
GeoPosition[
CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
"Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}];
usol = NDSolveValue[{op ==
150 Exp[-Norm[{x, y} - paris]^2/parisradius],
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh,
PlotPoints -> 100, ColorFunction -> "Temperature"],
bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PlotPoints -> 50]