LKM-code-base / RotatingMHD

GNU General Public License v3.0
1 stars 0 forks source link

[master] Convergence tests do not produce the expected results #114

Closed j507 closed 3 years ago

j507 commented 3 years ago

I wrote a python script to automatically run all convergence tests. While troubleshooting I kept getting erroneous results for the hydrodynamic problems. Turns out that problem is in the code and not in the script as I ran the .sh files manually obtaining for Guermond.cc

                               Velocity convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                H1               Linfty       
    0 1.65e-02  1024 8450 4.42e-02 6.320056e-04     - 1.096780e-02    - 1.100207e-03     - 
    0 2.11e-02  1024 8450 4.42e-02 7.055821e-04 -0.45 1.041052e-02 0.21 1.216645e-03 -0.41 

                               Pressure convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    0 1.65e-02  1024 1089 4.42e-02 3.085360e-04     - 1.289136e-02     - 2.038607e-03     - 
    0 2.11e-02  1024 1089 4.42e-02 3.330317e-04 -0.31 1.300297e-02 -0.04 2.872494e-03 -1.40 

and for TGV.cc

                               Velocity convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 7.99e-02  1024 8450 4.42e-02 1.200102e-02     - 1.492662e-01     - 1.846496e-02     - 
    5 8.80e-02  1024 8450 4.42e-02 1.284660e-02 -0.70 1.580094e-01 -0.59 1.954694e-02 -0.59 

                               Pressure convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 7.99e-02  1024 1089 4.42e-02 3.837631e-03     - 9.512578e-02     - 1.240798e-02     - 
    5 8.80e-02  1024 1089 4.42e-02 4.013308e-03 -0.46 9.668196e-02 -0.17 1.280077e-02 -0.32 

The problem seems to be that the time step is not being read. As a reference the results of the christensen_benchmark branch are

                               Velocity convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 1.00e-01  1024 8450 4.42e-02 2.597065e-02     - 3.127355e-01     - 3.945044e-02     - 
    5 5.00e-02  1024 8450 4.42e-02 7.677519e-03 -1.76 8.534220e-02 -1.87 1.093958e-02 -1.85 

                               Pressure convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 1.00e-01  1024 1089 4.42e-02 8.936978e-03     - 1.493329e-01     - 2.731630e-02     - 
    5 5.00e-02  1024 1089 4.42e-02 2.565349e-03 -1.80 8.426195e-02 -0.83 8.218220e-03 -1.73 

and

                               Velocity convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 1.00e-01  1024 8450 4.42e-02 8.700332e-03     - 1.407325e-01     - 2.771303e-02     - 
    5 5.00e-02  1024 8450 4.42e-02 2.105930e-03 -2.05 2.201024e-02 -2.68 2.659283e-03 -3.38 

                               Pressure convergence table
==============================================================================================
level    dt    cells dofs   hmax           L2                 H1               Linfty       
    5 1.00e-01  1024 1089 4.42e-02 7.336106e-03     - 6.912819e-02     - 3.570945e-02     - 
    5 5.00e-02  1024 1089 4.42e-02 1.889239e-03 -1.96 1.760381e-02 -1.97 1.044779e-02 -1.77 

respectively.

As for ThermalTGV.cc the results of master are

| level | dt       | cells | dofs  | hmax     | L2           | L2...red.rate.log2 | H1           | H1...red.rate.log2 | Linfty       | Linfty...red.rate.log2 |
| 8     | 1.00e-01 | 65536 | 66049 | 5.52e-03 | 7.253380e-04 | -                  | 1.573704e-02 | -                  | 1.497214e-03 | -                      |
| 8     | 5.00e-02 | 65536 | 66049 | 5.52e-03 | 1.685061e-04 | -2.11              | 1.438242e-02 | -0.13              | 3.838579e-04 | -1.96                  |
| 8     | 2.50e-02 | 65536 | 66049 | 5.52e-03 | 3.540994e-05 | -2.25              | 1.429836e-02 | -0.01              | 1.171757e-04 | -1.71                  |
| 8     | 1.25e-02 | 65536 | 66049 | 5.52e-03 | 7.805197e-06 | -2.18              | 1.429303e-02 | -0.00              | 5.626794e-05 | -1.06                  |
| 8     | 6.25e-03 | 65536 | 66049 | 5.52e-03 | 1.130275e-05 | 0.53               | 1.429255e-02 | -0.00              | 3.726674e-05 | -0.59                  |
| 8     | 3.12e-03 | 65536 | 66049 | 5.52e-03 | 1.108751e-05 | -0.03              | 1.429248e-02 | -0.00              | 2.817476e-05 | -0.40                  |
| 8     | 1.56e-03 | 65536 | 66049 | 5.52e-03 | 5.806826e-05 | 2.39               | 1.430603e-02 | 0.00               | 1.892189e-04 | 2.75                   |
| 8     | 7.81e-04 | 65536 | 66049 | 5.52e-03 | 2.878814e-04 | 2.31               | 1.454336e-02 | 0.02               | 5.660994e-04 | 1.58                   |

