~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/modules/navierstokes/ProblemNS.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

files added:

.hgignore

AUTHORS

COPYING

ChangeLog

INSTALL

Makefile

Makefile.am

Makefile.in

NEWS

README

aclocal.m4

configure

configure.in

data

data/grids

data/grids/dolfin-1.inp

data/grids/dolfin-2.inp

data/grids/tetgrid_16_16_16.inp

data/grids/tetgrid_1_1_1.inp

data/grids/tetgrid_24_8_8.inp

data/grids/tetgrid_4_4_4.inp

data/grids/tetgrid_8_8_8.inp

data/grids/tetgrid_LES_10_2_7.inp

data/grids/tetgrid_backward_facing_step_32_8_8.inp

data/grids/tetgrid_bluff_body_32_8_8.inp

data/grids/trigrid-1.inp

data/grids/trigrid-2.inp

data/grids/trigrid-3.inp

data/grids/trigrid-4_4_4.inp

data/grids/trigrid.gid

data/grids/trigrid_1_1_1.inp

doc/policy

doc/releasing

gendoc.sh

include

include/README

include/dolfin.h

install-sh

lib/README

makedist.sh

missing

mkinstalldirs

reconf.sh

src/Makefile

src/Makefile.am

src/Makefile.in

src/common

src/common/Function.hh

src/common/Makefile

src/common/Makefile.am

src/common/Makefile.in

src/common/Parameter.hh

src/common/ParameterList.C

src/common/ParameterList.hh

src/common/Settings.C

src/common/Settings.hh

src/common/SysInfo.C

src/common/SysInfo.hh

src/common/kw_constants.h

src/common/utils.C

src/common/utils.h

src/config

src/config/Makefile.am

src/config/Makefile.in

src/config/config.sh

src/config/dolfin-config.template

src/config/var.tmp

src/fem

src/fem/Discretiser.C

src/fem/Discretiser.hh

src/fem/Equation.C

src/fem/Equation.hh

src/fem/EquationSystem.C

src/fem/EquationSystem.hh

src/fem/FiniteElement.C

src/fem/FiniteElement.hh

src/fem/FunctionSpace.C

src/fem/FunctionSpace.hh

src/fem/GlobalField.C

src/fem/GlobalField.hh

src/fem/LocalField.C

src/fem/LocalField.hh

src/fem/Makefile

src/fem/Makefile.am

src/fem/Makefile.in

src/fem/Problem.C

src/fem/Problem.hh

src/fem/ShapeFunction.C

src/fem/ShapeFunction.hh

src/fem/TetLinFunction.C

src/fem/TetLinFunction.hh

src/fem/TetLinSpace.C

src/fem/TetLinSpace.hh

src/fem/TriLinFunction.C

src/fem/TriLinFunction.hh

src/fem/TriLinSpace.C

src/fem/TriLinSpace.hh

src/grid

src/grid/Cell.C

src/grid/Cell.hh

src/grid/CellType.hh

src/grid/Grid.C

src/grid/Grid.hh

src/grid/Makefile

src/grid/Makefile.am

src/grid/Makefile.in

src/grid/Node.C

src/grid/Node.hh

src/grid/Point.C

src/grid/Point.hh

src/grid/Tetrahedron.C

src/grid/Tetrahedron.hh

src/grid/Triangle.C

src/grid/Triangle.hh

src/io

src/io/Curses.C

src/io/Curses.hh

src/io/DataInfo.C

src/io/DataInfo.hh

src/io/Display.C

src/io/Display.hh

src/io/FileType.C

src/io/FileType.hh

src/io/Input.C

src/io/Input.hh

src/io/Makefile

src/io/Makefile.am

src/io/Makefile.in

src/io/MatlabOld.C

src/io/MatlabOld.hh

src/io/Output.C

src/io/Output.hh

src/io/Terminal.C

src/io/Terminal.hh

src/io/Value.C

src/io/Value.hh

src/io/gid.C

src/io/gid.h

src/io/inp.C

src/io/inp.h

src/io/matlab.C

src/io/matlab.h

src/io/opendx.C

src/io/opendx.h

src/la

src/la/DenseMatrix.C

src/la/DenseMatrix.hh

src/la/DirectSolver.C

src/la/DirectSolver.hh

src/la/KrylovSolver.C

src/la/KrylovSolver.hh

src/la/Makefile

src/la/Makefile.am

src/la/Makefile.in

src/la/SISolver.C

src/la/SISolver.hh

src/la/SparseMatrix.C

src/la/SparseMatrix.hh

src/la/Vector.C

src/la/Vector.hh

src/main

src/main/Makefile

src/main/Makefile.am

src/main/Makefile.in

src/main/dolfin.C

src/main/dolfin.h

src/main/dolfin_modules.h

src/modules

src/modules/Makefile

src/modules/Makefile.am

src/modules/Makefile.in

src/modules/convdiff

src/modules/convdiff/EquationConvDiff2d_cG1cG1.hh

src/modules/convdiff/EquationConvDiff2d_cG1dG0.hh

src/modules/convdiff/EquationConvDiff_cG1cG1.hh

src/modules/convdiff/EquationConvDiff_cG1dG0.hh

src/modules/convdiff/Makefile

src/modules/convdiff/Makefile.am

src/modules/convdiff/Makefile.in

src/modules/convdiff/ProblemConvDiff.C

src/modules/convdiff/ProblemConvDiff.hh

src/modules/convdiff/SettingsConvDiff.hh

src/modules/convdiff/convdiff.module

src/modules/modules.list

src/modules/navierstokes

src/modules/navierstokes/EquationDualPre_cG1cG1.hh

src/modules/navierstokes/EquationDualVel_cG1cG1.hh

src/modules/navierstokes/EquationPre2d_cG1cG1.hh

src/modules/navierstokes/EquationPre_cG1cG1.hh

src/modules/navierstokes/EquationVel2d_cG1cG1.hh

src/modules/navierstokes/EquationVel_cG1cG1.hh

src/modules/navierstokes/Makefile

src/modules/navierstokes/Makefile.am

src/modules/navierstokes/Makefile.in

src/modules/navierstokes/ProblemNS.C

src/modules/navierstokes/ProblemNS.hh

src/modules/navierstokes/SettingsNavierStokes.hh

src/modules/navierstokes/navierstokes.module

src/modules/poisson

src/modules/poisson/EquationPoisson.hh

src/modules/poisson/EquationPoisson2d.hh

src/modules/poisson/Makefile

src/modules/poisson/Makefile.am

src/modules/poisson/Makefile.in

src/modules/poisson/ProblemPoisson.C

src/modules/poisson/ProblemPoisson.hh

src/modules/poisson/SettingsPoisson.hh

src/modules/poisson/poisson.module

src/modules/scanmodules.sh

src/modules/template

src/modules/template/EquationTemplate.hh

src/modules/template/Makefile

src/modules/template/Makefile.am

src/modules/template/Makefile.in

src/modules/template/ProblemTemplate.C

src/modules/template/ProblemTemplate.hh

src/modules/template/README

src/modules/template/SettingsTemplate.hh

src/modules/template/template.module

src/problems

src/problems/Makefile

src/problems/Makefile.am

src/problems/Makefile.in

src/problems/convdiff

src/problems/convdiff/Makefile

src/problems/convdiff/main.C

src/problems/navierstokes

src/problems/navierstokes/Makefile

src/problems/navierstokes/Makefile.am

src/problems/navierstokes/Makefile.in

src/problems/navierstokes/benchmark

src/problems/navierstokes/benchmark/Makefile

src/problems/navierstokes/benchmark/README

src/problems/navierstokes/benchmark/main.C

src/problems/poisson

src/problems/poisson/Makefile

src/problems/poisson/main.C

src/problems/schroedinger

src/problems/schroedinger/Makefile

src/problems/schroedinger/main.C

src/problems/template

src/problems/template/Makefile

src/problems/template/main.C

src/test

src/test/Makefile

src/test/main.C

src/utils

src/utils/Makefile

src/utils/Makefile.am

src/utils/Makefile.in

src/utils/inp2dx

src/utils/inp2dx/Makefile

src/utils/inp2dx/Makefile.am

src/utils/inp2dx/Makefile.in

src/utils/inp2dx/inp2dx.c

src/utils/inp2dx/inp2dxanim.c

src/utils/matlab

src/utils/matlab/writegid.m

src/utils/matlab/writeinp.m

src/utils/opendx

src/utils/opendx/dolfin.cfg

src/utils/opendx/dolfin.net

todo

Show diffs side-by-side

added added

removed removed

src/modules/navierstokes/ProblemNS.C

// Licensed under the GNU GPL Version 2.

#include "ProblemNS.hh"

#include "EquationVel_cG1cG1.hh"

#include "EquationPre_cG1cG1.hh"

#include "EquationVel2d_cG1cG1.hh"

#include "EquationPre2d_cG1cG1.hh"

#include <Grid.hh>

#include <Vector.hh>

#include <SparseMatrix.hh>

#include <MatlabOld.hh>

#include <unistd.h>

//----------------------------------------------------------------------------------------------------

ProblemNS::ProblemNS(Grid *grid) : Problem(grid)

