| Author | Message | ||
     Andrea Cova (andrea25) Member Username: andrea25 Post Number: 7 Registered: 05-2006 |
Hello,I have written a code to simulate the problem I'm studying but I still have some problems,since the monitors during the execution don't seem to show anything.The domain seems to be correctly defined and a correct mesh has been done.Could you give me some help?Is the problem with the code or with the way in which I request monitors and plots? Thank you very much Andrea __________________________________________________ title 'Modellizzazione del disastro ecologico di Bhopal' coordinates cartesian3 variables c(threshold=0.1) definitions u=4 Ky=[0.0121*X^2*(1+0.0004*X)^(-1)*u]/(2*X) Kz=[0.0064*X^2*(1+0.00015*X)^(-1)*u]/(2*X) a=90 Q0=40 hmin=39.6 Hmax=200 ampx=100 semiampy=100 cmax=100 initial values c = 0 equations Ky*DYY(c)+Kz*DZZ(c)-u*DX(c)-Dt(c)=0 extrusion surface 'livello terrestre' Z=0 layer 'spessore del rilievo orografico' surface 'livello superiore del rilievo' Z=80 layer 'spessore dell"aria sovrastante' surface 'estremo superiore del dominio in esame' Z=Hmax boundaries region 1 'intero dominio spaziale' start(0,-semiampy) value(c)=0 line to (ampx,-semiampy) natural(c)=0 line to (ampx,semiampy) value(c)=0 line to (0,semiampy) natural(c)=0 line to close region 2 'regione del rilievo orografico' start(50,-20) line to (90,-20) to (90,20) to (50,20) to close time 0 by 1 to a monitors for cycle=1 contour(c) on z=0 as "XY Concentration" range=(0,cmax) contour(c) on x=0 as "YZ Concentration" range=(0,cmax) contour(c) on y=0 as "XZ Concentration" range=(0,cmax) elevation(c) from (0,0,0) to (ampx,0,0) as "X-Axis Concentration" range=(0,cmax) elevation(c) from (0,-semiampy,0) to (0,semiampy,0) as "Y-Axis Concentration" range=(0,cmax) elevation(c) from (0,0,0) to (0,0,Hmax) as "Z-Axis Concentration" range=(0,cmax) plots for t = endtime contour(c) on z=0 as "XY Concentration" range=(0,a) contour(c) on x=0 as "YZ Concentration" range=(0,a) contour(c) on y=0 as "XZ Concentration" range=(0,a) end | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 8 Registered: 05-2006 |
I also added a condition at the end of the boundaries section for a source in a specified point: boundaries region 1 'intero dominio spaziale' start(0,-semiampy) value(c)=0 line to (ampx,-semiampy) natural(c)=0 line to (ampx,semiampy) value(c)=0 line to (0,semiampy) natural(c)=0 line to close region 2 'regione del rilievo orografico' start(50,-20) line to (90,-20) to (90,20) to (50,20) to close fixed point(0,0,hmin) point value(c)=Q0 | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 633 Registered: 06-2003 |
Neither of your setups include any source of C. There are no source terms in the PDE, and the boundary conditions are either zero flux or zero value. So nothing is there to cause C to differ from zero. There must be a source of C somewhere in your system. Either it flows in through a boundary, or it is created internally. You must define what this source is. You don't want to use a fixed point. This is a 2D construct. | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 9 Registered: 05-2006 |
I thought that the problem was the lack of a source for the concentration c.But in the problem I'm trying to analyze the sorce is localized only in a fixed point,since the emission of pollutant happens exclusively at the point x=0,y=0,z=hmin. How can I introduce this souce condition in my flexpde code?Thank you very much! | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 635 Registered: 06-2003 |
Point values are non-physical, and in an adaptive-mesh system like FlexPDE and lead to intense regridding. You should program a source with some physical extent on the bottom surface. There are several ways to do this. 1) define a LIMITED REGION on the bottom surface, give it a small box area and assign the desired VALUE condition on the box. This will still cause dense regridding, because of the discontinuous boundary condition (jumping from zero to nonzero at the box edge). 2) define a narrow gaussian or super-gaussian VALUE condition on the bottom surface. This shape has a smooth transition and allows the error estimator to see the source even if it does not fall directly on a mesh node. 3) define the box of (1) with the source distribution of (2). | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 10 Registered: 05-2006 |
I have tried to follow your hints and I have redifined the BOUNDARIES section in the following way: _________________________________________________ extrusion surface 'livello terrestre' Z=0 layer 'spessore del rilievo orografico' surface 'livello superiore del rilievo' Z=80 layer "spessore dell'aria compresa tra sommità del rilievo e sorgente di emissione' surface 'livello della sorgente di emissione' Z=hmin layer 'spessore dell"aria sovrastante alla sorgente di emissione' surface 'estremo superiore del dominio in esame' Z=Hmax boundaries region 1 'intero dominio spaziale' start(0,-semiampy) value(c)=0 line to (ampx,-semiampy) natural(c)=0 line to (ampx,semiampy) value(c)=0 line to (0,semiampy) natural(c)=0 line to close region 2 'regione del rilievo orografico' start(50,-20) line to (90,-20) to (90,20) to (50,20) to close limited region 3 'regione della sorgente di emissione' start(0,-0.5) line to (1,-0.5) to (1,0.5) to (0,0.5) to close surface 3 value(c)=Q0/u _________________________________________________When I try to run the code,after some minutes of computation, flexpde gives the following error: Ambiguous common face(207,206,209) Ambiguous common face(207,209,206) Ambiguous common face(209,207,208) Ambiguous common face(208,207,209) Ambiguous common face(207,209,206) and then: Memory protection fault ---called from domain::find_face_bc ---called from operate::tabulate_bcs ---called from loading::map_matrix ---called from timsolve::evolve ---called from control::do_runjob What's the problem?How can I resolve it? Thank you very much! | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 639 Registered: 06-2003 |
Have you used the "Domain Review" button to watch the domain construction process? | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 11 Registered: 05-2006 |
Yes,and it seems correctly done!Even when I run the script,the domain is correctly realized and then begins the computation.Frequently during the first cycle on the left it appears the writing:STEP FAILED.Then it seems to work but,after 10 or 11 cycles the program stops ang gives the error I wrote in the previous post | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 12 Registered: 05-2006 |
I have tried to solve the problem but every attempt seems to be useless.Even the technique of defining a gaussian value doesn't seem to work,since I obtain the same error message I wrote some posts above. Can you please help me in defining the sorce condition for this problem? Thank you very much! | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 640 Registered: 06-2003 |
The extrusion specification you posted does not work, because hmin is less than 80 and surface 3 is below surface 2. What values of hmin and hmax did you use? | ||
| Andrea Cova (andrea25) Member Username: andrea25 Post Number: 13 Registered: 05-2006 |
Sorry Robert but I changed those parameters in the last code I am attempting to modify.Here's the code I'm actually working on: title 'Modellizzazione del disastro ecologico di Bhopal' coordinates cartesian3 variables c(threshold=0.1) definitions u=4 { velocità di trasporto del vento pari a 4m s^-1} Ky=[0.0121*X^2*(1+0.0004*X)^(-1)*u]/(2*X) { coefficiente di diffusione lungo l'asse y} Kz=[0.0064*X^2*(1+0.00015*X)^(-1)*u]/(2*X) /* coefficiente di diffusione lungo l'asse z */ a=90 /* durata dell'emissione pari a 90 min*/ Q0=40 /* intensità della sorgente emissiva pari a 40 tons in 90 min */ hmin=100 /* effettiva altezza del punto di emissione delle sostanze inquinanti pari a 39,6 m*/ Hmax=200 /* altezza massima raggiunta dagli agenti inquinanti pari a 200 m*/ ampx=100 /* NON RISPONDENTEampiezza del dominio spaziale lungo l'asse x pari ad 1 Km*/ semiampy=100 /* semiampiezza del dominio spaziale lungo l'asse y pari a 100 m*/ cmax=100 initial values c = 0 equations Ky*DYY(c)+Kz*DZZ(c)-u*DX(c)-Dt(c)=0 extrusion surface 'livello terrestre' Z=0 layer 'spessore del rilievo orografico' surface 'livello superiore del rilievo' Z=80 layer "spessore dell'aria compresa tra sommità del rilievo e sorgente di emissione" surface 'livello della sorgente di emissione' Z=hmin layer 'spessore dell"aria sovrastante alla sorgente di emissione' surface 'estremo superiore del dominio in esame' Z=Hmax boundaries region 1 'intero dominio spaziale' start(0,-semiampy) value(c)=0 line to (ampx,-semiampy) natural(c)=0 line to (ampx,semiampy) value(c)=0 line to (0,semiampy) natural(c)=0 line to close region 2 'regione del rilievo orografico' start(50,-20) line to (90,-20) to (90,20) to (50,20) to close limited region 3 'regione della sorgente di emissione' start(0,-0.2) line to (0.4,-0.2) to (0.4,0.2) to (0,0.2) to close surface 3 value(c)=Q0/u*exp(-x^2-y^2-(z-hmin)^2) time 0 by 10 to a monitors for cycle=1 contour(c) on z=0 as "XY Concentration" range=(0,cmax) contour(c) on x=0 as "YZ Concentration" range=(0,cmax) contour(c) on y=0 as "XZ Concentration" range=(0,cmax) elevation(c) from (0,0,0) to (ampx,0,0) as "X-Axis Concentration" range=(0,cmax) elevation(c) from (0,-semiampy,0) to (0,semiampy,0) as "Y-Axis Concentration" range=(0,cmax) elevation(c) from (0,0,0) to (0,0,Hmax) as "Z-Axis Concentration" range=(0,cmax) plots for t = endtime contour(c) on z=0 as "XY Concentration" range=(0,a) contour(c) on x=0 as "YZ Concentration" range=(0,a) contour(c) on y=0 as "XZ Concentration" range=(0,a) end So the problem shouldn't be the one you pointed out.Have you got any idea about what could be? Thank you very much | ||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 641 Registered: 06-2003 |
You have constructed a rather strange thing here. Your only x-transport term is a convection term U*dx(c). This condition requires that you specify a value of c at the start of each flow path; ie, along the x=0 face. But you have specified natural(c)=0 on this face, which is undefined for first-order terms (see Natural Boundary Conditions in the Help Index). Furthermore, you have constructed a small piece of a gaussian as a value condition on a sheet in x,y at z=hmin. Since this is a sheet at fixed z, it is invisible in the flow direction (x) and provides no information to the flow values. The sheet (delta-function) value condition creates a discontinuous initial value, which FlexPDE attempts to regularize by intense mesh refinement and diffusive smearing. Later when this non-source is smoothed out, the dense mesh is unrefined again. Unfortunately, in this intense refining and unrefining the unrefiner gets confused and fails. This is the diagnostic you see. The way to correct this is to create a well-posed problem. In the attached, I have removed the value condition from the patch at z=hmin and applied it throughout the incoming face (x=0). This provides a starting value for all flow paths entering the system. I have left the patch on surface 3 as a mesh control, and specified a cell size of 0.2/5 to give some resolution to the incoming narrow Gaussian shape. This formulation seems to work, at least as far as t=10, which is as far as I have run it. I had to make several other changes to make it useful: 1) You have specified an initial timestep of 10 minutes, which is wildly optimistic. Because of the very small patch source applied to a sheet at z=100, the solution will require very small timesteps to start, in order to accurately track the steep initial front. By specifying a 10 minute initial stepsize, you force FlexPDE to repeatedly recompute the first step, each time cutting the stepsize until a suitable stepsize is reached. You should either delete the initial stepsize specification (in which case FlexPDE will pick 1e-4 times the end time), or, better, to specify a very small number, realizing that a small step size will be required. This last approach is the most efficient, because no time is wasted re-cutting the initial step. I have used 1e-6. As long as the solution becomes smooth after some time, the stepsize will automatically be increased. 2. Your plot at z=0 will see nothing in the beginning, because the source is at z=100. Your plot at x=0 will see very little, because in this cut the initial sheet source on surface 3 looks like a delta function. 3. You have specified a fixed range for your plots as (0,100), while the source value you impose is only 10. So all the data are contained in the first two contour intervals. 4. Your plots at endtime have specified a data range of (0,a). "a" is the end time. It is unrelated to the data values, and may cause the plots to be meaningless. My modified script follows: title 'Modellizzazione del disastro ecologico di Bhopal' coordinates cartesian3 variables c(threshold=0.1) definitions u=4 { velocità di trasporto del vento pari a 4m s^-1} Ky=[0.0121*X^2*(1+0.0004*X)^(-1)*u]/(2*X) { coefficiente di diffusione lungo l'asse y} Kz=[0.0064*X^2*(1+0.00015*X)^(-1)*u]/(2*X) /* coefficiente di diffusione lungo l'asse z */ a=90 /* durata dell'emissione pari a 90 min*/ Q0=40 /* intensità della sorgente emissiva pari a 40 tons in 90 min */ hmin=100 /* effettiva altezza del punto di emissione delle sostanze inquinanti pari a 39,6 m*/ Hmax=200 /* altezza massima raggiunta dagli agenti inquinanti pari a 200 m*/ ampx=100 /* NON RISPONDENTEampiezza del dominio spaziale lungo l'asse x pari ad 1 Km*/ semiampy=100 /* semiampiezza del dominio spaziale lungo l'asse y pari a 100 m*/ cmax=100 initial values c =Q0/u*exp(-x^2-y^2-(z-hmin)^2) equations Ky*DYY(c)+Kz*DZZ(c)-u*DX(c)-Dt(c)=0 extrusion surface 'livello terrestre' Z=0 layer 'spessore del rilievo orografico' surface 'livello superiore del rilievo' Z=80 layer "spessore dell'aria compresa tra sommità del rilievo e sorgente di emissione" surface 'livello della sorgente di emissione' Z=hmin layer 'spessore dell"aria sovrastante alla sorgente di emissione' surface 'estremo superiore del dominio in esame' Z=Hmax boundaries region 1 'intero dominio spaziale' start(0,-semiampy) value(c)=0 line to (ampx,-semiampy) natural(c)=0 line to (ampx,semiampy) value(c)=0 line to (0,semiampy) value(c)=Q0/u*exp(-x^2-y^2-(z-hmin)^2) line to close region 2 'regione del rilievo orografico' start(50,-20) line to (90,-20) to (90,20) to (50,20) to close limited region 3 'regione della sorgente di emissione' surface 3 mesh_spacing=0.2/5 start(0,-0.2) line to (0.4,-0.2) to (0.4,0.2) to (0,0.2) to close time 0 to a by 1e-6 monitors for cycle=1 grid(x,y) on surface 3 zoom(0,-1,2,2) grid(y,z) on x=0 zoom(-1,hmin-1, 2,2) contour(c) on z=hmin as "XY Concentration" !range=(0,cmax) contour(c) on x=0 as "YZ Concentration" !range=(0,cmax) contour(c) on y=0 as "XZ Concentration" !range=(0,cmax) elevation(c) from (0,0,hmin) to (ampx,0,hmin) as "X-Axis Concentration" !range=(0,cmax) elevation(c) from (0,0,hmin) to (2,0,hmin) as "X-Axis Concentration" !range=(0,cmax) elevation(c) from (0,-semiampy,hmin) to (0,semiampy,hmin) as "Y-Axis Concentration" !range=(0,cmax) elevation(c) from (0,-2,hmin) to (0,2,hmin) as "Y-Axis Concentration" !range=(0,cmax) elevation(c) from (0,0,0) to (0,0,Hmax) as "Z-Axis Concentration" !range=(0,cmax) elevation(c) from (0,0,hmin-2) to (0,0,hmin+2) as "Z-Axis Concentration" !range=(0,cmax) plots for t = endtime contour(c) on z=hmin as "XY Concentration" !range=(0,a) contour(c) on x=0 as "YZ Concentration" !range=(0,a) contour(c) on y=0 as "XZ Concentration" !range=(0,a) end |