and those of christensen_benchmark are

                               Temperature convergence table
==============================================================================================
level    dt    cells dofs    hmax           L2                 H1               Linfty       
    8 1.00e-01 65536 66049 5.52e-03 7.253564e-04     - 1.573710e-02     - 1.496996e-03     - 
    8 5.00e-02 65536 66049 5.52e-03 1.685344e-04 -2.11 1.438244e-02 -0.13 3.842268e-04 -1.96 
    8 2.50e-02 65536 66049 5.52e-03 3.532344e-05 -2.25 1.429834e-02 -0.01 1.165773e-04 -1.72 
    8 1.25e-02 65536 66049 5.52e-03 7.556621e-06 -2.22 1.429297e-02 -0.00 5.138846e-05 -1.18 
    8 6.25e-03 65536 66049 5.52e-03 1.223268e-05  0.69 1.429256e-02 -0.00 3.186282e-05 -0.69 
    8 3.12e-03 65536 66049 5.52e-03 1.106205e-05 -0.15 1.429249e-02 -0.00 2.793086e-05 -0.19 

The results are similar but I don't see why there should be a difference at all. Furthermore the L2 norm behaves weirdly. There is a positive convergence rate all of the sudden in both instances. The error field from the 6.25e-2 cycle of the christensen_benchmark branch looks like this:

Screenshot from 2021-05-18 14-23-10

while the rest looks like this:

Screenshot from 2021-05-18 14-27-12

The solution itself still looks physically sound.

To spare time I had until now only checked the first two cycle when comparing between branches so I missed this entirely. I will have to look where it started because the original convergence tests (Those seen in the presentations) all look good with the expected convergence order.

sebglane commented 3 years ago

Could you send me the input file for TGV.cc which produces the convergence tables we discussed?

sebglane commented 3 years ago

Useful commands:

  1. List all branches merged into master
    git branch --remotes --merged master
  2. Get hash of most recent commit contained in master and other_branch
    git rev-list --max-count=1 master other_branch
j507 commented 3 years ago

Could you send me the input file for TGV.cc which produces the convergence tables we discussed?

All temporal convergence tests are to be done with 8 cycles and with a spatial discretization of 8 global refinements. For TGV.cc:

# Listing of Parameters
# ---------------------
set FE's polynomial degree - Pressure (Taylor-Hood) = 1
set Mapping - Apply to interior cells               = false
set Mapping - Polynomial degree                     = 1
set Problem type                                    = hydrodynamic
set Spatial dimension                               = 2
set Verbose                                         = false

subsection Convergence test parameters
  set Convergence test type                 = temporal
  set Number of spatial convergence cycles  = 2
  set Number of temporal convergence cycles = 8
  set Time-step reduction factor            = 0.5
end

subsection Dimensionless numbers
  set Reynolds number         = 100
end