{

// Get space dimension

settings->Get("space dimension",&nsd);

if ( (nsd!=2) && (nsd!=3) )

display->InternalError("ProblemNS::ProblemNS()","Only implemented for 2d and 3d");

// Get inner variables

settings->Get("write residuals", &write_residuals);

settings->Get("write data", &write_data);

settings->Get("compute projections",&compute_projections);

settings->Get("write couette pertubations",&pert_couette);

settings->Get("write poiseuille pertubations",&pert_poiseuille);

pertubation = false;

if ( pert_poiseuille || pert_couette ) pertubation = true;

// Initialize Opendx objects

char output_file[DOLFIN_LINELENGTH];

settings->Get("output file primal",output_file);

//opendx_residuals = new OpenDX(output_file,4,nsd,1,1,1);

//opendx_residuals->SetLabel(0,"R1");

//opendx_residuals->SetLabel(1,"R2");

//opendx_residuals->SetLabel(1,"R3");

//opendx_residuals->SetLabel(1,"R4");

// Initialise Output for 2d problems

if ( nsd == 2 ){

output = new Output("navier_stokes_2d.m",2,nsd,1);

}

// Initialize scalar output data

if ( nsd!=3 ){

write_data = false;

write_residuals = false;

pertubation = false;

write_reynolds_stresses = false;

}

no_data_entries = 0;

no_residual_entries = 0;

no_pertubation_entries = 0;

no_re_stress_entries = 0;

if ( write_data ) no_data_entries = 39;

if ( write_residuals) no_residual_entries = 12;

if ( pertubation ) no_pertubation_entries = 24;

if ( write_reynolds_stresses ) no_re_stress_entries = 27;

matlab_data = NULL;

if ( write_data || write_residuals || pertubation || write_reynolds_stresses ){

matlab_data = new MatlabOld(no_data_entries + no_residual_entries + no_pertubation_entries + no_re_stress_entries);

}

if ( write_data ){

matlab_data->SetLabel(0,"|| u ||_1");

matlab_data->SetLabel(1,"|| u ||_2");

matlab_data->SetLabel(2,"|| u ||_{max}");

matlab_data->SetLabel(3,"| u |_{1,1}");

matlab_data->SetLabel(4,"| u |_{1,2}");

matlab_data->SetLabel(5,"| u |_{1,max}");

matlab_data->SetLabel(6,"|| u ||_{1,1}");

matlab_data->SetLabel(7,"|| u ||_{1,2}");

matlab_data->SetLabel(8,"|| u ||_{1,max}");

matlab_data->SetLabel(9,"| u |_{2,1}");

matlab_data->SetLabel(10,"| u |_{2,2}");

matlab_data->SetLabel(11,"| u |_{2,max}");

matlab_data->SetLabel(12,"|| u ||_{2,1}");

matlab_data->SetLabel(13,"|| u ||_{2,2}");

matlab_data->SetLabel(14,"|| u ||_{2,max}");

matlab_data->SetLabel(15,"|| p ||_1");

matlab_data->SetLabel(16,"|| p ||_2");

matlab_data->SetLabel(17,"|| p ||_{max}");

matlab_data->SetLabel(18,"| p |_{1,1}");

matlab_data->SetLabel(19,"| p |_{1,2}");

matlab_data->SetLabel(20,"| p |_{1,max}");

matlab_data->SetLabel(21,"|| p ||_{1,1}");

matlab_data->SetLabel(22,"|| p ||_{1,2}");

matlab_data->SetLabel(23,"|| p ||_{1,max}");

matlab_data->SetLabel(24,"| p |_{2,1}");

matlab_data->SetLabel(25,"| p |_{2,2}");

matlab_data->SetLabel(26,"| p |_{2,max}");

matlab_data->SetLabel(27,"|| p ||_{2,1}");

matlab_data->SetLabel(28,"|| p ||_{2,2}");

matlab_data->SetLabel(29,"|| p ||_{2,max}");

matlab_data->SetLabel(30,"|| div u ||_1");

100

matlab_data->SetLabel(31,"|| div u ||_2");

101

matlab_data->SetLabel(32,"|| div u ||_{max}");

102

matlab_data->SetLabel(33,"|| du/dt ||_1");

103

matlab_data->SetLabel(34,"|| du/dt ||_2");

104

matlab_data->SetLabel(35,"|| du/dt ||_{max}");

105

matlab_data->SetLabel(36,"|| dp/dt ||_1");

106

matlab_data->SetLabel(37,"|| dp/dt ||_2");

107

matlab_data->SetLabel(38,"|| dp/dt ||_{max}");

108

}

109

110

if ( write_residuals ){

111

matlab_data->SetLabel(no_data_entries+0,"|| R1 ||_1");

112

matlab_data->SetLabel(no_data_entries+1,"|| R1 ||_2");

113

matlab_data->SetLabel(no_data_entries+2,"|| R1 ||_{max}");

114

matlab_data->SetLabel(no_data_entries+3,"|| R2 ||_1");

115

matlab_data->SetLabel(no_data_entries+4,"|| R2 ||_2");

116

matlab_data->SetLabel(no_data_entries+5,"|| R2 ||_{max}");

117

matlab_data->SetLabel(no_data_entries+6,"|| R3 ||_1");

118

matlab_data->SetLabel(no_data_entries+7,"|| R3 ||_2");

119

matlab_data->SetLabel(no_data_entries+8,"|| R3 ||_{max}");

120

matlab_data->SetLabel(no_data_entries+9,"|| R4 ||_1");

121

matlab_data->SetLabel(no_data_entries+10,"|| R4 ||_2");

122

matlab_data->SetLabel(no_data_entries+11,"|| R4 ||_{max}");

123

}

124

125

if ( pertubation ){

126

matlab_data->SetLabel(no_data_entries+no_residual_entries+0,"|| u1_{pert} ||_2");

127

matlab_data->SetLabel(no_data_entries+no_residual_entries+1,"|| u1_{pert} ||_{max}");

128

matlab_data->SetLabel(no_data_entries+no_residual_entries+2,"|| u2_{pert} ||_2");

129

matlab_data->SetLabel(no_data_entries+no_residual_entries+3,"|| u2_{pert} ||_{max}");

130

matlab_data->SetLabel(no_data_entries+no_residual_entries+4,"|| u3_{pert} ||_2");

131

matlab_data->SetLabel(no_data_entries+no_residual_entries+5,"|| u3_{pert} ||_{max}");

132

matlab_data->SetLabel(no_data_entries+no_residual_entries+6,"|| du_1/dx_1 ||_2");

133

matlab_data->SetLabel(no_data_entries+no_residual_entries+7,"|| du_1/dx_1 ||_{max}");

134

matlab_data->SetLabel(no_data_entries+no_residual_entries+8,"|| du_1/dx_2 ||_2");

135

matlab_data->SetLabel(no_data_entries+no_residual_entries+9,"|| du_1/dx_2 ||_{max}");

136

matlab_data->SetLabel(no_data_entries+no_residual_entries+10,"|| du_1/dx_3 ||_2");

137

matlab_data->SetLabel(no_data_entries+no_residual_entries+11,"|| du_1/dx_3 ||_{max}");

138

matlab_data->SetLabel(no_data_entries+no_residual_entries+12,"|| du_2/dx_1 ||_2");

139

matlab_data->SetLabel(no_data_entries+no_residual_entries+13,"|| du_2/dx_1 ||_{max}");

140

matlab_data->SetLabel(no_data_entries+no_residual_entries+14,"|| du_2/dx_2 ||_2");

141

matlab_data->SetLabel(no_data_entries+no_residual_entries+15,"|| du_2/dx_2 ||_{max}");

142

matlab_data->SetLabel(no_data_entries+no_residual_entries+16,"|| du_2/dx_3 ||_2");

143

matlab_data->SetLabel(no_data_entries+no_residual_entries+17,"|| du_2/dx_3 ||_{max}");

144

matlab_data->SetLabel(no_data_entries+no_residual_entries+18,"|| du_3/dx_1 ||_2");

145

matlab_data->SetLabel(no_data_entries+no_residual_entries+19,"|| du_3/dx_1 ||_{max}");

146

matlab_data->SetLabel(no_data_entries+no_residual_entries+20,"|| du_3/dx_2 ||_2");

147

matlab_data->SetLabel(no_data_entries+no_residual_entries+21,"|| du_3/dx_2 ||_{max}");

148

matlab_data->SetLabel(no_data_entries+no_residual_entries+22,"|| du_3/dx_3 ||_2");

149

matlab_data->SetLabel(no_data_entries+no_residual_entries+23,"|| du_3/dx_3 ||_{max}");

150

}

151

152

if ( write_reynolds_stresses ){

153

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+0,"|| tau_11 ||_1");

154

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+1,"|| tau_12 ||_1");

155

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+2,"|| tau_13 ||_1");

156

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+3,"|| tau_22 ||_1");

157

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+4,"|| tau_23 ||_1");

158

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+5,"|| tau_33 ||_1");

159

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+6,"|| tau_11 ||_2");

160

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+7,"|| tau_12 ||_2");

161

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+8,"|| tau_13 ||_2");

162

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+9,"|| tau_22 ||_2");

163

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+10,"|| tau_23 ||_2");

164

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+11,"|| tau_33 ||_2");

165

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+12,"|| tau_11 ||_{max}");

166

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+13,"|| tau_12 ||_{max}");

167

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+14,"|| tau_13 ||_{max}");

168

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+15,"|| tau_22 ||_{max}");

169

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+16,"|| tau_23 ||_{max}");

170

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+17,"|| tau_33 ||_{max}");

171

172

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+18,"|| div tau_1 ||_1");

173

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+19,"|| div tau_2 ||_1");

174

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+20,"|| div tau_3 ||_1");

175

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+21,"|| div tau_1 ||_2");

176

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+22,"|| div tau_2 ||_2");

177

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+23,"|| div tau_3 ||_2");

178

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+24,"|| div tau_1 ||_{max}");

179

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+25,"|| div tau_2 ||_{max}");

180

matlab_data->SetLabel(no_data_entries+no_residual_entries+no_pertubation_entries+26,"|| div tau_3 ||_{max}");

181

}

182

}

183

//----------------------------------------------------------------------------------------------------

184

ProblemNS::~ProblemNS()

185

{

186

}

187

//-----------------------------------------------------------------------------

188

const char *ProblemNS::Description()

189

{

190

return "Navier-Stokes equations";

191

}

192

//----------------------------------------------------------------------------------------------------

193

void ProblemNS::Solve()

194

