[aspect-devel] Problems with periodic boundary conditions

[Max Rudolph]

Ian,
Thanks for your help. Another way of dealing with the null space
associated with a net translation (or net rotation in 3D spherical
geometry) would be to fix the horizontal velocity for one vx degree of
freedom (essentially a dirichlet boundary condition for one dof). I
suppose that it would be key to do this in such a way as to keep the
stiffness matrix symmetric. Then, calculate the net translation or
linear momentum in the x-direction and subtract it after solving the
momentum and continuity equations. Would this solve the problem?

Max

On Thu, Nov 6, 2014 at 3:23 PM, Ian Rose <ian.rose at berkeley.edu> wrote:
> Hi Max,
>
> I suspect that you are running afoul of the nullspace in your problem.  I
> played with the parameter file a bit, and the solver failures seem to be
> very dependent on mesh resolution, refinements, and domain geometry (though
> not Ra directly, as far as I can tell).
>
> The periodic box with free slip boundaries has the compatibility condition
>
> \int_V \rho g_x = 0
>
> Of course, a vertical gravity field should satisfy this automatically, but
> it might be the case that roundoff errors could cause problems.  I still
> suspect the nullspace because I can make the crashes go a way by (a)
> reducing the solver tolerance or (b) making one of the boundaries no slip.
> What do others think, am I completely off base?  If this is the problem,
> what would be the best way to address it?
>
> Ian
>
>
> On Wed, Nov 5, 2014 at 8:58 PM, Max Rudolph <maxwellr at gmail.com> wrote:
>>
>> Dear All,
>> I'm trying to set up a simple 2D calculation (in a box) with periodic
>> boundary conditions. I keep running into the error message below
>> (about numerical breakdown). A solution is obtained for the first
>> timestep (0) but for the second timestep the 'numerical breakdown'
>> error appears. It seems that this error is related to Ra - aspect runs
>> fine for the low-Ra case but when I increase the layer depth so that
>> Ra~10^8, there are errors. I tried re-meshing every timestep and
>> reducing the CFL number to 0.01 to no avail. Is this a known problem?
>> Either way, do you have any suggestions about how to work around this?
>>
>> Thanks!
>> Max
>>
>> ----------------------------------------------------
>> Exception on processing:
>>
>> --------------------------------------------------------
>> An error occurred in line <246> of file
>> </opt/deal.II-8.1-src/source/lac/trilinos_solver.cc> in function
>>     void dealii::TrilinosWrappers::SolverBase::execute_solve(const
>> dealii::TrilinosWrappers::PreconditionBase &)
>> The violated condition was:
>>     false
>> The name and call sequence of the exception was:
>>     ExcMessage("AztecOO::Iterate error code -2: " "numerical breakdown")
>> Additional Information:
>> AztecOO::Iterate error code -2: numerical breakdown
>> --------------------------------------------------------
>>
>> Aborting!
>> ----------------------------------------------------
>>
>> input file used (left file is my modified input file, right file is
>> cookbook.
>>
>> max at harding:~/aspect_runs$ diff periodic_box.prm
>> /opt/aspect-dev/cookbooks/future/periodic_box.prm
>> 2c2
>> < set CFL number                             = 0.01
>> ---
>> > set CFL number                             = 1.0
>> 4c4
>> < set Output directory                       = output_pb
>> ---
>> > set Output directory                       = output
>> 39,41c39,41
>> <     set X extent = 4.2e6
>> <     set Y extent = 3.0e6
>> < #    set Z extent = 5.e5
>> ---
>> >     set X extent = 1.e6
>> >     set Y extent = 5.e5
>> >     set Z extent = 5.e5
>> 58,60c58
>> <     set Function constants = Ttop=0.0, Tbottom=1000.0, p=0.01,
>> H=3.0E6, W=4.2e6, pi=3.1415926
>> <     set Function expression = Tbottom - (Tbottom-Ttop)*y/H +
>> p*sin(pi*y/H)*cos(2.0*pi*x/W)
>> <
>> ---
>> >     set Function expression = if((sqrt((x-1.e5)^2+(y-4.0e5)^2)<5.0e4) |
>> > (sqrt((x-3.e5)^2+(y-2.e5)^2)<1.0e5) , 800.0, 0)
>> 85c83
>> <   set Time steps between mesh refinement = 1                       #
>> default: 10
>> ---
>> >   set Time steps between mesh refinement = 10                       #
>> > default: 10
>> 100c98
>> <   set List of postprocessors = velocity statistics, visualization,
>> temperature statistics, heat flux statistics, basic statistics
>> ---
>> >   set List of postprocessors = visualization
>> 104,105c102,103
>> <     set Output format                 = hdf5
>> <     set Time between graphical output = 0.0
>> ---
>> >     set Output format                 = vtu
>> >     set Time between graphical output = 1.e3
>> max at harding:~/aspect_runs$
>> _______________________________________________
>> Aspect-devel mailing list
>> Aspect-devel at geodynamics.org
>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
>