subsection Navier-Stokes solver parameters
  set Convective term time discretization    = semi-implicit
  set Convective term weak form              = skew-symmetric
  set Incremental pressure-correction scheme = rotational
  set Preconditioner update frequency        = 1
  set Verbose                                = false

  subsection Linear solver parameters - Correction step
    set Absolute tolerance           = 1e-9
    set Maximum number of iterations = 1000
    set Relative tolerance           = 1e-6

    subsection Preconditioner parameters
      set Preconditioner type      = ILU
      set Fill-in level            = 2
      set Absolute tolerance       = 1e-5
      set Relative tolerance       = 1.01
      set Overlap                  = 2
    end

  end

  subsection Linear solver parameters - Diffusion step
    set Absolute tolerance           = 1e-9
    set Maximum number of iterations = 1000
    set Relative tolerance           = 1e-6

    subsection Preconditioner parameters
      set Preconditioner type      = ILU
      set Fill-in level            = 2
      set Absolute tolerance       = 1e-5
      set Relative tolerance       = 1.01
      set Overlap                  = 2
    end

  end

  subsection Linear solver parameters - Poisson pre-step
    set Absolute tolerance           = 1e-9
    set Maximum number of iterations = 1000
    set Relative tolerance           = 1e-6

    subsection Preconditioner parameters
      set Preconditioner type      = ILU
      set Fill-in level            = 2
      set Absolute tolerance       = 1e-5
      set Relative tolerance       = 1.01
      set Overlap                  = 2
    end

  end

  subsection Linear solver parameters - Projection step
    set Absolute tolerance           = 1e-9
    set Maximum number of iterations = 1000
    set Relative tolerance           = 1e-6

    subsection Preconditioner parameters
      set Preconditioner type      = ILU
      set Fill-in level            = 2
      set Absolute tolerance       = 1e-5
      set Relative tolerance       = 1.01
      set Overlap                  = 2
    end

  end

end

subsection Output control parameters
  set Graphical output directory = TGVResults/
  set Graphical output frequency = 10000
  set Terminal output frequency  = 1
end

subsection Refinement control parameters
  set Adaptive mesh refinement               = false
  set Maximum number of levels               = 10
  set Number of initial global refinements   = 8
end

subsection Time stepping parameters
  set Adaptive time stepping        = false
  set Adaptive timestepping barrier = 2
  set Final time                    = 1.0
  set Initial time step             = 0.1
  set Maximum number of time steps  = 10
  set Maximum time step             = 2e-1
  set Minimum time step             = 9e-2
  set Start time                    = 0.0
  set Time stepping scheme          = BDF2
  set Verbose                       = false
end
j507 commented 3 years ago
2\. `git rev-list --max-count=1 master other_branch`

We ought go through all the pull requests and restore any deleted branches. At least until mit_benchmark.

sebglane commented 3 years ago

We ought go through all the pull requests and restore any deleted branches. At least until mit_benchmark.

I will take care of that.

sebglane commented 3 years ago

Could you send me the python script for running the convergence test?

j507 commented 3 years ago

Could you send me the python script for running the convergence test?

Still working in the script. Here is the "old" version (For some reason I keep failing the authentication needed to push)

import argparse     # Module to parse arguments
import subprocess   # Module to execute terminal commands
import enum         # Module to enable a enum class
import fileinput    # Module to modify .txt files
import numpy

class Test_type(enum.Enum):
  """Test_type Enum class for the type of convergence test to be done
  """
  spatial   = 1
  temporal  = 2
  both      = 3

class Dimensionles_number(enum.Enum):
  """Dimensionles_number Enum class for the relevant dimensionless
  numbers
  """
  Reynolds  = 1
  Peclet    = 2

