DetailedModel.m

classdef DetailedModel < Greedy.User.FVDetailedModelDefault
  % Detailed model for a linear evolution problem with a @ref fv discretization
  %
  % This is a detailed model for a problem of the type
  % @f{align}
  %   \partial_t u(x;t,\mu) - {\cal L}(t,\mu)[u(x,t;\mu)] &= 0 \qquad &\text{in }& \Omega \times [0, T], \\\\
  %   u(x;0,\mu) &= u_0(x;\mu) \qquad &\text{in }& \Omega \times \{0\}
  % @f}
  %  with suitable boundary conditions and linear operators
  %  `{\cal L}(t,\mu) : {\cal W} \to {\cal W}` parametrized over the
  %  time-parameter domain `[0,T] \times {\cal M}`.
  %
  % The numerical scheme for the numerical solution of this problem is assumed
  % to be of the following type: 
  % @f{align*}
  % u_h(\cdot;t^0,\mu) &= {\cal M}_h[u_0] \\\\
  % \left(\text{Id} - \Delta t^k {\cal L}_{h, I}(t^{k+1},\mu)\right)[u_h(\cdot;t^{k+1},\mu)]
  %    &=
  % \left(\text{Id} + \Delta t^k {\cal L}_{h, E}(t^{k},\mu)\right)[u_h(\cdot;t^{k},\mu)]
  % @f}
  % with the following ingredients: 
  % - a time discretization `0 = t^0 \leq t^1 = t^0 + \Delta t^0 \leq \ldots
  %   \leq t^{K+1} = t^{K} + \Delta t^{K} = T`,
  % - solution snapshots `u_h(\cdot;t,\mu) \in { cal W }_h ` for all
  % parameters and time instances `(t,\mu) \in [0,T] \times { \cal M }`
  % - separable operators for the discretization `{\cal L}_{h, I}(t,\mu):{\cal
  % W}_h \to {\cal W}_h`, `{\cal L}_{h, E}(t,\mu):{\cal W}_h \to {\cal W}_h`
  % and separable initial data `u_0` in the parameter argument.  That is these
  % must be decomposable like
  % @f{align*}
  % {\cal L}_{h, *}(t,\mu) &= \sum_{q=1}^{q_*} \sigma_*^q(\mu) {\cal L}^q_{h, *}(t) \\\\
  % u_0(x;\mu)   &= \sum_{q=1}^{u_0} \sigma_{u_0}^q(\mu) u_0(x)
  % @f} and
  % - a projection operator `{\cal M}_h : {\cal W} \to {\cal W}_h`


  properties(Constant)
    % This constant is used for a consistency check of the model descr with
    % help of Greedy.Remove.IDetailedModelBase.struct_check()
    ddescr_checks = struct(...
      'mass_matrix',           {{@(x) isequal(class(x), 'function_handle')}}...
      );
  end

  methods
    function dm = DetailedModel(descr)
      % function dm = DetailedModel(descr)
      % constructor based on problem description
      %
      % Parameters:
      %  descr: structure describing the problem and the descretization
      dm = dm@Greedy.User.FVDetailedModelDefault(descr);
      if ~IDetailedModel.struct_check(descr, ...
                                      LinEvol.DetailedModel.ddescr_checks)
        error('Consistency check for LinEvol.DetailedModel failed.');
      end
    end

    function this         = set_mu(this, mu)
      % function this = set_mu(this, mu, propagate)
      % sets the active parameter vector `\mu \in {\cal M}`
      %
      % The parameter set here, is used by the detailed_simulation() function.
      %
      % The default implementation sets all fieldnames 'descr.mu_names' in the
      % #descr struct.
      %
      % Parameters:
      %   mu: The parameter vector `\boldsymbol\mu`.
      %
      % Return values:
      %  this: handle of type IDetailedModel to the changed DetailedModel

      assert(length(mu) == length(this.descr.mu_names));

      this.descr.mexptr('set_mu', mu);
    end

    function mu           = get_mu(this)
      % function mu = get_mu(this)
      % returns the active parameter vector `\boldsymbol\mu \in {\cal M}`
      %
      % The default implementation returns a vector of the values of the fields
      % of the #descr structure with names 'descr.mu_names'.
      %
      % Return values:
      %   mu: The parameter vector `\boldsymbol\mu`

      mu = this.descr.mexptr('get_mu');
    end
    function sim_data   = detailed_simulation(this, model_data)
      sim_data = this.descr.mexptr('detailed_simulation');
    end

    function model_data = gen_model_data(this)
      model_data.id     = this.descr.mexptr('gen_model_data');
      model_data.mexptr = this.descr.mexptr;
    end

    function p          = plot_sim_data(this, model_data, sim_data, plot_params)

      if ~isfield(plot_params, 'gcf') || isempty(plot_params.gcf)
        ud.descr = this.descr;
        ud.dmodel = this;
        ud.sim_data = sim_data;
        fig_handle = default_plot_control;
        set(fig_handle, 'UserData', ud);
      else
        fig_handle = plot_params.gcf;
      end
      fig_ud     = get(fig_handle, 'UserData');
      if ~isempty(fig_ud) && isfield(fig_ud, 'plot_dir')
        plot_dir = fig_ud.plot_dir;
      else
        args = ['-c ', sim_data.cell_array_name, ' ', sim_data.outdir, ' ', sim_data.prefix];
        [status, res] = system([fullfile(rbmatlabhome, 'dune', 'time_scene_vis.sh'), ' ', args]);
        plot_dir = strtrim(res);
        set(fig_handle, 'CloseRequestFcn', @LinEvolDune.DetailedModel.close_figure_and_vtk_plot);
        fig_ud.plot_dir = plot_dir;
        set(fig_handle, 'UserData', fig_ud);
        pause(1);
      end

      if ~isfield(plot_params, 'timestep')
        plot_params.timestep = 0;
      end

      if isfield(plot_params, 'reload_required') && plot_params.reload_required
        LinEvolDune.DetailedModel.issue_plot_command(plot_dir, ['reload ', num2str(plot_params.timestep)]);
      else
        LinEvolDune.DetailedModel.issue_plot_command(plot_dir, ['time ', num2str(plot_params.timestep)]);
      end
      p = fig_handle;
    end
  end

  methods (Static)

    function U = get_dofs_from_sim_data(sim_data)

      U = sim_data.a;
    end

    function U = get_dofs_at_time(sim_data, time)
      U = sim_data.a;
    end
  end

  methods (Static, Hidden=true)

    function issue_plot_command(plot_dir, command)
      f = fopen(fullfile(plot_dir, 'control.cfg'), 'w');
      fprintf(f, command);
      fclose(f);
      system(['kill -SIGINT $(cat ', fullfile(plot_dir, 'tsv.pid'), ')']);
    end

    function close_figure_and_vtk_plot(src, event)
      fig_ud = get(gcf, 'UserData');
      if ~isempty(fig_ud) && isfield(fig_ud, 'plot_dir')
        plot_dir = fig_ud.plot_dir;
        try
          LinEvolDune.DetailedModel.issue_plot_command(plot_dir, 'exit');
          pause(0.5);
        catch me
          disp('Could not issue the exit command to the Paraview renderer');
          disp(me.message);
        end
        if exist(fullfile(plot_dir, 'tsv.pid'), 'file')
          system(['rm -rf ', plot_dir]);
        end
      end
      closereq
    end
  end

end