[aspect-devel] VectorTools::integrate_difference fails to return an error that decreases under grid refinement (Elbridge Gerry Puckett)

Harsha Lokavarapu lokavarapuh at gmail.com
Tue Oct 25 18:01:47 PDT 2016


Hi,

After further testing, we have discovered that for the exact same BCs 
detailed in the previous post but for a 2D box, [-0.5,0.5] x [-0.5,0.5], 
and a different forcing term as specified in attached parameter file (as 
opposed to 0 in annulus); integrate_difference converges beautifully. 
The L2 norm of the velocity converges at 3rd order for Q2_Q1 elements as 
expected for a 2D box but not for a 2D shell geometry model.

The call to VectorTools::integrate_difference, exactly as in SolKz and 
SolCx returns

                                                     TABLE 01

      Refinement
         Level       |  L1 Velocity      |    Rate   |   L2 Velocity     
|   Rate
-----------------------------------------------------------------------------------------------------
           2           |  6.273785e-05  |              | 6.213717e-05 |
           3           |  6.581690e-06  |    3.25   |  8.542105e-06 |   
2.863
           4           |  5.474571e-07  |    3.59   |  1.080996e-06 |   
2.982
           5           |  4.264332e-08  |    3.68   |  1.355643e-07 |   
2.995
           6           |  3.192472e-09  |    3.74   |  1.695992e-08 |   
2.999

Can the 2D shell model be used in conjunction with the function 
VectorTools::integrate_difference 
<https://dealii.org/8.4.0/doxygen/deal.II/namespaceVectorTools.html#a01174a2a7e2ee8fa6abdfdd93ac7a317> 
to measure L1 and L2 norms?

Thanks,
Harsha

On Tuesday 25 October 2016 01:18 PM, 
aspect-devel-request at geodynamics.org wrote:
> Send Aspect-devel mailing list submissions to
> 	aspect-devel at geodynamics.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> 	http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> or, via email, send a message with subject or body 'help' to
> 	aspect-devel-request at geodynamics.org
>
> You can reach the person managing the list at
> 	aspect-devel-owner at geodynamics.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of Aspect-devel digest..."
>
>
> Today's Topics:
>
>     1. VectorTools::integrate_difference fails to return an error
>        that decreases under grid refinement (Elbridge Gerry Puckett)
>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20161025/0443df83/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: egp_hvl_box.cc
Type: text/x-c++src
Size: 22543 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20161025/0443df83/attachment-0001.cc>
-------------- next part --------------
set Additional shared libraries            = ./libegp_hvl_box.so

############### Global parameters

set Dimension                              = 2

set Start time                             = 0

set End time                               = 0

set Output directory                       = output_egp_hvl_box_2

#set Pressure normalization                 = volume
set Nonlinear solver scheme                = Stokes only

set Use years in output instead of seconds     = false

set CFL number                             = 0.8

set Linear solver tolerance                    = 1e-14

#set Resume computation = true
############### Parameters describing the model

subsection Geometry model
  set Model name = box

  subsection Box
    set X extent = 1
    set X repetitions = 1
    set Y extent = 1

    set Box origin X coordinate = -0.5
    set Box origin Y coordinate = -0.5
  end
end


subsection Model settings

  set Prescribed velocity boundary indicators = 0: function, 1: function, 2: function, 3: function
  set Tangential velocity boundary indicators = 
  set Zero velocity boundary indicators       =
end

subsection Boundary velocity model
  subsection Function
    set Variable names      = x,y
    set Function constants  = pi=6.283185307179586
    set Function expression = -y;x
  end
end

subsection Material model
  set Model name = EGPHVLMaterial

  subsection EGPHVL
    # Do not change these set of parameters! EGP && HL
    set Viscosity jump = 1
    set Reference density = 1
  end
end


subsection Gravity model
  set Model name = function
  subsection Radial constant
    set Magnitude = 0
  end

  subsection Function
    set Variable names      = x,y
    set Function constants = PI=3.141592653589793
    set Function expression = (1/(sqrt(x*x+y*y)))*x;(1/(sqrt(x*x+y*y)))*y
  end
end


############### Parameters describing the temperature field

subsection Boundary temperature model
  set Model name = box
end


subsection Initial conditions
  set Model name = function
  subsection Function
    set Function expression = 0
  end
end

subsection Adiabatic conditions model
  set Model name = function

  subsection Function
    set Function constants  =
    set Function expression = 0; 0
    set Variable names      = depth
  end

end


############### Parameters describing the discretization

subsection Discretization
  set Stokes velocity polynomial degree       = 2
  set Use locally conservative discretization = false
end


subsection Mesh refinement
  set Initial adaptive refinement              = 0
  set Initial global refinement = 6
  set Time steps between mesh refinement       = 0
  set Coarsening fraction                      = 0
end


############### Parameters describing what to do with the solution

subsection Postprocess
  set List of postprocessors = visualization, EGPHVLPostprocessor
  
  subsection Visualization
    set Output format                 = vtu
    set Number of grouped files       = 1
    set Time between graphical output = 0
#    set List of output variables = density, viscosity, partition, gravity
  end
end

subsection Termination criteria
  set Checkpoint on termination = true
  set Termination criteria      = user request

  subsection User request
    set File name = end
  end
end

subsection Checkpointing
  #set Steps between checkpoint = 5
end

-------------- next part --------------
# Copyright (C) 2015 by the authors of the ASPECT code.
#
# This file is part of ASPECT.
#
# ASPECT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2, or (at your option)
# any later version.
#
# ASPECT is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ASPECT; see the file doc/COPYING.  If not see
# <http://www.gnu.org/licenses/>.

CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)

FIND_PACKAGE(Aspect 1.5.0 QUIET HINTS ${Aspect_DIR} ../ ../../ $ENV{ASPECT_DIR})

IF (NOT Aspect_FOUND)
  MESSAGE(FATAL_ERROR "\n"
	"Could not find a valid ASPECT build/installation directory. "
	"Please specify the directory where you are building ASPECT by passing\n"
	"   -D Aspect_DIR=<path to ASPECT>\n"
	"to cmake or by setting the environment variable ASPECT_DIR in your shell "
	"before calling cmake. See the section 'How to write a plugin' in the "
        "manual for more information.")
ENDIF ()

DEAL_II_INITIALIZE_CACHED_VARIABLES()

SET(TARGET1 "egp_hvl_box")
#SET(TARGET2 "egp_hvl_annulus")

PROJECT(${TARGET1})

ADD_LIBRARY(${TARGET1} SHARED egp_hvl_box.cc)
#ADD_LIBRARY(${TARGET2} SHARED egp_hvl_annulus.cc)

IF(DEBUG)
  set_target_properties(${TARGET1} PROPERTIES LINK_FLAGS "-DDEBUG")
#  set_target_properties(${TARGET2} PROPERTIES LINK_FLAGS "-DDEBUG")
ENDIF(DEBUG)

ASPECT_SETUP_PLUGIN(${TARGET1})
#ASPECT_SETUP_PLUGIN(${TARGET2})


More information about the Aspect-devel mailing list