class ConvergenceTest:
  """Class encompasing all the data related to the convergence test.

  The class stores the convergence tests parameters passed on to the
  constructor and the pertinent line numbers using the `get_line_number`
  function

  :param source_file_name: The name of the source file
  :type source_file_name: str
  :param dimensionless_number_values: A numpy array containing the values
    of the dimensionless number for which a convergence test is
    to be performed, defaults to numpy.array([100])
  :type dimensionless_number_values: numpy.array, optional
  :param dimensionless_number_type: The pertinent dimensionless number
    of the problem (Reynolds or Peclet), defaults to
    Dimensionles_number.Reynolds
  :type dimensionless_number_type: class.Dimensionless_number, optional
  :param convergence_test_type: The type of convergence test to be done,
    defaults to Test_type.temporal
  :type convergence_test_type: class.Test_type, optional
  :param n_spatial_cycles: The number of cycles to be performed in the
    spatial convergence test, defaults to 2
  :type n_spatial_cycles: int, optional
  :param n_initial_global_refinements: The initial number of global
    refinements, defaults to 5
  :type n_initial_global_refinements: int, optional
  :param time_step: The fixed time step size used throughout the spatial
    convergence test, defaults to 1e-3
  :type time_step: float, optional
  :param n_temporal_cycles: The number of cycles to be performed in the
    temporal convergence test, defaults to 2
  :type n_temporal_cycles: int, optional
  :param n_global_refinements: The fixed number of global refinements
    used throughout the temporal convergence test, defaults to 8
  :type n_global_refinements: int, optional
  :param initial_time_step: The initial time step size, defaults to 1e-1
  :type initial_time_step: float, optional
  :param exec_file: The shell script which actually runs the convergence
    test.
  :type exec_file: str
  :param source_file: The source file
  :type source_file: str
  :param prm_file: The .prm file
  :type prm_file: str
  :param line_dimensionless_number: The line number of the .prm file
    at which the pertinent dimensionless number is set
  :type line_dimensionless_number: int
  :param line_convergence_test_type: The line number of the .prm file
    at which the convergence test type is set
  :type line_convergence_test_type: int
  :param line_n_spatial_cycles: The line number of the .prm file
    at which number of spatial cycles is set
  :type line_n_spatial_cycles: int
  :param line_n_initial_global_refinements: The line number of the .prm
    file at which the initial number of global refinements is set
  :type line_n_initial_global_refinements: int
  :param line_n_temporal_cycles: The line number of the .prm file
    at which the number of temporal cycles is set
  :type line_n_temporal_cycles: int
  :param line_initial_time_step: The line number of the .prm file
    at which the initial time step size is set
  :type line_initial_time_step: int
  :param line_graphical_output_directory: The line number of the .prm
    file at which graphical output directory is set
  :type line_graphical_output_directory: int
  :param line_graphical_output_frequency: The line number of the .prm
    file at which graphical output frequency is set
  :type line_graphical_output_frequency: int

  """
  def __init__(self,
               source_file_name,
               dimensionless_number_values  = numpy.array([100]),
               dimensionless_number_type    = Dimensionles_number.Reynolds,
               convergence_test_type        = Test_type.temporal,
               n_spatial_cycles             = 2,
               n_initial_global_refinements = 5,
               time_step                    = 1e-3,
               n_temporal_cycles            = 2,
               n_global_refinements         = 8,
               initial_time_step            = 1e-1):
    """Contructor method

    Stores the arguments passed on to the constructor in the class'
    attributes. Determines the line numbers of the relevant data inside
    the .prm file using the `get_line_number` function and stores them
    in the class' attributes.

    """
    self.source_file_name             = source_file_name
    self.dimensionless_number_values  = dimensionless_number_values
    self.dimensionless_number_type    = dimensionless_number_type
    self.convergence_test_type        = convergence_test_type
    self.n_spatial_cycles             = n_spatial_cycles
    self.n_initial_global_refinements = n_initial_global_refinements
    self.time_step                    = time_step
    self.n_temporal_cycles            = n_temporal_cycles
    self.initial_time_step            = initial_time_step
    self.n_global_refinements         = n_global_refinements
    self.exec_file                    = "./" + source_file_name + ".sh"
    self.source_file                  = source_file_name + ".cc"
    self.prm_file                     = "applications/" + source_file_name + ".prm"
    self.line_dimensionless_number          = self.get_line_number(
                                              "subsection Dimensionless numbers") + 1
    self.line_convergence_test_type         = self.get_line_number(
                                              "Convergence test type")
    self.line_n_spatial_cycles              = self.get_line_number(
                                              "Number of spatial convergence cycles")
    self.line_n_initial_global_refinements  = self.get_line_number(
                                              "Number of initial global refinements")
    self.line_n_temporal_cycles             = self.get_line_number(
                                              "Number of temporal convergence cycles")
    self.line_initial_time_step             = self.get_line_number(
                                              "Initial time step")
    self.line_graphical_output_directory    = self.get_line_number(
                                              "Graphical output directory")
    self.line_graphical_output_frequency    = self.get_line_number(
                                              "Graphical output frequency")
    self.replace_line(self.line_graphical_output_frequency,
                      "  set Graphical output frequency = 10000")

  def get_line_number(self, string_to_search):
    """get_line_number Locates the `string_to_search` in the `prm_file`
    and returns its line number

    :param string_to_search: String to search in the `prm_file`
    :type string_to_search: str
    :return: Line number where `string_to_search` is found
    :rtype: int

    """
    line_number = 0
    with open(self.prm_file, 'r') as read_obj:
      for line in read_obj:
        line_number += 1
        if string_to_search in line:
          break
    return line_number

  def replace_line(self, line_number, string_to_replace):
    """replace_line Replaces the `line_number`-th line of the `prm_file`
    with `string_to_replace`

    :param line_number: The number of the line to be overwritten
    :type line_number: int
    :param string_to_replace: String which overwrites the
      `line_number`-th line
    :type string_to_replace: str

    """
    for line in fileinput.input(self.prm_file, inplace=True):
      if fileinput.filelineno() == line_number:
        print(string_to_replace)
      else:
        print(line, end="")

