We start by creating a 1 D mesh
In[]:=
Needs["NDSolve`FEM`"]region=Line[{{0},{1}}];includePoints={{1/3},{2/3}};mesh=ToElementMesh[region,"IncludePoints"->includePoints,"MaxCellMeasure"->0.0008]
Out[]=
ElementMesh[{{0.,1.}},{LineElement[<1248>]}]
The mesh has points included at 1/3 and 2/3. This is where we will position the point sources later . I also refined the mesh, because I am going to add a convection term which will make this PDE convection dominant (large Peclet number) . We set up the dependent variable, the mass concentration c, time variable t and the spatial independent variable x . Also some material properties are defined .
In[]:=
vars={c[t,x],t,{x}};pars=<|"DiffusionCoefficient"->1,"MassConvectionVelocity"->{100}|>;
For the point sources the idea is to use a regularized delta function . The Dirac delta function poses a problem in numerical simulations with FEM implementation that is not prepared to take it as an input . This is because the Dirac delta function is singular at the source location and needs special treatment that differs from how FEM treats other terms . This can be avoided with an approximation to the Dirac delta function . The process of approximating the Dirac delta function is called regularization .
X
s
In[]:=
RegularizedDeltaPoint[g_,X_List,Xs_List]:=Piecewise[{{Times@@Thread[1/(4g)(1+Cos[π/(2g)(X-Xs)])],And@@Thread[RealAbs[X-Xs]<=2g]},{0,True}}]
More information about this can be found in the Heat Transfer tutorial (www.reference.wolfram.com/language/PDEModels/tutorial/HeatTransfer/HeatTransfer.html#945943234) . The process for mass transport is the same. So if you also want to volume sources you can proceed in the same way as mentioned in that tutorial .
In[]:=
Subscript[h,mesh]=Sqrt[Min[mesh["MeshElementMeasure"]]];Subscript[gamma,reg]=Subscript[h,mesh]/2;
The regularization function needs a parameter gamma, which is based on the minimal element diameter . Let' s construct the point source at the first included point and visualize it .
In[]:=
temp=RegularizedDeltaPoint[Subscript[gamma,reg],{x},includePoints[[1]]];Plot[temp,{x,0,1}]
Out[]=
You see that this is an idealized approximation to the delta distribution . Note that
In[]:=
Integrate[temp,{x,0,1}]
Out[]=
1.
The Comsol approach uses a weak formulation . That is a different approach . I am trying to show that (point) sources can be set up in different ways .
In[]:=
Qp=10;pars["MassSource"]=RegularizedDeltaPoint[Subscript[gamma,reg],{x},includePoints[[1]]]*Qp+RegularizedDeltaPoint[Subscript[gamma,reg],{x},includePoints[[2]]]*If[t>1/2*10^-3,5*Qp,0];
Qp is the point source strength . The unit would be [mol/s]. The first point source is active all the time while the second only starts at a specific time but is stronger . We set up the PDE and solve it .
In[]:=
pde={MassTransportPDEComponent[vars,pars]==0,c[0,x]==0};tEnd=10^-3;cfun=NDSolveValue[pde,c,{t,0,tEnd},{x}∈mesh];
Visualize/explore the result .
In[]:=
Manipulate[Plot[cfun[t,x],{x}∈region,PlotRange->{{0,1},{0,1/2}}],{{t,2tEnd/3},0,tEnd},SaveDefinitionsTrue]
Out[]=
