| Author | Message | ||||
     Mohammad Rahmani (mrahmani) Member Username: mrahmani Post Number: 8 Registered: 10-2004 |
Suppose a set of two first order PDEs as following: f: dx(f)+dy(f)=(x+y)*(x*y-g) g: dx(g)+dy(g)=(x+y)*(1-f) subjected to the following BCs: @ x=0, f=0, g=-1 @ y=0, f=0, g=-1 where x: [0,1], y: [0,1] We know the exact solution of the above set is: fex=sinh(x*y) !Exact Solution of f gex=-cosh(x*y)+x*y !Exact solution of g Describing this set in FlexPDE 4.2.3s gives wrong solution. The problem is with BCs. How we can treat such problems in FlexPDE? I should note that if we use the exact solution and define all the four BCs, we will get the acurate solution by FlexPDE but practically we have only two BCs. On the other boundaries FlexPDE assume a natural BCs which is not correct. Thanks /Mohammad ---------------------------------------------------- The following code describe the above Problem TITLE 'Test of 1st order PDEs in two dimensions ' { This is only a test to ses how FlexPDE treats a set of coupled first order PDEs in two dimensins. we have the exact solution as f=sinh(xy) g=-cosh(xy)+xy } VARIABLES { system variables } f g DEFINITIONS fex=sinh(x*y) !Exact Solution of f gex=-cosh(x*y)+x*y !Exact solution of g EQUATIONS { PDE's, one for each variable } f: dx(f)+dy(f)=(x+y)*(x*y-g) g: dx(g)+dy(g)=(x+y)*(1-f) BOUNDARIES REGION 1 START(0,0) value(f)=fex value(g)=gex line to (0,1) nobc(f) nobc(g) line to (1,1) nobc(f) nobc(g) line to (1,0) value(f)=fex value(g)=gex line to finish MONITORS contour(f) contour(g) PLOTS contour(f) as "Calculated f " contour(g) as "Calculated g" contour(fex) as "Exact f" contour(gex) as "Exact g" contour(g-gex) as "Error in g" contour(f-fex) as "Error in f" elevation(f-fex) from (0,0.5) to (1,0.5) as "Error in f" elevation(g-gex) from (0,0.5) to (1,0.5) as "Error in g " END | ||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 262 Registered: 06-2003 |
I don't understand your remarks. I ran this problem as you sent it and got solutions with maximum errors of 1e-5. Is this what you are calling the "wrong solution"? Given that the default error tolerance in FlexPDE is 1e-3, this is a better answer than you had a right to expect. In regard to the Natural BC, FlexPDE integrates only second-order terms by parts, which is the operation that gives rise to the surface terms defined by the Natural BC. In your system, there are no second-order terms, so the Natural is undefined. Using NoBC, as you have done, is the correct treatment. This applies no conditions on the affected boundaries. First-order equations require the specification of a single boundary condition along each characteristic path. In your system, the characteristics are given by dx/dy=1, so the BC's you have specified are exactly the right ones. And, as I have said, the correct solution results. Where's the beef? | ||||
| Mohammad Rahmani (mrahmani) Member Username: mrahmani Post Number: 9 Registered: 10-2004 |
Dear Robert, I have ran the above code and got a very different solutions. Please see attached the exact solution and the calculated one by FlexPDE. Don't you think there is a big difference between exact values of f and g and their calculated values? ********************************************** Results, using FlexPDE 4.2.3s       | ||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 263 Registered: 06-2003 |
You didn't tell me you were using the student version, so I tested it on the professional, which uses slightly different default settings. There seems to be some kind of strange resonance between the equations and the mesh size. If you run with NGRID=11, it solves nicely. See attached. Interestingly, it seems to work with any ODD gridsize over 10. Probably selecting higher density along the bottom and left boundaries would work, too.
| ||||
| Mohammad Rahmani (mrahmani) Member Username: mrahmani Post Number: 10 Registered: 10-2004 |
Dear Robert, Thanks for your guidelines. I ran the code with NGRID=11 and got the right solutions but I also got a very good results for NGRID=1,2,3,4,5. The error is in the range of 1e-3 or smaller. I implemented another 1st order PDE set from NAG toolbox in FlexPDE 4.2.4s and encountered the same problem i.e there is no right solution for NGRID=10 or 9 but very good results for NGRID=1, 5 ,7. What's wrong here? Is there any other hidden settings in student version which cause this behaviour. /Mohammad | ||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 266 Registered: 06-2003 |
Pure convection equations when discretized have typically very sparse coupling, and can sometimes admit oscillatory solutions. The classic (if oversimplified) example is dx(u)+ k*dxx(u)=0 with value conditions at both ends. The linear finite element discretization of this problem is oscillatory as k becomes small, finally becoming insoluble for k=0 and odd number of cells. While quadratic triangles with nonzero source terms are much more complex than this simple example, clearly some similar process is going on in small meshes. FlexPDE implements an upwind weighting shift, which is supposed to cure this problem. This essentially adds a diffusion term in the direction of flow to shift the differencing to the upwind side. It is apparently insufficient in some grid combinations for your problem. You might try increasing the upwind weight: SELECT UPFACTOR=4 This is effective for your original problem with NGRID=9. |