def nproc_type(x):
  """nproc_type A function acting as a type for admissible number of
  processors, i.e. only positive integers bigger than zero

  :param x: Desired number of processors
  :type x: int
  :raises argparse.ArgumentTypeError: Error if `x` is smaller than 1
  :return: Admissible number of proccesors
  :rtype: int

  """
  x = int(x)
  if x < 1:
    raise argparse.ArgumentTypeError("Minimum number of processors is 1")
  return x

def main(nproc):
  branch  = subprocess.run(["git", "rev-parse", "--abbrev-ref", "HEAD"],
                           stdout=subprocess.PIPE,
                           universal_newlines=True)
  hash    = subprocess.run(["git", "rev-parse", "--short", "HEAD"],
                           stdout=subprocess.PIPE,
                           universal_newlines=True)
  print("Convergence tests script running with " + str(nproc)
        + " processors")
  print(" Branch: " + str(branch.stdout).strip("\n"))
  print(" Hash:   " + str(hash.stdout))

  # Load all the convergence tests and its setting into a list
  tests = [ConvergenceTest("ThermalTGV",
                           numpy.array([100]),
                           Dimensionles_number.Peclet,
                           Test_type.temporal,
                           2,
                           4,
                           1e-3,
                           2,
                           5,
                           0.1),
           ConvergenceTest("TGV",
                           numpy.array([100]),
                           Dimensionles_number.Reynolds,
                           Test_type.temporal,
                           2,
                           4,
                           1e-3,
                           2,
                           5,
                           0.1)]

  # Loop over the list
  for test in tests:
    print("  Source file: " + test.source_file)

    # Set the tests' number of cycles in the .prm file
    test.replace_line(test.line_n_spatial_cycles,
                      "  set Number of spatial convergence cycles  = %s"
                        % test.n_spatial_cycles)
    test.replace_line(test.line_n_temporal_cycles,
                      "  set Number of temporal convergence cycles = %s"
                        % test.n_temporal_cycles)

    # Loops over the values of the dimensionless number
    for dimensionless_number_value in test.dimensionless_number_values:
      # Set the test's dimensionless value in the .prm file
      if test.dimensionless_number_type == Dimensionles_number.Reynolds:
        test.replace_line(test.line_dimensionless_number,
                          "  set Reynolds number         = %s"
                            % dimensionless_number_value)
        print("    Reynolds number: " + str(dimensionless_number_value))
      elif test.dimensionless_number_type == Dimensionles_number.Peclet:
        test.replace_line(test.line_dimensionless_number,
                          "  set Peclet number           = %s"
                            % dimensionless_number_value)
        print("    Peclet number: " + str(dimensionless_number_value))
      else:
        raise Exception("Invalid enum value")

      # if-branches depending on the type of convergence test. Only the
      # first one is commented as the all follow the same pattern
      if test.convergence_test_type == Test_type.spatial:
        # Terminal output
        print("      Running:")
        print("       Convergence test type                 = spatial")
        print("       Number of cycles                      = "
                      + str(test.n_spatial_cycles))
        print("       Number of initial global refinements  = "
                      + str(test.n_initial_global_refinements))
        print("       Time step                             = "
                      + str(test.time_step))

        # Set the pertinent values in the .prm file
        test.replace_line(
              test.line_convergence_test_type,
              "  set Convergence test type                 = spatial")
        test.replace_line(
              test.line_n_initial_global_refinements,
              "  set Number of initial global refinements   = %s"
                % test.n_initial_global_refinements)
        test.replace_line(
              test.line_initial_time_step,
              "  set Initial time step             = %s"
                % test.time_step)
        test.replace_line(
              test.line_graphical_output_directory,
              "   set Graphical output directory = %s_Spatial_%s/"
                % (test.source_file_name, dimensionless_number_value))

        # Run the shell script
        with open("Dump_%s_Spatial_%s.txt"
                    % (test.source_file_name,
                       dimensionless_number_value),
                  "w") as dump_file:
          process = subprocess.run([test.exec_file, str(nproc)],
                                  stdout=dump_file,
                                  stderr=subprocess.STDOUT)

      elif test.convergence_test_type == Test_type.temporal:
        print("      Running:")
        print("       Convergence test type                 = temporal")
        print("       Number of cycles                      = "
                      + str(test.n_temporal_cycles))
        print("       Number of global refinements          = "
                      + str(test.n_global_refinements))
        print("       Initial time step                     = "
                      + str(test.initial_time_step))

        test.replace_line(
              test.line_convergence_test_type,
              "  set Convergence test type                 = temporal")
        test.replace_line(
              test.line_n_initial_global_refinements,
              "  set Number of initial global refinements   = %s"
                % test.n_global_refinements)
        test.replace_line(
              test.line_initial_time_step,
              "  set Initial time step             = %s"
                % test.initial_time_step)
        test.replace_line(
              test.line_graphical_output_directory,
              "  set Graphical output directory = %s_Temporal_%s/"
                % (test.source_file_name, dimensionless_number_value))

        with open("Dump_%s_Temporal_%s.txt"
                    % (test.source_file_name,
                       dimensionless_number_value),
                  "w") as dump_file:
          process = subprocess.run([test.exec_file, str(nproc)],
                                  stdout=dump_file,
                                  stderr=subprocess.STDOUT)

      elif test.convergence_test_type == Test_type.both:
        print("      Running:")
        print("       Convergence test type                 = spatial")
        print("       Number of cycles                      = "
                      + str(test.n_spatial_cycles))
        print("       Number of initial global refinements  = "
                      + str(test.n_initial_global_refinements))
        print("       Time step                             = "
                      + str(test.time_step))

        test.replace_line(
              test.line_convergence_test_type,
              "  set Convergence test type                 = spatial")
        test.replace_line(
              test.line_n_initial_global_refinements,
              "  set Number of initial global refinements   = %s"
                % test.n_initial_global_refinements)
        test.replace_line(
              test.line_initial_time_step,
              "  set Initial time step             = %s"
                % test.time_step)
        test.replace_line(
              test.line_graphical_output_directory,
              "  set Graphical output directory = %s_Spatial_%s/"
                % (test.source_file_name, dimensionless_number_value))

        with open("Dump_%s_Spatial_%s.txt"
                    % (test.source_file_name,
                       dimensionless_number_value),
                  "w") as dump_file:
          process = subprocess.run([test.exec_file, str(nproc)],
                                  stdout=dump_file,
                                  stderr=subprocess.STDOUT)

        print("      Running:")
        print("       Convergence test type                 = temporal")
        print("       Number of cycles                      = "
                      + str(test.n_temporal_cycles))
        print("       Number of global refinements          = "
                      + str(test.n_global_refinements))
        print("       Initial time step                     = "
                      + str(test.initial_time_step))

        test.replace_line(
                test.line_convergence_test_type,
                "  set Convergence test type                 = temporal")
        test.replace_line(
              test.line_n_initial_global_refinements,
              "  set Number of initial global refinements   = %s"
                % test.n_global_refinements)
        test.replace_line(
              test.line_initial_time_step,
              "  set Initial time step             = %s"
                % test.initial_time_step)
        test.replace_line(
                test.line_graphical_output_directory,
                "  set Graphical output directory = %s_Temporal_%s/"
                  % (test.source_file_name, dimensionless_number_value))

        with open("Dump_%s_Temporal_%s.txt"
                    % (test.source_file_name,
                       dimensionless_number_value),
                  "w") as dump_file:
          process = subprocess.run([test.exec_file, str(nproc)],
                                  stdout=dump_file,
                                  stderr=subprocess.STDOUT)

      else:
        raise Exception("Invalid enum value")
    print("")
  print("Tests completed!")
  print(" The convergence tables are in the applications/*.txt files")
  print(" The terminal output are in the Dump_*.txt files")

