Assigning homogeneous inflow Dirichlet boundary conditions

[rajeevrajaram]

Hello - I am simulating the linear advection equation with a nonlinear autonomous vector field f(x,y) (as opposed to a constant velocity which is usually assumed). For this situation, I need to assign zero Dirichlet boundary values at inflow points on the boundary i.e., those points on the boundary where the vector field points inwards. Mathematically, the inflow points can be found out by picking those points on the boundary that have a negative vector dot product between the vector field f and the unit normal (unormal). So logically, we need an if statement inside the boundary conditions i.e., if the dot product of the two vectors is < 0 for the point under consideration, then assign value(v)=0, otherwise do not assign any boundary conditions. The outflow points for this problem (points where the dot product is > 0) is automatically determined by the solution (since the solution flows out eventually, we cannot specify values on the outflow points, otherwise the problem will be ill-posed).

FlexPDE does not allow the use of if statements withing the boundary condition part of the code. So I would really appreciate it if anybody could help me with this special type of boundary assignment.

Currently, I calculate the inflow points separately on a sheet of paper and then assign zero values to those points. I would like to automate the assignment using an if then else statement.

Thanks.
Rajeev

[moderator]

FlexPDE WILL allow the use of an IF statement within a boundary condition. What it does NOT allow is the BC type to be changed, i.e. from a VALUE to a NATURAL or vice versa. What you will need to do is formulate a NATURAL BC that is equivalent to a VALUE. Then you can use an IF statement to switch between the two BCs.

You can create a NATURAL BC that acts like a VALUE BC by defining a large penalty factor, and applying a NATURAL flux that is proportional to the difference between the current value and the desired value. This will drive the boundary value to the desired value:
NATURAL(v) = penalty*(v0-v)

Then you can switch between these two based on whatever criteria you like:
NATURAL(v) = IF (condition) THEN penalty*(v0-v) ELSE 0

Here is what a script might look like (switching BC at time 1) :

DEFINITIONS
v0 = 0
penalty = 100

BOUNDARIES
START(0,0)
NATURAL(v) = IF (t>1) THEN penalty*(v0-v) ELSE 0
LINE TO ...

Note : the required sign of the penalty term will depend on how your equation is formulated. See the definition of NATURAL flux in the Help documentation. Also, the appropriate size of the penalty term will have to be discovered through trial and error.

[rajeevrajaram]

Thanks. That worked. I have a different question. I was trying out a 3D example and all of a sudden I got an error message that said that Memory allocation was insufficient. I am attaching the code.

{ Fill in the following sections (removing comment marks ! if necessary),
and delete those that are unused.}
TITLE '3D example' { the problem identification }
COORDINATES cartesian3 { coordinate system, 1D,2D,3D, etc }

VARIABLES { system variables }
v (threshold=0.01){ choose your own names }
!Rho (threshold=0.001)
L(threshold=0.01)
SELECT { method controls }
contourgrid=60
!contours = 100
!errlim=1e-3
DEFINITIONS { parameter definitions }
f1=y
f2=z
f3 = -x -2*y-z+x^3
R=0.25
Zsphere=sqrt(max(R^2-x^2-y^2,0))
INITIAL VALUES
v=swage(2.25-x^2-y^2-z^2,0,1,0.001)
L=0
!Rho=1
EQUATIONS { PDE's, one for each variable }
v: dt(v) =- div(f1*v,f2*v,f3*v)+0.01*div(grad(v)){ one possibility }
!Rho:dt(Rho)=v
L: dt(L)=exp(-L)*v
! CONSTRAINTS { Integral constraints }
Extrusion
Surface 'bottom of the cube' z=-1
Layer 'below the sphere'
Surface 'bottom of the sphere' z=-max(Zsphere,0)
Layer 'inside the sphere'
Surface 'top of the sphere' z=max(Zsphere,0)
Layer 'above the sphere'
Surface 'top of the cube' z=1
BOUNDARIES { The domain definition }
REGION 1 'cube' { For each material region }
START (-1,-1) value(v)=0 line to (1,-1)
value(v)=0 line to (1,1)
value(v)=0 line to (-1,1)
value(v)=0 line to close
Surface 'bottom of the cube' value(v)=0
Surface 'top of the cube' value(v)=0
exclude limited region 2 'sphere'
Layer 'inside the sphere'
start (0.25,0)
nobc(v)
arc(center=0,0) angle=360
TIME 0 by 0.01 TO 1{ if time dependent }
MONITORS { show progress }
PLOTS { save result displays }
FOR T = 0 by 0.01 to 1
contour(arctan(L)) on 'cube' on z=0 painted contours=100
!contour(if (Rho>0) then log10(Rho) else 0) on 'outside' notags contours=2000 painted range=(globalmin(log10(Rho)),globalmax(log10(Rho)))
!surface(if Rho>0 then log10(Rho) else 0) painted
END

[rajeevrajaram]

I forgot to mention that I have 232 GB of free memory on my computer. In the example, I am trying to simulate the advection equation on a cube centered at origin and exclude a small sphere centered at the origin (because the solution is singular at the origin).
Thanks.
Rajeev

[moderator]

Because of your discontinuous initial condition on V around the sidewalls (initial value=1, BC=0), FlexPDE refines the grid on these sidewalls quite aggressively, resulting in about 600,000 nodes in the mesh. FlexPDE reports a memory usage of 2.2GB, but apparently because of memory space fragmentation, the actual memory space allocation is closer to 5GB (as reported by the Windows Task Manager) . Your quoted 232GB free memory refers to disk space, not RAM space. If you are running a 32-bit version of FlexPDE, the maximum RAM space available is 2GB (this is a hardware imposed limit). If you are running a 64-bit version of FlexPDE, then the limit will depend on how much physical RAM you have installed and whether you have virtual memory enabled.
The best thing to do is NOT to impose discontinuous initial conditions. You have constructed a large sphere of initial value 1, but this sphere is larger than your box, so the initial value is 1 almost everywhere. The SWAGE you have used is very sharp (0.001 transition width), so even if it were contained in the box, you would still have a sharp demarcation and a cause for dense gridding. You should re-think what your initial conditions really are.

[rajeevrajaram]

THank you so much for that. I changed my initial conditions, and now it runs (albeit taking a long time to complete). I am guessing that I have to expect this for 3D simulations depending on how hard the problem is.
Rajeev

[rajeevrajaram]

Dear moderator - I saw your reply to another post regarding 3D visualizations. I am able to successfully obtain contour plots on cut planes of a cube. I know that there is a way to export data from FlexPDE to other visualization software. I am guessing it is possible to view the flow on the whole cube?? I am reading up on the export command from the manual, and I think I know how to export the coordinates of points and the value of the solution to a file. But I have never used the visulaization softwares that are mentioned in the manual. I know how to use MATLAB. Is it possible to use the plot routine in MATLAB to see the whole cube after exporting data from FlexPDE?
Thanks.
Rajeev