{

195

display->Status(0,"Solving Navier-Stokes primal problem");

196

197

// Set Reynolds number

198

settings->Get("reynolds number",&Re);

199

200

// Reset variables

201

Reset();

202

203

// Set FEM discretisation

204

Equation *eq_vel;

205

Equation *eq_pre;

206

207

if ( nsd == 3 ){

208

eq_vel = new EquationVel_cG1cG1;

209

eq_pre = new EquationPre_cG1cG1;

210

}

211

else{

212

eq_vel = new EquationVel2d_cG1cG1;

213

eq_pre = new EquationPre2d_cG1cG1;

214

}

215

216

Discretiser discrete_vel(grid,eq_vel);

217

Discretiser discrete_pre(grid,eq_pre);

218

219

// Set number of equations

220

noeq_vel = eq_vel->GetNoEq();

221

noeq_pre = eq_pre->GetNoEq();

222

noeq = noeq_vel + noeq_pre;

223

224

// Initialize solution variables

225

InitSolutionVariables();

226

227

GlobalField upTS_field(grid,upTS,noeq_vel);

228

GlobalField ppTS_field(grid,ppTS,noeq_pre);

229

GlobalField upNL_field(grid,upNL,noeq_vel);

230

GlobalField ppNL_field(grid,ppNL,noeq_pre);

231

GlobalField u_field(grid,u,noeq_vel);

232

GlobalField p_field(grid,p,noeq_pre);

233

234

GlobalField uc_field(grid,u_coarse,noeq_vel);

235

GlobalField pc_field(grid,u_coarse,noeq_pre);

236

GlobalField uf_field(grid,u_fine,noeq_vel);

237

GlobalField pf_field(grid,u_fine,noeq_pre);

238

239

GlobalField *fx;

240

GlobalField *fy;

241

GlobalField *fz;

242

243

// A global field for saving the solution

244

up_field = new GlobalField(&u_field,&p_field);

245

if ( nsd == 3 )

246

up_field->SetSize(2,3,1);

247

else

248

up_field->SetSize(2,2,1);

249

up_field->SetLabel("u","Velocity",0);

250

up_field->SetLabel("p","Pressure",1);

251

252

// A global field for the right-hand side

253

fx = new GlobalField(grid,"fx");

254

fy = new GlobalField(grid,"fy");

255

if (nsd==3) fz = new GlobalField(grid,"fz");

256

257

GlobalField tau_field(grid,tau,6);

258

259

int field_no = 0;

260

for (int i=0;i<noeq_vel;i++) eq_vel->AttachField(field_no++,&upTS_field,i);

261

eq_vel->AttachField(field_no++,&ppTS_field);

262

for (int i=0;i<noeq_vel;i++) eq_vel->AttachField(field_no++,&upNL_field,i);

263

eq_vel->AttachField(field_no++,&ppNL_field);

264

for (int i=0;i<noeq_vel;i++) eq_vel->AttachField(field_no++,&u_field,i);

265

eq_vel->AttachField(field_no++,&p_field);

266

for (int i=0;i<noeq_vel;i++) eq_vel->AttachField(field_no++,&uf_field,i);

267

eq_vel->AttachField(field_no++,&pf_field);

268

for (int i=0;i<noeq_vel;i++) eq_vel->AttachField(field_no++,&uc_field,i);

269

eq_vel->AttachField(field_no++,&pc_field);

270

eq_vel->AttachField(field_no++,fx);

271

eq_vel->AttachField(field_no++,fy);

272

if (nsd==3) eq_vel->AttachField(field_no++,fz);

273

274

field_no = 0;

275

for (int i=0;i<noeq_vel;i++) eq_pre->AttachField(field_no++,&upTS_field,i);

276

eq_pre->AttachField(field_no++,&ppTS_field);

277

for (int i=0;i<noeq_vel;i++) eq_pre->AttachField(field_no++,&upNL_field,i);

278

eq_pre->AttachField(field_no++,&ppNL_field);

279

for (int i=0;i<noeq_vel;i++) eq_pre->AttachField(field_no++,&u_field,i);

280

eq_pre->AttachField(field_no++,&p_field);

281

for (int i=0;i<noeq_vel;i++) eq_pre->AttachField(field_no++,&uf_field,i);

282

eq_pre->AttachField(field_no++,&pf_field);

283

for (int i=0;i<noeq_vel;i++) eq_pre->AttachField(field_no++,&uc_field,i);

284

eq_pre->AttachField(field_no++,&pc_field);

285

eq_pre->AttachField(field_no++,fx);

286

eq_pre->AttachField(field_no++,fy);

287

if (nsd==3) eq_pre->AttachField(field_no++,fz);

288

289

// Reassemble the momentum matrix or not

290

ReAssembleVelMatrix(true);

291

292

// Set initial data

293

//initial_data = zero_initial_data;

294

initial_data = couette_flow;

295

SetInitialData();

296

297

// Set time t to starting time T0

298

t = T0;

299

300

// Write data to file

301

WriteDataToFile();

302

303

// Start time stepping

304

display->Status(0,"Start computation at t = %.1f",t);

305

while ( t <= T ){

306

307

// Set time step

308

GetMeshSize();

309

dt = GetTimeStep();

310

t += dt;

311

312

// Set time and time step for equations

313

eq_vel->SetTime(t);

314

eq_pre->SetTime(t);

315

eq_vel->SetTimeStep(dt);

316

eq_pre->SetTimeStep(dt);

317

318

// Save solution if we should

319

if ( SaveSample() ) WriteSolutionToFile();

320

321

// Set turbulent inflow condition

322

SetTurbulentInflow(0.04);

323

324

// Write some output for gdisp

325

progress = (t-T0) / (T-T0);

326

display->Progress(0,progress,"Primal solution: %.1f %%. Currently at t = %.4f, dt = %.4f",progress*100.0,t,dt);

327

328

// Set solution from previous time step

329

upTS->CopyFrom(u);

330

// Set solution from previous outer non linear iteration

331

upNL->CopyFrom(u);

332

333

// Solve non linear problem

334

SolveNonLinearProblemUsingFixPoint(&discrete_vel,&discrete_pre,eq_vel,eq_pre);

335

//SolveNonLinearProblemUsingGMRES();

336

337

// Write data to file

338

WriteDataToFile();

339

340

// Check if the solution is stationary

341

if ( CheckIfStationary() ) break;

342

}

343

344

// Solution computed to end time

345

display->Status(0,"Solution has been computed to endtime T = %f",t);

346

347

// Delete solution vectors

348

DeleteSolutionVariables();

349

350

delete eq_vel;

351

delete eq_pre;

352

353

delete fx;

354

delete fy;

355

delete fz;

356

}

357

//----------------------------------------------------------------------------------------------------

358

void ProblemNS::SolveNonLinearProblemUsingFixPoint(Discretiser *discrete_vel, Discretiser *discrete_pre,

359

Equation *eq_vel, Equation *eq_pre)

360

{

361

real pre_res_norm_2, pre_res_norm_inf, vel_res_norm_2, vel_res_norm_inf;

362

363

// Start non linear iteration

364

for (NL_it=0; NL_it < max_no_NL_iter; NL_it++){

365

366

// Display information

367

display->Status(0,"Starting non-linear iteration %d.",NL_it+1);

368

369

// Compute Reynolds stresses - for subgrid model

370

if ( compute_reynolds_stresses ){

371

display->Status(5,"Estimate Reynolds stresses (Update subgrid model)");

372

//scale_expol.EstimateReynoldsStresses( u, tau, &(*GRIDREL)(0), &(*GRIDREL)(1),

373

// &(*GRIDREL)(2), &(*GRIDREL)(3), &(*GRIDREL)(4) );

374

}

375

376

// Compute projections - for subgrid model

377

if ( compute_projections ){

378

display->Status(5,"Extract coarse and fine scale field");

379

//(*GRIDREL)(0).projNodalData(u,u_coarse);

380

//(*GRIDREL)(0).projNodalData(u_coarse);

381

//for (int i=0; i < no_nodes*noeq_vel; i++ ) u_fine(i) = u(i) - u_coarse(i);

382

}

383

384

// Start inner iteration

385

for (UP_it=0; UP_it < max_no_UP_iter; UP_it++){

386

387

// Display information

388

display->Status(1,"Starting u-p iteration %d.",UP_it+1);

389

390

// Solve continuity equation

391

display->Status(2,"Solving for the pressure.");

392

TimestepPressure(discrete_pre,eq_pre);

393

394

// Solve momentum equation

395

display->Status(2,"Solving for the velocity.");

396

TimestepVelocity(discrete_vel,eq_vel);

397

398

// Re assemble the momentum matrix or not

399

ReAssembleVelMatrix(false);

400

401

// Compute residuals for inner iteration

402

ComputeDiscreteResidual(discrete_pre,eq_pre,A_pre,p,pre_res_norm_2, pre_res_norm_inf);

403

ComputeDiscreteResidual(discrete_vel,eq_vel,A_vel,u,vel_res_norm_2, vel_res_norm_inf);

404

405

// Check if inner loop has converged

406

if ( CheckUPConvergence(pre_res_norm_2, pre_res_norm_inf, vel_res_norm_2, vel_res_norm_inf) ) break;

407

}

408

409

// Re assemble the momentum matrix or not

410

ReAssembleVelMatrix(true);

411

412

// Compute residuals for non linear outer iteration

413

upNL->CopyFrom(u);

414

ComputeDiscreteResidual(discrete_pre,eq_pre,A_pre,p,pre_res_norm_2, pre_res_norm_inf);

415

ComputeDiscreteResidual(discrete_vel,eq_vel,A_vel,u,vel_res_norm_2,vel_res_norm_inf);

416

417

// Re assemble the momentum matrix or not

418

ReAssembleVelMatrix(false);

419

420

// Check if outer loop has converged

421

if ( CheckNonLinearConvergence(pre_res_norm_2, pre_res_norm_inf, vel_res_norm_2, vel_res_norm_inf) ) break;

422

}

423

}

424

//----------------------------------------------------------------------------------------------------

425

void SolveNonLinearProblemUsingGMRES(Discretiser *discrete_vel, Discretiser *discrete_pre,

426

Equation *eq_vel, Equation *eq_pre)

427

{

428

display->InternalError("ProblemNS::SolveNonLinearProblemUsingGMRES()","Not implemented yet");

429

}

430

//----------------------------------------------------------------------------------------------------

431

void ProblemNS::TimestepVelocity(Discretiser *discrete, Equation *equation)

432

