Add ascii data class for heating model

[prajnapdas]

I am trying to add a new class for the heating model, where it can read an ascii data (just like for initial temperature). I have tried to make the files I would need to add this class (attached below). However, even though it compiles without an error and even reads the file, it is not actually using the heat sources that I provide in the ascii file.
Could you please let me know if the files I created have some errors and why aspect cannot use it?

P.S.: I do not code in C++, so there might be very obvious problems which I missed.

file: ascii_data.cc

/*
  Copyright (C) 2011 - 2023 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 LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
 */


#include <aspect/global.h>
#include <aspect/heating_model/ascii_data.h>
#include  <aspect/material_model/interface.h>

namespace aspect
{
  namespace HeatingModel
  {
    
    template <int dim>
    AsciiData<dim>::AsciiData () = default;


    template <int dim>
    void
    AsciiData<dim>::initialize ()
    {
      Utilities::AsciiDataInitial<dim>::initialize(1);
    }


    template <int dim>
    double AsciiData<dim>::heating_model (const Point<dim> &position) const
    {
      return Utilities::AsciiDataInitial<dim>::get_data_component(position,0);
    }
    
    template <int dim>
    void
    AsciiData<dim>::
    evaluate( const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
              const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
              HeatingModel::HeatingModelOutputs &heating_model_outputs) const

     {
     std::cout << "evaluating heating model";
     for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
        {
         const Point<dim> position = material_model_inputs.position[q];
         heating_model_outputs.heating_source_terms[q]= heating_model(position) 
		 					* material_model_outputs.densities[q]; 
         heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
	 std::cout << heating_model_outputs.heating_source_terms[q];
         }
     std::cout << std::endl;
     }

    template <int dim>
    void
    AsciiData<dim>::update ()
    {
        return;
    }


    template <int dim>
    void
    AsciiData<dim>::declare_parameters (ParameterHandler &prm)
    {
      prm.enter_subsection ("Heating model");
      {
        Utilities::AsciiDataInitial<dim>::declare_parameters(prm,
							     "$ASPECT_SOURCE_DIR/data/heating-model/ascii-data/test",
                                                             "box_2d.txt");
      }
      prm.leave_subsection();
    }


    template <int dim>
    void
    AsciiData<dim>::parse_parameters (ParameterHandler &prm)
    {
      prm.enter_subsection ("Heating model");
      {
        Utilities::AsciiDataInitial<dim>::parse_parameters(prm);
      }
      prm.leave_subsection();
    }
  }
}

// explicit instantiations
namespace aspect
{
  namespace HeatingModel
  {
    ASPECT_REGISTER_HEATING_MODEL(AsciiData,
                                          "ascii data",
                                          "Implementation of a model in which the initial "
                                          "temperature is derived from files containing data "
                                          "in ascii format. Note the required format of the "
                                          "input data: The first lines may contain any number of comments "
                                          "if they begin with `#', but one of these lines needs to "
                                          "contain the number of grid points in each dimension as "
                                          "for example `# POINTS: 3 3'. "
                                          "The order of the data columns "
                                          "has to be `x', `y', `Heating[W/m^3]' in a 2d model and "
                                          " `x', `y', `z', `Heating[W/m^3]' in a 3d model, which means that "
                                          "there has to be a single column "
                                          "containing the heating. "
                                          "Note that the data in the input "
                                          "files need to be sorted in a specific order: "
                                          "the first coordinate needs to ascend first, "
                                          "followed by the second and the third at last in order to "
                                          "assign the correct data to the prescribed coordinates. "
                                          "If you use a spherical model, "
                                          "then the assumed grid changes. `x' will be replaced by "
                                          "the radial distance of the point to the bottom of the model, "
                                          "`y' by the azimuth angle and `z' by the polar angle measured "
                                          "positive from the north pole. The grid will be assumed to be "
                                          "a latitude-longitude grid. Note that the order "                                              
					                                "of spherical coordinates is `r', `phi', `theta' "
                                          "and not `r', `theta', `phi', since this allows "
                                          "for dimension independent expressions.")
  }
}