if __name__ == "__main__":
  parser = argparse.ArgumentParser()
  parser.add_argument("--nproc",
                      "-np",
                      type = nproc_type,
                      default = 4,
                      help = "Number of processors")
  args = parser.parse_args()
  main(args.nproc)

This only runs for the latest .prm file format and for the branches that have the shell scripts to run the applications.

j507 commented 3 years ago

Update @sebglane : The results have indicated that the problematic branches were neumann_boundary_conditions, reintroduce_tolerance_assert and first_step_fix (Pull-Request #47, #68, #70 respectively). The problems introduced by #68 are easily corrected by choosing an absolute tolerance of 1e-11 and relative tolerance of 1e-10 on all solvers. Through this fix the expected convergence rates of the Guermond.cc and TGV.cc appear to be recovered.

The convergence rates of ThermalTGV are nonetheless still plagued by an increase of the L2 error norm after a certain time step size. The other norms behave as expected. As discussed, this could be either a problem inside ThermalTGV.cc, problem_class.* or in boundary_conditions.*. I would not discard the possibility that there is an error inside the convection diffusion solver itself as we observed an increase of the Linfty norm in the mit_benchmark branch as well.

As reference here are the resutls of mit_benchmark

                               Temperature convergence table
==============================================================================================
level    dt    cells dofs    hmax           L2                 H1               Linfty       
    8 1.00e-01 65536 66049 5.52e-03 3.422179e-04     - 1.459552e-02     - 6.388700e-04     - 
    8 5.00e-02 65536 66049 5.52e-03 9.724730e-05 -1.82 1.431300e-02 -0.03 1.485290e-04 -2.10 
    8 2.50e-02 65536 66049 5.52e-03 3.221060e-05 -1.59 1.429359e-02 -0.00 1.732616e-05 -3.10 
    8 1.25e-02 65536 66049 5.52e-03 1.626297e-05 -0.99 1.429249e-02 -0.00 1.651510e-05 -0.07 
    8 6.25e-03 65536 66049 5.52e-03 1.262790e-05 -0.36 1.429246e-02 -0.00 2.488773e-05  0.59 
    8 3.12e-03 65536 66049 5.52e-03 1.178418e-05 -0.10 1.429247e-02  0.00 2.694933e-05  0.11 
    8 1.56e-03 65536 66049 5.52e-03 1.157977e-05 -0.03 1.429247e-02  0.00 2.747482e-05  0.03 
    8 7.81e-04 65536 66049 5.52e-03 1.154990e-05 -0.00 1.429247e-02  0.00 2.754016e-05  0.00 

and here the results of christensen_benchmark with corrected tolerances

                               Temperature convergence table
==============================================================================================
level    dt    cells dofs    hmax           L2                 H1               Linfty
    8 1.00e-01 65536 66049 5.52e-03 7.254068e-04     - 1.573729e-02     - 1.496475e-03     -
    8 5.00e-02 65536 66049 5.52e-03 1.683148e-04 -2.11 1.438222e-02 -0.13 3.819421e-04 -1.97
    8 2.50e-02 65536 66049 5.52e-03 3.537623e-05 -2.25 1.429835e-02 -0.01 1.148627e-04 -1.73
    8 1.25e-02 65536 66049 5.52e-03 7.438373e-06 -2.25 1.429294e-02 -0.00 4.926018e-05 -1.22
    8 6.25e-03 65536 66049 5.52e-03 9.572181e-06  0.36 1.429253e-02 -0.00 3.298919e-05 -0.58
    8 3.12e-03 65536 66049 5.52e-03 1.101154e-05  0.20 1.429248e-02 -0.00 2.893446e-05 -0.19
    8 1.56e-03 65536 66049 5.52e-03 1.139844e-05  0.05 1.429248e-02 -0.00 2.792520e-05 -0.05
    8 7.81e-04 65536 66049 5.52e-03 1.149028e-05  0.01 1.429247e-02 -0.00 2.768662e-05 -0.01

Inside the tubcloud folder linked above there is a folder with the results of this last table. The odd frame numbers correspond to the graphical output at the end of each cycle. Interestingly enough, the error scalar field computed in ThermalTGV.cc seems to to monotonously decrease, contradicting the results seen in the L2 norm. Therefore it might be worth checking if the norms are being computed correctly.

As mentioned, I will push the tolerances in christensen_benchmark and leave a Pull-Request to correct them in the current 'master` branch.

sebglane commented 3 years ago

There is only a problem remaining in the thermal Taylor-Green benchmark.