{

433

Vector b(no_nodes*equation->GetNoEq());

434

435

// assemble matrix if not assembled, or reassemble if we should

436

if ( !A_vel_assembled ){

437

display->Status(3,"Assemble matrix A");

438

discrete->AssembleLHS(A_vel);

439

A_vel_assembled = true;

440

} else if ( A_vel_assembled && re_assemble_matrix ){

441

display->Status(3,"Re-assembling matrix A");

442

discrete->AssembleLHS(A_vel);

443

}

444

445

// assemble vector

446

display->Status(3,"Assemble vector b");

447

discrete->AssembleRHS(&b);

448

449

// set bc

450

display->Message(2,"Load norm: %f",b.Norm(2));

451

discrete->SetBoundaryConditions(A_vel,&b);

452

display->Message(2,"Load norm: %f (with b.c.)",b.Norm(2));

453

454

// solve discrete system

455

display->Status(3,"Solve algebraic system of equations");

456

457

KrylovSolver solver;

458

solver.SetMethod(gmres);

459

solver.Solve(A_vel,u,&b);

460

461

// compute norm of solution

462

display->Message(4,"Norm of u = %f",u->Norm(2));

463

display->Message(4,"Norm of A = %f",A_vel->Norm());

464

}

465

//----------------------------------------------------------------------------------------------------

466

void ProblemNS::TimestepPressure(Discretiser *discrete, Equation *equation)

467

{

468

Vector b(no_nodes*equation->GetNoEq());

469

470

if ( !A_pre_assembled ){

471

display->Status(3,"Assemble matrix A and vector b.");

472

discrete->AssembleLHS(A_pre);

473

A_pre_assembled = true;

474

}

475

476

display->Status(3,"Assemble vector b");

477

discrete->AssembleRHS(&b);

478

479

// set dirichlet bc

480

display->Message(2,"Load norm: %f",b.Norm(2));

481

discrete->SetBoundaryConditions(A_pre,&b);

482

display->Message(2,"Load norm: %f (with b.c.)",b.Norm(2));

483

484

display->Status(3,"Solve algebraic system of equations");

485

486

KrylovSolver solver;

487

solver.SetMethod(gmres);

488

solver.Solve(A_pre,p,&b);

489

490

// compute norm of solution

491

display->Message(4,"Norm of p = %f",p->Norm(2));

492

display->Message(4,"Norm of A = %f",A_pre->Norm());

493

}

494

//----------------------------------------------------------------------------------------------------

495

void ProblemNS::ComputeDiscreteResidual(Discretiser *discrete, Equation *equation,

496

SparseMatrix *A, Vector *U, real& norm_2, real& norm_inf)

497

{

498

Vector b(no_nodes*equation->GetNoEq());

499

500

// Only for momentum equation

501

if ( re_assemble_matrix && (equation->GetNoEq() == noeq_vel) ){

502

display->Status(3,"Assemble matrix A");

503

discrete->AssembleLHS(A);

504

}

505

506

display->Status(3,"Assemble vector b");

507

discrete->AssembleRHS(&b);

508

509

// set dirichlet bc

510

display->Message(2,"Load norm: %f",b.Norm(2));

511

discrete->SetBoundaryConditions(A,&b);

512

display->Message(2,"Load norm: %f (with b.c.)",b.Norm(2));

513

514

display->Message(4,"Norm of solution = %f",U->Norm(2));

515

display->Message(4,"Norm of matrix = %f",A->Norm());

516

517

display->Message(3,"Compute discrete residual");

518

Vector AU(no_nodes*equation->GetNoEq());

519

A->Mult(U,&AU);

520

521

Vector res(no_nodes*equation->GetNoEq());

522

for (int i=0; i < no_nodes*equation->GetNoEq(); i++) res.Set(i,b(i) - AU(i));

523

524

norm_2 = res.Norm(2);

525

norm_inf = res.Norm(0);

526

}

527

//----------------------------------------------------------------------------------------------------

528

void ProblemNS::WriteDataToFile()

529

{

530

display->Status(0,"Write data to file");

531

if ( write_data ) WriteData();

532

if ( write_residuals ) CompContResidual();

533

if ( pertubation ) WritePertubationData();

534

if ( write_reynolds_stresses ) WriteReynoldsStresses();

535

}

536

//----------------------------------------------------------------------------------------------------

537

void ProblemNS::WriteSolutionToFile()

538

{

539

display->Status(0,"Write solution to file");

540

541

up_field->Save(t);

542

}

543

//----------------------------------------------------------------------------------------------------

544

void ProblemNS::ReAssembleVelMatrix(bool re_assemble_matrix)

545

{

546

this->re_assemble_matrix = re_assemble_matrix;

547

}

548

//-----------------------------------------------------------------------------

549

bool ProblemNS::SaveSample()

550

{

551

// Check if we should save a sample

552

bool save_sample = false;

553

554

int ops;

555

settings->Get("output samples", &ops);

556

real M = real(ops);

557

558

real n = current_frame;

559

real sample_time;

560

561

if (dt < 0) {

562

sample_time = T - (T - T0) / M * n;

563

if (t <= (sample_time + DOLFIN_EPS))

564

save_sample = true;

565

} else {

566

sample_time = T0 + (T - T0) / M * n;

567

if (t >= (sample_time - DOLFIN_EPS))

568

save_sample = true;

569

}

570

571

if (save_sample)

572

current_frame += 1;

573

574

return (save_sample);

575

}

576

//----------------------------------------------------------------------------------------------------

577

void ProblemNS::InitSolutionVariables()

578

{

579

// Initialize matrices

580

A_pre = new SparseMatrix();

581

A_vel = new SparseMatrix();

582

583

// Initialize solution vectors

584

u = new Vector(no_nodes*noeq_vel);

585

p = new Vector(no_nodes*noeq_pre);

586

587

upTS = new Vector(no_nodes*noeq_vel);

588

ppTS = new Vector(no_nodes*noeq_pre);

589

upNL = new Vector(no_nodes*noeq_vel);

590

ppNL = new Vector(no_nodes*noeq_pre);

591

592

upTS->SetToConstant(0.0);

593

ppTS->SetToConstant(0.0);

594

upNL->SetToConstant(0.0);

595

upNL->SetToConstant(0.0);

596

597

print_vector = new Vector(print_no_nodes*noeq);

598

599

// Initialize projections (for subgrid modeling)

600

u_fine = new Vector(1);

601

p_fine = new Vector(1);

602

u_coarse = new Vector(1);

603

p_coarse = new Vector(1);

604

605

// FIXME: Should only be initialized if necessary

606

u_fine->Resize(no_nodes*noeq_vel);

607

p_fine->Resize(no_nodes*noeq_pre);

608

u_coarse->Resize(no_nodes*noeq_vel);

609

p_coarse->Resize(no_nodes*noeq_pre);

610

611

u_fine->SetToConstant(0.0);

612

p_fine->SetToConstant(0.0);

613

u_coarse->SetToConstant(0.0);

614

p_coarse->SetToConstant(0.0);

615

616

tau = new Vector;

617

if ( compute_reynolds_stresses ){

618

tau->Resize(no_nodes*6);

619

}

620

621

tau->SetToConstant(0.0);

622

}

623

//----------------------------------------------------------------------------------------------------

624

void ProblemNS::DeleteSolutionVariables()

625

{

626

// Delete matrices

627

delete A_pre;

628

delete A_vel;

629

630

// Delete solution vectors

631

delete u;

632

delete p;

633

delete upTS;

634

delete ppTS;

635

delete upNL;

636

delete ppNL;

637

delete print_vector;

638

639

if ( compute_projections ){

640

delete u_fine;

641

delete p_fine;

642

delete u_coarse;

643

delete p_coarse;

644

}

645

646

delete tau;

647

}

648

//-----------------------------------------------------------------------------

649

void ProblemNS::Reset()

650

{

651

// Reset gdisp

652

progress = 0.0;

653

654

// Reset OpenDX objects

655

//opendx_residuals->Reset();

656

657

// Set output frame to zero

658

current_frame = 0;

659

660

// Reset Matlab file

661

if ( matlab_data ) matlab_data->Reset();

662

663

// Set tolerances for solution method

664

SetTolerances();

665

666

// Initialize subgrid modeling utilities

667

SetSubgridModelingUtilities();

668

669

// Update output settings

670

settings->Get("write residuals", &write_residuals);

671

settings->Get("write data", &write_data);

672

673

settings->Get("start time", &T0);

674

settings->Get("final time", &T);

675

676

// Update internal variables

677

no_nodes = grid->GetNoNodes();

678

no_cells = grid->GetNoCells();

679

680

print_no_nodes = grid->GetNoNodes();

681

print_no_cells = grid->GetNoCells();

682

683

A_vel_assembled = false;

684

A_pre_assembled = false;

685

}

686

//-----------------------------------------------------------------------------

687

bool ProblemNS::CheckIfStationary()

688

{

689

// Check if the solution is stationary

690

bool stationary_solution = false;

691

692

real maxnorm_du = 0.0;

693

for (int i = 0; i < no_nodes; i++) {

694

if (fabs(u->Get(i)-upTS->Get(i)) > maxnorm_du)

695

maxnorm_du = fabs(u->Get(i) - upTS->Get(i));

696

}

697

if (maxnorm_du / dt < stat_tol) {

698

display->Message(0,"Stationary limit reached, at time t = %.4f, max du/dt = %1.4e",t,maxnorm_du/dt);

699

stationary_solution = true;

700

}

701

702

return stationary_solution;

703

}

704

//----------------------------------------------------------------------------------------------------

705

bool ProblemNS::CheckNonLinearConvergence(real pre_res_norm_2, real pre_res_norm_inf,

706

real vel_res_norm_2, real vel_res_norm_inf )

707