file: ascii_data.h

/*
  Copyright (C) 2011 - 2019 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 LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
*/


#ifndef _aspect_heating_model_ascii_data_h
#define _aspect_heating_model_ascii_data_h

#include <aspect/heating_model/interface.h>
#include <aspect/material_model/interface.h>
#include <aspect/utilities.h>

namespace aspect
{
  namespace HeatingModel
  {
    using namespace dealii;

    /**
     * A class that implements a prescribed heating field determined from
     * an AsciiData input file.
     *
     * @ingroup HeatingModels
     */
    template <int dim>
    class AsciiData : public Interface<dim>, public Utilities::AsciiDataInitial<dim>
    {
      public:
        /**
         * Constructor.
         */
        AsciiData ();

        /**
         * Initialization function. Called once at the beginning of the program.
         */
        void initialize () override;
        using Utilities::AsciiDataInitial<dim>::initialize;


	double heating_model (const Point<dim> &position) const; 
        /**
         * Evaluate the heating source term(s) at all quadrature points.
         * This fills outputs.heating_source_terms[q] for each quadrature point.
         */
	void evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
                       const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
                       HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;

        void
        update () override;
        /**
         * Declare the parameters this class takes through input files.
         */
        static void declare_parameters (ParameterHandler &prm);

        /**
         * Read the parameters this class declares from the parameter file.
         */
        void parse_parameters (ParameterHandler &prm) override;
    };
  }
}

#endif

[gassmoeller]

Hi @prajnapdas, that sounds like a reasonable and useful plugin, which we may also be happy to integrate into the main ASPECT version. I took a brief look over your code and it looks alright to me.
In order to help you debug the problem, can you:

• post the input file you use to test your heating plugin, also including the heating ascii data file you used
• also post a screenshot of the paraview output, in particular adding 'heating' to the parameter "List of output variables". Also add "heating statistics" to the "List of postprocessors" parameter and post the log.txt file. Both should show you the output of the heating plugin.
• also describe what heating output you would expect (i.e. if you compute the approximate temperature change from the heating you prescribe, how large is that compared to the temperature you have in your model over the run time of your model

[prajnapdas]

Hi @gassmoeller, thank you so much! I agree that this might be a good addition to the main ASPECT version.
Below, I am attaching the files you asked for. You will have to allow negative heat source to be able to run these.

Input parameter file:

set Dimension                              = 2
set Use years in output instead of seconds = true
set End time                               = 2.0e7
set Use conduction timestep 	         = true
#set Maximum first time step		= 1.0e3
set Pressure normalization 		= surface
set Surface pressure			= 0.0
set Resume computation			= false 
set Output directory                    = output


subsection Material model
  set Model name = simple 

  subsection Simple model
    set Reference density       = 12600 
    set Thermal conductivity    = 110.0
    set Viscosity               = 1.0e19
    set Reference specific heat = 800
    set Thermal expansion coefficient = 1.0e-5
    set Reference temperature   = 0.0
  end
end


subsection Geometry model
  set Model name = sphere

  subsection Sphere
    set Radius  = 1221000
  end
end


subsection Gravity model
  set Model name = radial linear 

  subsection Radial linear 
    set Magnitude at bottom = 0.0 
    set Magnitude at surface = 4.4 
  end
end


subsection Formulation
  set Formulation = Boussinesq approximation 
end


subsection Boundary temperature model
  set Fixed temperature boundary indicators = surface
  set List of model names = spherical constant

  subsection Spherical constant
    set Inner temperature = 0
    set Outer temperature = 0
  end
end


subsection Boundary velocity model
	set Prescribed velocity boundary indicators = surface: zero velocity
end


subsection Initial temperature model
  set Model name = function 

  subsection Function
  	set Variable names = x, y
        set Function constants = dT=-0, c=1221000, Ts=0, a=0.001
        set Function expression = Ts + dT*(1 - (x^2+y^2)/c^2) - a*(((x^2+y^2)^0.5)/c)*(2*(cos(atan2(y,x)))^2-1)
  end
end

subsection Mesh refinement
  set Initial global refinement	       = 5
  set Initial adaptive refinement        = 0
  set Time steps between mesh refinement = 0
end


####################################
####################################
subsection Heating model
   set List of model names = ascii data

   subsection Ascii data model
      set Data directory = /global/home/users/prajnapdas/ASPECT/asciiData_HM/
      set Data file name = heating_model.txt
   end 
end
####################################
####################################


subsection Postprocess
  set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics

  subsection Visualization
    set Time between graphical output = 1.0e6
    set Number of grouped files       = 1
    set List of output variables      = temperature anomaly, stress second invariant, spherical velocity components 
    set Interpolate output = false
  end

end

subsection Checkpointing
  set Time between checkpoint  = 100
end


### higher order polynomial for velocity
subsection Discretization
   set Stokes velocity polynomial degree = 3
end

heating model function as used previously:

subsection Heating model
    set List of model names = function
    subsection Function
        set Variable names = x, y, z
        set Function constants = dTc=145,dk=10,c=1221000.0,rho=12600.0,k=110.0,adT=0,b=3480000
        set Function expression =-(2*dTc*dk/(rho*c))*((x/c)^2 - (y/c)^2)/sqrt(x^2+y^2) - 4*k*adT/(rho*c*c)
    end
end

'heating_model.txt' (which I obtained using the same function):

# inner core heating model, no stratification
# POINTS: 6 6
# Columns: x y Heating[W/m^3]
-1.221000000000000000e+06 -1.221000000000000000e+06 -0.000000000000000000e+00
-7.326000000000000000e+05 -1.221000000000000000e+06 8.472403767896095942e-14
-2.442000000000000000e+05 -1.221000000000000000e+06 1.453284656175965674e-13
2.442000000000000000e+05 -1.221000000000000000e+06 1.453284656175965674e-13
7.326000000000000000e+05 -1.221000000000000000e+06 8.472403767896095942e-14
1.221000000000000000e+06 -1.221000000000000000e+06 -0.000000000000000000e+00
-1.221000000000000000e+06 -7.326000000000000000e+05 -8.472403767896095942e-14
-7.326000000000000000e+05 -7.326000000000000000e+05 -0.000000000000000000e+00
-2.442000000000000000e+05 -7.326000000000000000e+05 7.811170319825381910e-14
2.442000000000000000e+05 -7.326000000000000000e+05 7.811170319825381910e-14
7.326000000000000000e+05 -7.326000000000000000e+05 -0.000000000000000000e+00
1.221000000000000000e+06 -7.326000000000000000e+05 -8.472403767896095942e-14
-1.221000000000000000e+06 -2.442000000000000000e+05 -1.453284656175965674e-13
-7.326000000000000000e+05 -2.442000000000000000e+05 -7.811170319825381910e-14
-2.442000000000000000e+05 -2.442000000000000000e+05 -0.000000000000000000e+00
2.442000000000000000e+05 -2.442000000000000000e+05 -0.000000000000000000e+00
7.326000000000000000e+05 -2.442000000000000000e+05 -7.811170319825381910e-14
1.221000000000000000e+06 -2.442000000000000000e+05 -1.453284656175965674e-13
-1.221000000000000000e+06 2.442000000000000000e+05 -1.453284656175965674e-13
-7.326000000000000000e+05 2.442000000000000000e+05 -7.811170319825381910e-14
-2.442000000000000000e+05 2.442000000000000000e+05 -0.000000000000000000e+00
2.442000000000000000e+05 2.442000000000000000e+05 -0.000000000000000000e+00
7.326000000000000000e+05 2.442000000000000000e+05 -7.811170319825381910e-14
1.221000000000000000e+06 2.442000000000000000e+05 -1.453284656175965674e-13
-1.221000000000000000e+06 7.326000000000000000e+05 -8.472403767896095942e-14
-7.326000000000000000e+05 7.326000000000000000e+05 -0.000000000000000000e+00
-2.442000000000000000e+05 7.326000000000000000e+05 7.811170319825381910e-14
2.442000000000000000e+05 7.326000000000000000e+05 7.811170319825381910e-14
7.326000000000000000e+05 7.326000000000000000e+05 -0.000000000000000000e+00
1.221000000000000000e+06 7.326000000000000000e+05 -8.472403767896095942e-14
-1.221000000000000000e+06 1.221000000000000000e+06 -0.000000000000000000e+00
-7.326000000000000000e+05 1.221000000000000000e+06 8.472403767896095942e-14
-2.442000000000000000e+05 1.221000000000000000e+06 1.453284656175965674e-13
2.442000000000000000e+05 1.221000000000000000e+06 1.453284656175965674e-13
7.326000000000000000e+05 1.221000000000000000e+06 8.472403767896095942e-14
1.221000000000000000e+06 1.221000000000000000e+06 -0.000000000000000000e+00

'log.txt' file for ascii_data heating model (few initial lines):

-----------------------------------------------------------------------------
-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
--     . version 2.5.0
--     . using deal.II 9.5.1
--     .       with 32 bit indices and vectorization level 1 (128 bits)
--     . using Trilinos 14.4.0
--     . using p4est 2.3.2
--     . running in DEBUG mode
--     . running with 48 MPI processes
-----------------------------------------------------------------------------


   Loading Ascii data initial file /global/home/users/prajnapdas/ASPECT/asciiData_HM/heating_model.txt.

-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
--   https://aspect.geodynamics.org/citing.html?ver=2.5.0&sha=&src=code
-----------------------------------------------------------------------------
Number of active cells: 5,120 (on 6 levels)
Number of degrees of freedom: 133,764 (92,546+20,609+20,609)

*** Timestep 0:  t=0 years, dt=0 years
   Skipping temperature solve because RHS is zero.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 97+0 iterations.

   Postprocessing:
     Writing graphical output:           output/solution/solution-00000
     RMS, max velocity:                  5.57e-08 m/year, 1.01e-06 m/year
     Temperature min/avg/max:            0 K, 0 K, 0 K
     Heat fluxes through boundary parts: -75.27 W
     Heating rate (average/total):       -9.78e-14 W/kg, -5771 W

*** Timestep 1:  t=1.07246e+06 years, dt=1.07246e+06 years
   Solving temperature system... 37 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 92+0 iterations.

   Postprocessing:
     Writing graphical output:           output/solution/solution-00001
     RMS, max velocity:                  5.46e-08 m/year, 9.97e-07 m/year
     Temperature min/avg/max:            -0.005871 K, -0.00394 K, 0 K
     Heat fluxes through boundary parts: -344 W
     Heating rate (average/total):       -9.78e-14 W/kg, -5771 W

expected temperature anomaly (obtained using the function for the heating model):

temperature anomaly when using ascii_data for heating model:

[gassmoeller]

Ah, I think I see the problem, and it is not with the code, it is an assumption in your data file. If you use the ascii data plugin in a model with spherical geometry, it:

• automatically converts the position into spherical coordinates
• assumes that the data file you provide is also stored in terms of spherical coordinates.

Because the data file you use is using Cartesian coordinates, it defines the heating term in an area far outside your model domain.

Here is an example for an ASCII data file for spherical models (for initial temperature):

# Test data for ascii data initial conditions.
# Only next line is parsed in format: [nx] [ny] [nz] because of keyword "POINTS:"
# POINTS: 3 3 3
# Columns: r phi      temperature [K]
3481000    0          000.
4908500    0          000.
6336000    0          000.
3481000    0.7854     000.
4908500    0.7854     100.
6336000    0.7854     000.
3481000    1.5708     000.
4908500    1.5708     000.
6336000    1.5708     000.

Can you rework your data file into spherical coordinates and check your test case again?