← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/modules/convdiff/ProblemConvDiff.C

  • Committer: logg
  • Date: 2002-09-13 12:55:37 UTC
  • Revision ID: devnull@localhost-20020913125537-gz6ry1id9xsvu6np
Tailorized "2002-09-13 07:55:37 by logg"
Initial revision

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/modules/convdiff/ProblemConvDiff.C
 
1
// Copyright (C) 2002 Johan Hoffman and Anders Logg.
 
2
// Licensed under the GNU GPL Version 2.
 
3
 
 
4
#include "ProblemConvDiff.hh"
 
5
#include "EquationConvDiff_cG1dG0.hh"
 
6
#include "EquationConvDiff2d_cG1dG0.hh"
 
7
#include <stdlib.h>
 
8
 
 
9
//-----------------------------------------------------------------------------
 
10
const char *ProblemConvDiff::Description()
 
11
{
 
12
  if ( space_dimension == 2 )
 
13
         return "Convection-Diffusion equation (2D)";
 
14
  else
 
15
         return "Convection-Diffusion equation (3D)";
 
16
}
 
17
//-----------------------------------------------------------------------------
 
18
void ProblemConvDiff::Solve()
 
19
{
 
20
  Equation *equation;
 
21
  SparseMatrix A;
 
22
  Vector b;
 
23
  Vector u(no_nodes);
 
24
  Vector up(no_nodes);
 
25
  KrylovSolver solver;
 
26
  real T0,T;
 
27
  
 
28
  // Choose equation
 
29
  if ( space_dimension == 2 )
 
30
         equation = new EquationConvDiff2d_cG1dG0();
 
31
  else
 
32
         equation = new EquationConvDiff_cG1dG0();
 
33
  
 
34
  // Global fields
 
35
  GlobalField Up(grid,&up);
 
36
  GlobalField U(grid,&u);
 
37
  GlobalField f(grid,"source");
 
38
  GlobalField eps(grid,"diffusivity");
 
39
  GlobalField bx(grid,"x-convection");
 
40
  GlobalField by(grid,"y-convection");
 
41
  GlobalField bz(grid,"z-convection");
 
42
  
 
43
  // Attach fields to equation
 
44
  equation->AttachField(0,&Up);
 
45
  equation->AttachField(1,&eps);
 
46
  equation->AttachField(2,&f);
 
47
  equation->AttachField(3,&bx);
 
48
  equation->AttachField(4,&by);
 
49
  if ( space_dimension == 3 )
 
50
         equation->AttachField(5,&bz);
 
51
  
 
52
  // Set up the discretiser
 
53
  Discretiser discretiser(grid,equation);
 
54
 
 
55
  // Prepare time stepping
 
56
  settings->Get("start time",&T0);
 
57
  settings->Get("final time",&T);
 
58
  int time_step = 0;
 
59
  real t        = T0;
 
60
  real dt       = grid->GetSmallestDiameter();
 
61
  
 
62
  // Save initial value
 
63
  U.SetLabel("u","Temperature");
 
64
  U.Save(t);
 
65
  
 
66
  // Start time stepping
 
67
  while (t < T){
 
68
 
 
69
    time_step++;
 
70
         t += dt;
 
71
    
 
72
    // Set solution to previous solution
 
73
    up.CopyFrom(&u);
 
74
 
 
75
         // Write info
 
76
    display->Progress(0,(t-T0)/(T-T0),"i=%i t=%f",time_step,t);
 
77
    display->Message(0,"||u0||= %f ||u1|| = %f",up.Norm(),u.Norm());
 
78
 
 
79
         // Set time and time-step
 
80
         equation->SetTime(t);
 
81
         equation->SetTimeStep(dt);
 
82
         
 
83
    // Discretise equation
 
84
    discretiser.Assemble(&A,&b);
 
85
         
 
86
    // Solve the linear system
 
87
    solver.Solve(&A,&u,&b);
 
88
         
 
89
    // Save the solution
 
90
         U.Save(t);
 
91
  }
 
92
 
 
93
  // Clean up
 
94
  delete equation;
 
95
}
 
96
//-----------------------------------------------------------------------------