{

708

bool non_linear_convergence = false;

709

710

NL_res_norm_2 = sqrt( sqr(pre_res_norm_2) + sqr(vel_res_norm_2) );

711

if ( pre_res_norm_inf > vel_res_norm_inf ) NL_res_norm_inf = pre_res_norm_inf;

712

else NL_res_norm_inf = vel_res_norm_inf;

713

714

display->Message(4,"l2-norm of NL-residuals at NL iteration %i: ",NL_it+1);

715

display->Message(4,"%1.4e (momentum equation)",vel_res_norm_2);

716

display->Message(4,"%1.4e (continuity equation)",pre_res_norm_2);

717

display->Message(4,"%1.4e (total l2-residual)",NL_res_norm_2);

718

display->Message(4,"%1.4e (total max-residual)",NL_res_norm_inf);

719

720

//display->Value("Non linear residual",type_real,NL_res_norm_2);

721

display->Status(2,"Non linear residual: %1.4e",NL_res_norm_2);

722

723

if ( NL_res_norm_2 < NL_tol ) non_linear_convergence = true;

724

if ( (NL_res_norm_2 > NL_max_tol) && (NL_it == (max_no_NL_iter-1)) ){

725

display->Error("NL residual too large");

726

}

727

728

return non_linear_convergence;

729

}

730

//----------------------------------------------------------------------------------------------------

731

bool ProblemNS::CheckUPConvergence(real pre_res_norm_2, real pre_res_norm_inf,

732

real vel_res_norm_2, real vel_res_norm_inf )

733

{

734

bool UP_convergence = false;

735

736

UP_res_norm_2 = sqrt( sqr(pre_res_norm_2) + sqr(vel_res_norm_2) );

737

if ( pre_res_norm_inf > vel_res_norm_inf ) UP_res_norm_inf = pre_res_norm_inf;

738

else UP_res_norm_inf = vel_res_norm_inf;

739

740

display->Message(4,"l2-norm of UP-residuals at UP iteration %i: ",UP_it+1);

741

display->Message(4,"%1.4e (momentum equation)",vel_res_norm_2);

742

display->Message(4,"%1.4e (continuity equation)",pre_res_norm_2);

743

display->Message(4,"%1.4e (total l2-residual)",UP_res_norm_2);

744

display->Message(4,"%1.4e (total max-residual)",UP_res_norm_inf);

745

746

//display->Value("UP residual",type_real,NL_res_norm_2);

747

display->Status(2,"UP residual: %1.4e",UP_res_norm_2);

748

749

if ( UP_res_norm_2 < UP_tol ) UP_convergence = true;

750

if ( (UP_res_norm_2 > UP_max_tol) && (UP_it == (max_no_UP_iter-1)) ){

751

display->Error("UP residual too large");

752

}

753

754

return UP_convergence;

755

}

756

//----------------------------------------------------------------------------------------------------

757

void ProblemNS::SetSubgridModelingUtilities()

758

{

759

settings->Get("compute reynolds stresses",&compute_reynolds_stresses);

760

settings->Get("write reynolds stresses",&write_reynolds_stresses);

761

762

subgrid_model_on_divergence_form = true;

763

764

// N = (sqr(nsd)+nsd)/2;

765

if ( subgrid_model_on_divergence_form ) N = 6;

766

// N = sqr(nsd);

767

else N = 9;

768

// 3 levels

769

// N *= 3;

770

}

771

//----------------------------------------------------------------------------------------------------

772

void ProblemNS::SetInitialData()

773

{

774

display->Status(1,"Set initial conditions");

775

776

initial_data = zero_initial_data;

777

778

Point pnt;

779

u->SetToConstant(0.0);

780

p->SetToConstant(0.0);

781

782

char tmp[256];

783

real x,y,z,pert,noRolls;

784

real pi = DOLFIN_PI;

785

786

switch ( initial_data ){

787

case zero_initial_data:

788

break;

789

case white_noise:

790

for (int i=0; i < no_nodes*noeq_vel; i++ ){

791

u->Set(i,0.1 * 2.0*(drand48() - 0.5));

792

}

793

break;

794

case poiseuille_flow:

795

if (nsd != 3) display->InternalError("ProlemNS::SetInitialData()","Only implemented for 3d");

796

pert = 0.0;

797

noRolls = 1.0;

798

for (int i=0; i < no_nodes; i++ ){

799

pnt = grid->GetNode(i)->GetCoord();

800

x = real(pnt.x);

801

y = real(pnt.y);

802

z = real(pnt.z);

803

u->Set(i*noeq_vel+0,16.0*y*(1.0-y)*z*(1-z));

804

u->Set(i*noeq_vel+1,pert * ( sin(noRolls*2.0*pi*y) * cos(noRolls*pi*z) ));

805

u->Set(i*noeq_vel+2,pert * ( - cos(noRolls*2.0*pi*y) * sin(noRolls*pi*z) ));

806

}

807

break;

808

case couette_flow:

809

if (nsd != 3) display->InternalError("ProlemNS::SetInitialData()","Only implemented for 3d");

810

pert = 0.0;

811

for (int i=0; i < no_nodes; i++ ){

812

pnt = grid->GetNode(i)->GetCoord();

813

x = real(pnt.x);

814

y = real(pnt.y);

815

z = real(pnt.z);

816

u->Set(i*noeq_vel+0,2.0*y-1.0);

817

u->Set(i*noeq_vel+1,pert * ( sin(2.0*pi*y) * cos(pi*z) ));

818

u->Set(i*noeq_vel+2,pert * ( - cos(2.0*pi*y) * sin(pi*z) ));

819

}

820

break;

821

default:

822

display->Error("The initial condition is not implemented.");

823

}

824

}

825

//----------------------------------------------------------------------------------------------------

826

void ProblemNS::WriteData()

827

{

828

829

int el,nod,i;

830

831

// Initializing Fields

832

Field Ufld_1(grid);

833

Field Ufld_2(grid);

834

Field Ufld_3(grid);

835

Field Pfld(grid);

836

837

Ufld_1.InitByNodes(U,0);

838

Ufld_2.InitByNodes(U,1);

839

Ufld_3.InitByNodes(U,2);

840

Pfld.InitByNodes(U,3);

841

842

// Compute norms of solution

843

real L1_norm_U = Ufld_1.Lp_norm(1) + Ufld_2.Lp_norm(1) + Ufld_3.Lp_norm(1);

844

real L2_norm_U = sqrt( sqr(Ufld_1.Lp_norm(2)) + sqr(Ufld_2.Lp_norm(2)) + sqr(Ufld_3.Lp_norm(2)) );

845

real Lmax_norm_U = max( max(Ufld_1.Lmax_norm(),Ufld_2.Lmax_norm()) , Ufld_3.Lmax_norm() );

846

847

real W11_seminorm_U = Ufld_1.Wkp_seminorm(1,1) + Ufld_2.Wkp_seminorm(1,1) + Ufld_3.Wkp_seminorm(1,1);

848

real W12_seminorm_U = sqrt( sqr(Ufld_1.Wkp_seminorm(1,2)) + sqr(Ufld_2.Wkp_seminorm(1,2)) + sqr(Ufld_3.Wkp_seminorm(1,2)) );

849

real W1max_seminorm_U = max( max(Ufld_1.Wkmax_seminorm(1),Ufld_2.Wkmax_seminorm(1)), Ufld_3.Wkmax_seminorm(1) );

850

851

real W11_norm_U = L1_norm_U + W11_seminorm_U;

852

real W12_norm_U = sqrt( sqr(L2_norm_U) + sqr(W12_seminorm_U) );

853

real W1max_norm_U = max( Lmax_norm_U, W1max_seminorm_U );

854

855

real W21_seminorm_U,W22_seminorm_U,W2max_seminorm_U,W21_norm_U,W22_norm_U,W2max_norm_U;

856

W21_seminorm_U = W22_seminorm_U = W2max_seminorm_U = W21_norm_U = W22_norm_U = W2max_norm_U = 0.0;

857

858

// W21_seminorm_U = Ufld_1.Wkp_seminorm(2,1) + Ufld_2.Wkp_seminorm(2,1) + Ufld_3.Wkp_seminorm(2,1);

859

// W22_seminorm_U = sqrt( sqr(Ufld_1.Wkp_seminorm(2,2)) + sqr(Ufld_2.Wkp_seminorm(2,2)) + sqr(Ufld_3.Wkp_seminorm(2,2)) );

860

// W2max_seminorm_U = max( max(Ufld_1.Wkmax_seminorm(2),Ufld_2.Wkmax_seminorm(2)), Ufld_3.Wkmax_seminorm(2) );

861

862

// W21_norm_U = W11_norm_U + W21_seminorm_U;

863

// W22_norm_U = sqrt( sqr(W12_norm_U) + sqr(W22_seminorm_U) );

864

// W2max_norm_U = max( W1max_norm_U, W2max_seminorm_U );

865

866

real L1_norm_P = Pfld.Lp_norm(1);

867

real L2_norm_P = Pfld.Lp_norm(2);

868

real Lmax_norm_P = Pfld.Lmax_norm();

869

870

real W11_seminorm_P = Pfld.Wkp_seminorm(1,1);

871

real W12_seminorm_P = Pfld.Wkp_seminorm(1,2);

872

real W1max_seminorm_P = Pfld.Wkmax_seminorm(1);

873

874

real W11_norm_P = L1_norm_P + W11_seminorm_P;

875

real W12_norm_P = sqrt( sqr(L2_norm_P) + sqr(W12_seminorm_P) );

876

real W1max_norm_P = max( Lmax_norm_P, W1max_seminorm_P );

877

878

real W21_seminorm_P,W22_seminorm_P,W2max_seminorm_P,W21_norm_P,W22_norm_P,W2max_norm_P;

879

W21_seminorm_P = W22_seminorm_P = W2max_seminorm_P = W21_norm_P = W22_norm_P = W2max_norm_P = 0.0;

880

881

// W21_seminorm_P = Pfld.Wkp_seminorm(2,1);

882

// W22_seminorm_P = Pfld.Wkp_seminorm(2,2);

883

// W2max_seminorm_P = Pfld.Wkmax_seminorm(2);

884

885

// W21_norm_P = W11_norm_P + W21_seminorm_P;

886

// W22_norm_P = sqrt( sqr(W12_norm_P) + sqr(W22_seminorm_P) );

887

// W2max_norm_P = max( W1max_norm_P, W2max_seminorm_P );

888

889

// Compute the divergence of the velocity field

890

MV_Vector<real> divU(no_cells);

891

for ( el=0; el < no_cells; el++ ) divU(el) = Ufld_1.EvalGradient(el,0) + Ufld_2.EvalGradient(el,1) + Ufld_3.EvalGradient(el,2);

892

893

Field Ufld_div(grid);

894

Ufld_div.InitByElms(&divU);

895

896

real L1_norm_divU = Ufld_div.Lp_norm(1);

897

real L2_norm_divU = Ufld_div.Lp_norm(2);

898

real Lmax_norm_divU = Ufld_div.Lmax_norm();

899

900

// Compute time derivative

901

MV_ColMat<real> Udot(no_nodes,(*U).size(1));

902

for ( nod=0; nod < no_nodes; nod++ ){

903

for ( i=0; i < (*U).size(1); i++ ) Udot(nod,i) = ( (*U)(nod,i) - (*Uprev)(nod,i) ) / dt;

904

}

905

Field Udotfld_1(grid);

906

Field Udotfld_2(grid);

907

Field Udotfld_3(grid);

908

Field Pdotfld(grid);

909

Udotfld_1.InitByNodes(&Udot,0);

910

Udotfld_2.InitByNodes(&Udot,1);

911

Udotfld_3.InitByNodes(&Udot,2);

912

Pdotfld.InitByNodes(&Udot,3);

913

914

real L1_norm_Udot = Udotfld_1.Lp_norm(1) + Udotfld_2.Lp_norm(1) + Udotfld_3.Lp_norm(1);

915

real L2_norm_Udot = sqrt( sqr(Udotfld_1.Lp_norm(2)) + sqr(Udotfld_2.Lp_norm(2)) + sqr(Udotfld_3.Lp_norm(2)) );

916

real Lmax_norm_Udot = max( max(Udotfld_1.Lmax_norm(),Udotfld_2.Lmax_norm()) , Udotfld_3.Lmax_norm() );

917

918

real L1_norm_Pdot = Pdotfld.Lp_norm(1);

919

real L2_norm_Pdot = Pdotfld.Lp_norm(2);

920

real Lmax_norm_Pdot = Pdotfld.Lmax_norm();

921

922

// Write norms of solution to matlab file

923

matlab_data->SetTime(t);

924

925

matlab_data->Set(0,L1_norm_U);

926

matlab_data->Set(1,L2_norm_U);

927

matlab_data->Set(2,Lmax_norm_U);

928

matlab_data->Set(3,W11_seminorm_U);

929

matlab_data->Set(4,W12_seminorm_U);

930

matlab_data->Set(5,W1max_seminorm_U);

931

matlab_data->Set(6,W11_norm_U);

932

matlab_data->Set(7,W12_norm_U);

933

matlab_data->Set(8,W1max_norm_U);

934

matlab_data->Set(9,W21_seminorm_U);

935

matlab_data->Set(10,W22_seminorm_U);

936

matlab_data->Set(11,W2max_seminorm_U);

937

matlab_data->Set(12,W21_norm_U);

938

matlab_data->Set(13,W22_norm_U);

939

matlab_data->Set(14,W2max_norm_U);

940

matlab_data->Set(15,L1_norm_P);

941

matlab_data->Set(16,L2_norm_P);

942

matlab_data->Set(17,Lmax_norm_P);

943

matlab_data->Set(18,W11_seminorm_P);

944

matlab_data->Set(19,W12_seminorm_P);

945

matlab_data->Set(20,W1max_seminorm_P);

946

matlab_data->Set(21,W11_norm_P);

947

matlab_data->Set(22,W12_norm_P);

948

matlab_data->Set(23,W1max_norm_P);

949

matlab_data->Set(24,W21_seminorm_P);

950

matlab_data->Set(25,W22_seminorm_P);

951

matlab_data->Set(26,W2max_seminorm_P);

952

matlab_data->Set(27,W21_norm_P);

953

matlab_data->Set(28,W22_norm_P);

954

matlab_data->Set(29,W2max_norm_P);

955

matlab_data->Set(30,L1_norm_divU);

956

matlab_data->Set(31,L2_norm_divU);

957

matlab_data->Set(32,Lmax_norm_divU);

958

matlab_data->Set(33,L1_norm_Udot);

959

matlab_data->Set(34,L2_norm_Udot);

960

matlab_data->Set(35,Lmax_norm_Udot);

961

matlab_data->Set(36,L1_norm_Pdot);

962

matlab_data->Set(37,L2_norm_Pdot);

963

matlab_data->Set(38,Lmax_norm_Pdot);

964

matlab_data->Save();

965

966

967

}

