| Author | Message | |||||
     Kyle (kyle) Member Username: kyle Post Number: 5 Registered: 08-2006 |
Okay, I've tried everything I can think of on this script, but I cannot get past the startup error "timestep has fallen below 1e-9." I do have Value BC's (which are bad), but I made a surface-fit for the initial condition in an effort to match the BC's, and then I even used SWAGE. I also added extra dummy regions to force local mesh refinement near the Value BC's, but this didn't help either. I think this should have a smooth startup given: (a) BC matching initial condition, (b) use of SWAGE, and (c) local mesh refinement, but No Luck!!! I give up - I'm ready for help!
| |||||
| Kyle (kyle) Member Username: kyle Post Number: 7 Registered: 08-2006 |
Due to trouble uploading attachment, here is the FlexPDE script in text form. TITLE 'Kyles Problem' { the problem identification } COORDINATES cartesian2 { coordinate system, 1D,2D,3D, etc } SELECT ERRLIM = 1e-3 CUBIC = ON SMOOTHINIT = ON VARIABLES { system variables } c (THRESHOLD = 1e-7) ! [=] (grams H2O/cm^3 "of Region 1"). Region 1 is Nitrogen Gas @ 333.15 Kelvin and 1 atm. SCALAR VARIABLES Flux_into_Solid2 (THRESHOLD = 1e-7) ! [g/(cm-depth)*sec] Mass_into_Solid2 (THRESHOLD = 1e-7) ! [g/(cm-depth)] DEFINITIONS { parameter definitions } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. ! Problem Geometry X1 = 0.00250 ! [=] cm. Distance to Gas1/Solid1 interface. X2 = X1 + 0.0100 ! [=] cm. Distance to Solid1/Gas2 interface. X3 = X2 + 0.00250 ! [=] cm. Distance to Gas2/Solid2 interface. X4 = X3 + 0.01000 ! [=] cm. Distance to axis of symmetry. Y1 = 0.0025 ! [=] cm. Distance to mid-height Solid1. Y2 = 0.0035 ! [=] cm. Distance to Solid2/Gas2 interface. Y3 = 0.00425 ! [=] cm. Distance to mid-height Gas2 ("Head-space" - main region of interest). y4 = 0.0050 ! [=] cm. Distance to top of cavity. Solid2_volume = (X4 - X3)*(Y2) ! [=] [(cm^3)/(cm-depth)]. Solid2_density = 1.50 ! [=] [g/(cm^3)]. Solid2_mass = (Solid2_volume)*(Solid2_density) ! [=] [g/(cm-depth)]. Initial_loading = 0.001 ! Dimensionless mass fraction. Solid2_loading = Initial_loading + ((Mass_into_Solid2)/(Solid2_mass)) ! Dimensionless mass fraction. p_over_Solid2 = (5.90180e-4)*(exp((47.41955)*(Solid2_loading))) ! [=] torr. ! Note: Initial value for p_over_Solid2 = 6.1884e-4 ! [=] torr. c_over_Solid2 = ((p_over_Solid2)*(18.015))/((6.2363e4)*(333.15)) ! [=] [g/(cm^3 Gas @ 60C)]. ! Note: Initial value for c_over_Solid2 = 5.36594e-10 ! [=] [g/(cm^3 Gas @ 60C)]. c_source = 1.1662e-4 ! [=] (g/cm^3 Gas @ 60C). aa = 1.166199999999750e-4 bb = 2.986215050484080e+8 cc = -9.810590733903220e+7 dd = 1.549303502201380e+18 ee = -8.412442749016930e+11 z_initial = (aa + bb*x + cc*y)/(1 + dd*x + ee*y) INITIAL VALUES c = z_initial ! [=] g/cm^3. Initial moisture concentration. EQUATIONS { PDE's, one for each variable } c: div(D*grad(K*c)) = dt(K*c) Flux_into_Solid2: Flux_into_Solid2 = (line_integral(-D*normal(grad(K*c)),"Witness_boundary_Solid2")) ! [g/(cm-depth)*sec)] Mass_into_Solid2: dt(Mass_into_Solid2) = (Flux_into_Solid2) ! [g/(cm-depth)] BOUNDARIES { The domain definition } REGION 1 { Gas1 and Gas2 at 333.15 Kelvin and 1 atm } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. START(0, 0) { Walk the domain boundary } NATURAL(c) = 0 LINE TO (X4, 0) NATURAL(c) = 0 LINE TO (X4, Y4) NATURAL(c) = 0 LINE TO (0, Y4) VALUE(c) = SWAGE(t, z_initial, c_source, 60) LINE TO FINISH REGION 2 { Solid1 } D = 7.383e-9 ! [=] (cm^2/sec). K = 215.06 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. START(X1, 0) { Walk the domain boundary } LINE TO (X2, 0) TO (X2, Y4) TO (X1, Y4) TO FINISH REGION 3 { Solid2 } D = 0.606 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. START(X3, 0) { Walk the domain boundary } LINE TO (X4, 0) TO (X4, Y2) VALUE(c) = SWAGE(t, z_initial, c_over_Solid2, 60) LINE TO (X3, Y2) VALUE(c) = SWAGE(t, z_initial, c_over_Solid2, 60) LINE TO FINISH REGION 4 { Added for mesh control only } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. mesh_spacing = (1e-4) START(0, 0) { Walk the domain boundary } LINE TO (0.0001, 0) TO (0.0001, Y4) TO (0, Y4) TO FINISH REGION 5 { Added for mesh control only } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. mesh_spacing = (1e-4) START((X1 - 0.0001), 0) { Walk the domain boundary } LINE TO (X1, 0) TO (X1, Y4) TO ((X1 - 0.0001), Y4) TO FINISH REGION 6 { Added for mesh control only } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. mesh_spacing = (1e-4) START((X2 + 0.0001), 0) { Walk the domain boundary } LINE TO (X2, 0) TO (X2, Y4) TO ((X2 + 0.0001), Y4) TO FINISH REGION 7 { Added for mesh control only } D = 0.303 ! [=] (cm^2/sec). K = 1 ! [=] [(grams H2O/cm^3 "of present Region")/(grams H2O/cm^3 "of Region 1")]. mesh_spacing = (1e-4) START((X3 - 0.0001), 0) { Walk the domain boundary } LINE TO (X3, 0) TO (X3, Y2) TO (X4, Y2) TO (X4, (Y2 + 0.0001)) TO ((X3 - 0.0001), (Y2 + 0.0001)) TO FINISH FEATURE START "Witness_line_X0" (0, 0) LINE TO (0, Y4) START "Witness_line_X1" (X1, 0) LINE TO (X1, Y4) START "Witness_line_X2" (X2, 0) LINE TO (X2, Y4) START "Witness_line_Y1" (0, Y1) LINE TO (X4, Y1) START "Witness_line_Y3" (0, Y3) LINE TO (X4, Y3) START "Witness_boundary_Solid2" (X3, 0) LINE TO (X3, Y2) LINE TO (X4, Y2) TIME { if time dependent } FROM 0 TO 1080000 MONITORS { show progress } For t = 0 BY 1 to 30 by 30 to 3600 by 600 TO 1080000 grid(x,y) ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" HISTORY(Flux_into_Solid2) AS "[g/(cm-depth)*sec]" HISTORY(Mass_into_Solid2) AS "[g/(cm-depth)]" HISTORY(Solid2_loading) AS "Dimensionless mass fraction" HISTORY(p_over_Solid2) AS "[torr]" HISTORY(c_over_Solid2) AS "[g/(cm^3 Gas @ 60C)]" PLOTS { save result displays } For t = 0 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 1 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 30 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 3600 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 72000 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 144000 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 288000 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" For t = 1080000 ELEVATION(c) ON "Witness_line_Y1" ELEVATION(c) ON "Witness_line_Y3" ELEVATION(K*c) ON "Witness_line_Y1" ELEVATION(K*c) ON "Witness_line_Y3" HISTORIES HISTORY(Flux_into_Solid2) AS "[g/(cm-depth)*sec]" HISTORY(Mass_into_Solid2) AS "[g/(cm-depth)]" HISTORY(Solid2_loading) AS "Dimensionless mass fraction" HISTORY(p_over_Solid2) AS "[torr]" HISTORY(c_over_Solid2) AS "[g/(cm^3 Gas @ 60C)]" END | |||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 689 Registered: 06-2003 |
I ran this problem to time 600 with both version 5.0.11 and 5.0.12 on Windows with no sign of trouble. What version and platform are you using? | |||||
| Kyle (kyle) Member Username: kyle Post Number: 8 Registered: 08-2006 |
That's actually great news! I am running FlexPDE Professional Version 4.2.16 2D on Windows XP. What extra startup features does version 5.0.12 have that allows it to run this script? Best regards, Kyle | |||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 690 Registered: 06-2003 |
This appears to be related to the discontinuous initial conditions. Version 5 does an initial diffusion pass to smear discontinuities in initial conditions. Using version 4, if you specify an initial timestep of 1e-5, it will cut back to about 1e-12 or less and grind a while, but it will eventually (after some 300 timesteps) get the timestep up and continue at a reasonable rate. All this implies that version 4 is struggling with the discontinuous initial values. If you can come up with a way to smear these out a bit, it should help. I tried just ramping the c boundary condition in time, but this isn't sufficient, because the Zinitial has some fine structure to it that causes the same problem. Your Swage doesn't do it because the swage centers the smear about the specified switchpoint. So if you swage about the start time, you still have a jump of half the target value at t=0. Maybe a ramp together with a bigger Threshold? The version 5 smearing seems to work for this problem, but if the structure of zinitial is important (which it must not be, or you wouldn't be trying for a large timestep) you should start with a small step anyway and turn off the smear. | |||||
| Kyle (kyle) Member Username: kyle Post Number: 9 Registered: 08-2006 |
I am definitely not interested in the fine structure of the initial condition. This was just my attempt to smooth out the startup from a "zero" initial condition to VALUE BCs on a 2D geometry. I am open to other forms of an initial condition that might satisfy the startup. I still do not understand why SWAGE doesn't smear out the start sufficiently. I actually thought this was indeed the purpose of this function. If the command is SWAGE(t, initial_condition, VALUE_BC, 60) aren’t the SWAGE function and its derivative supposed to be smooth and continuous from zero to 60 seconds with no steps? Also, in Version 4, how can you specify an initial time step without activation of FIXDT? Thank you, Kyle |