968

//----------------------------------------------------------------------------------------------------

969

void ProblemNS::CompContResidual()

970

{

971

972

int nod,el,sd;

973

974

// R1 = | u_t + u*Du + Dp - f |

975

// R2 = max jump in gradient over element boundaries divided by h

976

// R3 = time discontinuity (=0 for continuous time integration)

977

// R4 = | div u |

978

979

// Computing the R2 residual is expensive due to the approximation of second derivatives

980

bool compute_R2 = false;

981

982

//neighbor = &((*GRID)(compGrid).getNeighbor());

983

//if ( compute_R2 ){

984

// if ( neighbor->elementSize() == 0 ) neighbor->init( (*GRID)(compGrid), false, false, true );

985

// }

986

987

988

Grid *grid = grid;

989

990

Field Ufld_1(grid);

991

Field Ufld_2(grid);

992

Field Ufld_3(grid);

993

Field Pfld(grid);

994

995

Ufld_1.InitByNodes(U,0);

996

Ufld_2.InitByNodes(U,1);

997

Ufld_3.InitByNodes(U,2);

998

Pfld.InitByNodes(U,3);

999

1000

MV_ColMat<real> Udot(no_nodes,nsd);

1001

for ( nod=0; nod < no_nodes; nod++ ){

1002

for ( sd=0; sd < nsd; sd++ ) Udot(nod,sd) = ( (*U)(nod,sd) - (*Uprev)(nod,sd) ) / dt;

1003

}

1004

Field Udotfld_1(grid);

1005

Field Udotfld_2(grid);

1006

Field Udotfld_3(grid);

1007

Udotfld_1.InitByNodes(&Udot,0);

1008

Udotfld_2.InitByNodes(&Udot,1);

1009

Udotfld_3.InitByNodes(&Udot,2);

1010

1011

// R3=0 for continuous time integration

1012

MV_ColMat<real> R1(no_cells,nsd);

1013

MV_Vector<real> R2(no_cells);

1014

//MV_ColMat<real> R3(no_cells,nsd);

1015

MV_Vector<real> R4(no_cells);

1016

1017

Field Ffld_1(grid);

1018

Field Ffld_2(grid);

1019

Field Ffld_3(grid);

1020

real zero = 0.0;

1021

Ffld_1.InitByConstant(&zero);

1022

Ffld_2.InitByConstant(&zero);

1023

Ffld_3.InitByConstant(&zero);

1024

1025

for ( el=0; el < no_cells; el++ ){

1026

R1(el,0) = fabs( Udotfld_1.Eval(el) + Ufld_1.Eval(el)*Ufld_1.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_1.EvalGradient(el,1)

1027

+ Ufld_3.Eval(el)*Ufld_1.EvalGradient(el,2) + Pfld.EvalGradient(el,0) - Ffld_1.Eval(el) );

1028

R1(el,1) = fabs( Udotfld_2.Eval(el) + Ufld_1.Eval(el)*Ufld_2.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_2.EvalGradient(el,1)

1029

+ Ufld_3.Eval(el)*Ufld_2.EvalGradient(el,2) + Pfld.EvalGradient(el,1) - Ffld_2.Eval(el) );

1030

R1(el,2) = fabs( Udotfld_3.Eval(el) + Ufld_1.Eval(el)*Ufld_3.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_3.EvalGradient(el,1)

1031

+ Ufld_3.Eval(el)*Ufld_3.EvalGradient(el,2) + Pfld.EvalGradient(el,2) - Ffld_3.Eval(el) );

1032

}

1033

Field R1fld_1(grid);

1034

Field R1fld_2(grid);

1035

Field R1fld_3(grid);

1036

R1fld_1.InitByElms(&R1,0);

1037

R1fld_2.InitByElms(&R1,1);

1038

R1fld_3.InitByElms(&R1,2);

1039

1040

// Compute norms of R1 residual

1041

real L1_norm_R1 = R1fld_1.Lp_norm(1) + R1fld_2.Lp_norm(1) + R1fld_3.Lp_norm(1);

1042

real L2_norm_R1 = sqrt( sqr(R1fld_1.Lp_norm(2)) + sqr(R1fld_2.Lp_norm(2)) + sqr(R1fld_3.Lp_norm(2)) );

1043

real Lmax_norm_R1 = max( max(R1fld_1.Lmax_norm(),R1fld_2.Lmax_norm()) , R1fld_3.Lmax_norm() );

1044

1045

// Compute norms of R2 residual

1046

real L1_norm_R2,L2_norm_R2,Lmax_norm_R2;

1047

if ( compute_R2 ){

1048

1049

MV_Vector<real> locGrad(nsd);

1050

MV_Vector<real> neighborGrad(nsd);

1051

1052

int i,el_start,el_stop,noNeighbors;

1053

Cell *c;

1054

1055

for ( el=0; el < no_cells; el++ ){

1056

1057

c = grid->GetCell(el);

1058

1059

R2(el) = 0.0;

1060

Ufld_1.EvalGradient(el,locGrad);

1061

1062

for (int i=0;i<c->GetNoCellNeighbors();i++){

1063

Ufld_1.EvalGradient(c->GetCellNeighbor(i),neighborGrad);

1064

for (sd=0;sd<nsd;sd++)

1065

if ( fabs(locGrad(sd) - neighborGrad(sd)) > R2(el) )

1066

R2(el) = fabs(locGrad(sd) - neighborGrad(sd));

1067

}

1068

1069

Ufld_2.EvalGradient(el,locGrad);

1070

1071

for (int i=0;i<c->GetNoCellNeighbors();i++){

1072

Ufld_2.EvalGradient(c->GetCellNeighbor(i),neighborGrad);

1073

for ( sd=0; sd < nsd; sd++ ){

1074

if ( fabs(locGrad(sd) - neighborGrad(sd)) > R2(el) ) R2(el) = fabs(locGrad(sd) - neighborGrad(sd));

1075

}

1076

}

1077

1078

Ufld_3.EvalGradient(el,locGrad);

1079

1080

for (int i=0;i<c->GetNoCellNeighbors();i++){

1081

Ufld_3.EvalGradient(c->GetCellNeighbor(i),neighborGrad);

1082

for ( sd=0; sd < nsd; sd++ ){

1083

if ( fabs(locGrad(sd) - neighborGrad(sd)) > R2(el) ) R2(el) = fabs(locGrad(sd) - neighborGrad(sd));

1084

}

1085

}

1086

1087

R2(el) /= ( reynolds_number * LocalMeshSize(el) );

1088

}

1089

Field R2fld(grid);

1090

R2fld.InitByElms(&R2);

1091

1092

L1_norm_R2 = R2fld.Lp_norm(1);

1093

L2_norm_R2 = R2fld.Lp_norm(2);

1094

Lmax_norm_R2 = R2fld.Lmax_norm();

1095

}

1096

1097

// Compute norms of R4 residual

1098

for ( el=0; el < no_cells; el++ ) R4(el) = Ufld_1.EvalGradient(el,0) + Ufld_2.EvalGradient(el,1) + Ufld_3.EvalGradient(el,2);

1099

Field R4fld(grid);

1100

R4fld.InitByElms(&R4);

1101

1102

real L1_norm_R4 = R4fld.Lp_norm(1);

1103

real L2_norm_R4 = R4fld.Lp_norm(2);

1104

real Lmax_norm_R4 = R4fld.Lmax_norm();

1105

1106

if ( !compute_R2 ){

1107

L1_norm_R2 = 0.0;

1108

L2_norm_R2 = 0.0;

1109

Lmax_norm_R2 = 0.0;

1110

}

1111

real L1_norm_R3 = 0.0;

1112

real L2_norm_R3 = 0.0;

1113

real Lmax_norm_R3 = 0.0;

1114

1115

// Write norms to matlab file

1116

matlab_data->SetTime(t);

1117

1118

matlab_data->Set(no_data_entries+0,L1_norm_R1);

1119

matlab_data->Set(no_data_entries+1,L2_norm_R1);

1120

matlab_data->Set(no_data_entries+2,Lmax_norm_R1);

1121

matlab_data->Set(no_data_entries+3,L1_norm_R2);

1122

matlab_data->Set(no_data_entries+4,L2_norm_R2);

1123

matlab_data->Set(no_data_entries+5,Lmax_norm_R2);

1124

matlab_data->Set(no_data_entries+6,L1_norm_R3);

1125

matlab_data->Set(no_data_entries+7,L2_norm_R3);

1126

matlab_data->Set(no_data_entries+8,Lmax_norm_R3);

1127

matlab_data->Set(no_data_entries+9,L1_norm_R4);

1128

matlab_data->Set(no_data_entries+10,L2_norm_R4);

1129

matlab_data->Set(no_data_entries+11,Lmax_norm_R4);

1130

matlab_data->Save();

1131

1132

// Write element representation of residual

1133

contResidual.newsize(no_cells);

1134

for ( el=0; el < no_cells; el++ ) contResidual(el) = R1(el,0) + R1(el,1) + R1(el,2) + R4(el);

1135

if ( compute_R2 ){

1136

for ( el=0; el < no_cells; el++ ) contResidual(el) += R2(el);

1137

}

1138

1139

1140

// Write nodal representation of residual

1141

//MV_Vector<real> R1_1_vec;

1142

//MV_Vector<real> R1_2_vec;

1143

//MV_Vector<real> R1_3_vec;

1144

//MV_Vector<real> R2_vec;

1145

//MV_Vector<real> R4_vec;

1146

//R1fld_1.getNodalVector(R1_1_vec);

1147

//R1fld_2.getNodalVector(R1_2_vec);

1148

//R1fld_3.getNodalVector(R1_3_vec);

1149

//R2fld.getNodalVector(R2_vec);

1150

//R4fld.getNodalVector(R4_vec);

1151

//MV_ColMat<real> printRes(printNoNodes,5);

1152

//for ( i=0; i < printNoNodes; i++ ){

1153

//printRes(i,0) = R1_1_vec(i);

1154

//printRes(i,1) = R1_2_vec(i);

1155

//printRes(i,2) = R1_3_vec(i);

1156

//printRes(i,3) = R2_vec(i);

1157

//printRes(i,4) = R4_vec(i);

1158

//}

1159

//opendx_residuals->AddFrame(&(*GRID)(printGrid),&printRes,t);

1160

1161

}

1162

//----------------------------------------------------------------------------------------------------

1163

void ProblemNS::WritePertubationData()

1164

{

1165

1166

int nod,el;

1167

real x,y,z;

1168

Point p;

1169

Grid *grid = grid;

1170

1171

// Write pertubation data to file

1172

Field Ufld_1(grid);

1173

Field Ufld_2(grid);

1174

Field Ufld_3(grid);

1175

1176

Ufld_1.InitByNodes(U,0);

1177

Ufld_2.InitByNodes(U,1);

1178

Ufld_3.InitByNodes(U,2);

1179

1180

Field Upertfld_1(grid);

1181

Field Upertfld_2(grid);

1182

Field Upertfld_3(grid);

1183

1184

MV_Vector<real> Utmp(no_nodes);

1185

1186

for (int i=0;i<no_nodes;i++){

1187

p = grid->GetNode(i)->GetCoord();

1188

y = real(p.y);

1189

z = real(p.z);

1190

if ( pert_poiseuille ) Utmp(nod) = (*U)(nod,0) - ( 16.0*y*(1.0-y)*z*(1-z) );

1191

if ( pert_couette ) Utmp(nod) = (*U)(nod,0) - ( 2.0*y-1.0 );

1192

}

1193

Upertfld_1.InitByNodes(&Utmp);

1194

Upertfld_2.InitByNodes(U,1);

1195

Upertfld_3.InitByNodes(U,2);

1196

1197

real L2_norm_U1pert, L2_norm_U2pert, L2_norm_U3pert;

1198

real Lmax_norm_U1pert, Lmax_norm_U2pert, Lmax_norm_U3pert;

1199

real L2_norm_D1U1, L2_norm_D2U1, L2_norm_D3U1;

1200

real L2_norm_D1U2, L2_norm_D2U2, L2_norm_D3U2;

1201

real L2_norm_D1U3, L2_norm_D2U3, L2_norm_D3U3;

1202

real Lmax_norm_D1U1, Lmax_norm_D2U1, Lmax_norm_D3U1;

1203

real Lmax_norm_D1U2, Lmax_norm_D2U2, Lmax_norm_D3U2;

1204

real Lmax_norm_D1U3, Lmax_norm_D2U3, Lmax_norm_D3U3;

1205

1206

L2_norm_U1pert = Upertfld_1.Lp_norm(2);

1207

L2_norm_U2pert = Upertfld_2.Lp_norm(2);

1208

L2_norm_U3pert = Upertfld_3.Lp_norm(2);

1209

Lmax_norm_U1pert = Upertfld_1.Lmax_norm();

1210

Lmax_norm_U2pert = Upertfld_2.Lmax_norm();

1211

Lmax_norm_U3pert = Upertfld_3.Lmax_norm();

1212

1213

Field Ugradfld(grid);

1214

Utmp.newsize(no_cells);

1215

1216

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_1.EvalGradient(el,0);

1217

Ugradfld.InitByElms(&Utmp);

1218

L2_norm_D1U1 = Ugradfld.Lp_norm(2);

1219

Lmax_norm_D1U1 = Ugradfld.Lmax_norm();

1220

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_1.EvalGradient(el,1);

1221

Ugradfld.InitByElms(&Utmp);

1222

L2_norm_D2U1 = Ugradfld.Lp_norm(2);

1223

Lmax_norm_D2U1 = Ugradfld.Lmax_norm();

1224

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_1.EvalGradient(el,2);

1225

Ugradfld.InitByElms(&Utmp);

1226

L2_norm_D3U1 = Ugradfld.Lp_norm(2);

1227

Lmax_norm_D3U1 = Ugradfld.Lmax_norm();

1228

1229

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_2.EvalGradient(el,0);

1230

Ugradfld.InitByElms(&Utmp);

1231

L2_norm_D1U2 = Ugradfld.Lp_norm(2);

1232

Lmax_norm_D1U2 = Ugradfld.Lmax_norm();

1233

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_2.EvalGradient(el,1);

1234

Ugradfld.InitByElms(&Utmp);

1235

L2_norm_D2U2 = Ugradfld.Lp_norm(2);

1236

Lmax_norm_D2U2 = Ugradfld.Lmax_norm();

1237

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_2.EvalGradient(el,2);

1238

Ugradfld.InitByElms(&Utmp);

1239

L2_norm_D3U2 = Ugradfld.Lp_norm(2);

1240

Lmax_norm_D3U2 = Ugradfld.Lmax_norm();

1241

1242

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_3.EvalGradient(el,0);

1243

Ugradfld.InitByElms(&Utmp);

1244

L2_norm_D1U3 = Ugradfld.Lp_norm(2);

1245

Lmax_norm_D1U3 = Ugradfld.Lmax_norm();

1246

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_3.EvalGradient(el,1);

1247

Ugradfld.InitByElms(&Utmp);

1248

L2_norm_D2U3 = Ugradfld.Lp_norm(2);

1249

Lmax_norm_D2U3 = Ugradfld.Lmax_norm();

1250

for ( el=0; el < no_cells; el++ ) Utmp(el) = Ufld_3.EvalGradient(el,2);

1251

Ugradfld.InitByElms(&Utmp);

1252

L2_norm_D3U3 = Ugradfld.Lp_norm(2);

1253

Lmax_norm_D3U3 = Ugradfld.Lmax_norm();

1254

1255

// Write data to matlab file

1256

matlab_data->SetTime(t);

1257

1258

matlab_data->Set(no_data_entries+no_residual_entries+0,L2_norm_U1pert);

1259

matlab_data->Set(no_data_entries+no_residual_entries+1,Lmax_norm_U1pert);

1260

matlab_data->Set(no_data_entries+no_residual_entries+2,L2_norm_U2pert);

1261

matlab_data->Set(no_data_entries+no_residual_entries+3,Lmax_norm_U2pert);

1262

matlab_data->Set(no_data_entries+no_residual_entries+4,L2_norm_U3pert);

1263

matlab_data->Set(no_data_entries+no_residual_entries+5,Lmax_norm_U3pert);

1264

matlab_data->Set(no_data_entries+no_residual_entries+6,L2_norm_D1U1);

1265

matlab_data->Set(no_data_entries+no_residual_entries+7,Lmax_norm_D1U1);

1266

matlab_data->Set(no_data_entries+no_residual_entries+8,L2_norm_D2U1);

1267

matlab_data->Set(no_data_entries+no_residual_entries+9,Lmax_norm_D2U1);

1268

matlab_data->Set(no_data_entries+no_residual_entries+10,L2_norm_D3U1);

1269

matlab_data->Set(no_data_entries+no_residual_entries+11,Lmax_norm_D3U1);

1270

matlab_data->Set(no_data_entries+no_residual_entries+12,L2_norm_D1U2);

1271

matlab_data->Set(no_data_entries+no_residual_entries+13,Lmax_norm_D1U2);

1272

matlab_data->Set(no_data_entries+no_residual_entries+14,L2_norm_D2U2);

1273

matlab_data->Set(no_data_entries+no_residual_entries+15,Lmax_norm_D2U2);

1274

matlab_data->Set(no_data_entries+no_residual_entries+16,L2_norm_D3U2);

1275

matlab_data->Set(no_data_entries+no_residual_entries+17,Lmax_norm_D3U2);

1276

matlab_data->Set(no_data_entries+no_residual_entries+18,L2_norm_D1U3);

1277

matlab_data->Set(no_data_entries+no_residual_entries+19,Lmax_norm_D1U3);

1278

matlab_data->Set(no_data_entries+no_residual_entries+20,L2_norm_D2U3);

1279

matlab_data->Set(no_data_entries+no_residual_entries+21,Lmax_norm_D2U3);

1280

matlab_data->Set(no_data_entries+no_residual_entries+22,L2_norm_D3U3);

1281

matlab_data->Set(no_data_entries+no_residual_entries+23,Lmax_norm_D3U3);

1282

matlab_data->Save();

1283

1284

}

1285

//----------------------------------------------------------------------------------------------------

1286

void ProblemNS::WriteReynoldsStresses()

1287

{

1288

}

1289

//----------------------------------------------------------------------------------------------------

1290

void ProblemNS::ComputeFunctionals()

1291

{

1292

1293

int nod,el,sd;

1294

1295

Field Ufld_1(grid);

1296

Field Ufld_2(grid);

1297

Field Ufld_3(grid);

1298

Field Pfld(grid);

1299

1300

Ufld_1.InitByNodes(U,0);

1301

Ufld_2.InitByNodes(U,1);

1302

Ufld_3.InitByNodes(U,2);

1303

Pfld.InitByNodes(U,3);

1304

1305

MV_ColMat<real> Udot(no_nodes,nsd);

1306

for ( nod=0; nod < no_nodes; nod++ ){

1307

for ( sd=0; sd < nsd; sd++ ) Udot(nod,sd) = ( (*U)(nod,sd) - (*Uprev)(nod,sd) ) / dt;

1308

}

1309

Field Udotfld_1(grid);

1310

Field Udotfld_2(grid);

1311

Field Udotfld_3(grid);

1312

Udotfld_1.InitByNodes(&Udot,0);

1313

Udotfld_2.InitByNodes(&Udot,1);

1314

Udotfld_3.InitByNodes(&Udot,2);

1315

1316

Field Ffld_1(grid);

1317

Field Ffld_2(grid);

1318

Field Ffld_3(grid);

1319

real zero = 0.0;

1320

Ffld_1.InitByConstant(&zero);

1321

Ffld_2.InitByConstant(&zero);

1322

Ffld_3.InitByConstant(&zero);

1323

1324

Field Psi_1(grid);

1325

Field Psi_2(grid);

1326

Field Psi_3(grid);

1327

Psi_1.InitByConstant(&zero);

1328

Psi_2.InitByConstant(&zero);

1329

Psi_3.InitByConstant(&zero);

1330

1331

real drag = 0.0;

1332

for ( el=0; el < no_cells; el++ ){

1333

drag += ( Udotfld_1.Eval(el) + Ufld_1.Eval(el)*Ufld_1.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_1.EvalGradient(el,1)

1334

+ Ufld_3.Eval(el)*Ufld_1.EvalGradient(el,2) ) * Psi_1.Eval(el)

1335

+ ( Udotfld_2.Eval(el) + Ufld_1.Eval(el)*Ufld_2.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_2.EvalGradient(el,1)

1336

+ Ufld_3.Eval(el)*Ufld_2.EvalGradient(el,2) ) * Psi_2.Eval(el)

1337

+ ( Udotfld_3.Eval(el) + Ufld_1.Eval(el)*Ufld_3.EvalGradient(el,0) + Ufld_2.Eval(el)*Ufld_3.EvalGradient(el,1)

1338

+ Ufld_3.Eval(el)*Ufld_3.EvalGradient(el,2) ) * Psi_3.Eval(el)

1339

- Pfld.Eval(el) * (Psi_1.EvalGradient(el,0)+Psi_2.EvalGradient(el,1)+Psi_3.EvalGradient(el,2))

1340

+ (2.0/reynolds_number) * ( Ufld_1.EvalGradient(el,0)*Psi_1.EvalGradient(el,0)

1341

+ 0.5*(Ufld_1.EvalGradient(el,1)+Ufld_2.EvalGradient(el,0))*0.5*(Psi_1.EvalGradient(el,1)+Psi_2.EvalGradient(el,0))

1342

+ 0.5*(Ufld_1.EvalGradient(el,2)+Ufld_3.EvalGradient(el,0))*0.5*(Psi_1.EvalGradient(el,2)+Psi_3.EvalGradient(el,0))

1343

+ 0.5*(Ufld_2.EvalGradient(el,0)+Ufld_1.EvalGradient(el,1))*0.5*(Psi_2.EvalGradient(el,0)+Psi_1.EvalGradient(el,1))

1344

+ Ufld_2.EvalGradient(el,1)*Psi_2.EvalGradient(el,1)

1345

+ 0.5*(Ufld_2.EvalGradient(el,2)+Ufld_3.EvalGradient(el,1))*0.5*(Psi_2.EvalGradient(el,2)+Psi_3.EvalGradient(el,1))

1346

+ 0.5*(Ufld_3.EvalGradient(el,0)+Ufld_1.EvalGradient(el,2))*0.5*(Psi_3.EvalGradient(el,0)+Psi_1.EvalGradient(el,0))

1347

+ 0.5*(Ufld_3.EvalGradient(el,1)+Ufld_2.EvalGradient(el,2))*0.5*(Psi_3.EvalGradient(el,1)+Psi_2.EvalGradient(el,0))

1348

+ Ufld_3.EvalGradient(el,2)*Psi_3.EvalGradient(el,2) )

1349

- ( Ffld_1.Eval(el)*Psi_1.Eval(el) + Ffld_2.Eval(el)*Psi_2.Eval(el) + Ffld_3.Eval(el)*Psi_3.Eval(el) );

1350

}

1351

1352

cout << "Drag force = " << drag << endl;

1353

1354

1355

// Write norms to matlab file

1356

//matlab_data->SetTime(t);

1357

//matlab_data->Set(no_data_entries+0,drag);

1358

//matlab_data->Save();

1359

1360

}

1361

//----------------------------------------------------------------------------------------------------

1362

void ProblemNS::SetTolerances()

1363

{

1364

// Set tolerances for iterative methods

1365

NL_tol = 1.0e-4;

1366

UP_tol = 1.0e-1;

1367

stat_tol = 1.0e-10;

1368

1369

NL_max_tol = 1000.0;

1370

UP_max_tol = 1000.0;

1371

1372

max_no_UP_iter = 3;

1373

max_no_NL_iter = 100;

1374

1375

UP_conv = true;

1376

NL_conv = true;

1377

}

1378

//-----------------------------------------------------------------------------

1379

real ProblemNS::GetTimeStep()

1380

{

1381

return grid_h;

1382

}

1383

//-----------------------------------------------------------------------------

1384

void ProblemNS::GetMeshSize()

1385

{

1386

// Compute smallest element diameter in grid

1387

1388

real lms = LocalMeshSize(0);

1389

grid_h = lms;

1390

for (int el = 1; el < no_cells; el++) {

1391

lms = LocalMeshSize(el);

1392

if (lms < grid_h) grid_h = lms;

1393

}

1394

}

1395

//-----------------------------------------------------------------------------

1396

real ProblemNS::LocalMeshSize(int el)

1397

{

1398

// Compute element diameter

1399

1400

real gh = 2.0 * grid->GetCell(el)->ComputeCircumRadius(grid);

1401

1402

return gh;

1403

}

1404

//----------------------------------------------------------------------------------------------------

1405

void ProblemNS::SetTurbulentInflow( real trb )

1406

{

1407

1408

int ierr;

1409

int i,j;

1410

int *bndNodes;

1411

int *b_nodes;

1412

int noBndNodes;

1413

1414

b_nodes = new int[no_nodes];

1415

1416

ierr = getBndNodes(b_nodes, &noBndNodes, bndx0, grid);

1417

1418

turbulence.newsize(noBndNodes,noeq);

1419

for ( i=0; i < noBndNodes; i++ ){

1420

for ( j=0; j < noeq; j++ ){

1421

turbulence(i,j) = trb * 2.0*(drand48() - 0.5);

1422

}

1423

}

1424

1425

delete b_nodes

1426

1427

}

1428

//----------------------------------------------------------------------------------------------------

Older »