1 function sim_data = ldg_conv_detailed_simulation(model,model_data,dummy)
2 %
function sim_data = ldg_conv_detailed_simulation(model,model_data,dummy)
4 % required fields of model:
6 % nt : number of time-intervals until T, i.e. nt+1
7 % solution slices are computed
8 % init_values_algorithm: name of function for computing the
9 % initvalues-DOF with arguments (grid, params)
10 % example: init_values_cog
11 % operators_algorithm: name of function for computing the
12 % L_E,L_I,b-operators with arguments (grid, params)
13 % example operators_conv_diff
14 % data_const_in_time : if this optional field is 1, the time
15 % evolution is performed with constant operators,
16 % i.e. only the initial-time-operators are computed
17 % and used throughout the time simulation.
19 % fields of sim_data: U an array of columnwise degrees of freedom-vectors
21 % Bernard Haasdonk 27.8.2009
23 disp(' to be adjusted!! ');
27 disp([
'entered detailed simulation ']);
30 if isempty(model_data)
31 model_data = gen_model_data(model);
34 model.decomp_mode = 0; % == complete;
37 % initial values by midpoint evaluation
38 U0 = model.init_values_algorithm(model,model_data);
44 U = zeros(length(U0),model.nt+1);
46 model.dt = model.T/model.nt;
47 %operators_required = 1;
48 %if ~isfield(model,'data_const_in_time')
49 % model.data_const_in_time = 0;
53 params.axis_equal = 1;
54 model.plot(U0,model_data.grid,params);
57 % loop over time-steps
61 disp(['entered time-loop step ',num2str(t)]);
66 % get matrices and bias-vector
69 model.t = (t-1)*model.dt;
70 [L_I, L_E, b] = model.operators_ptr(model, model_data);
71 if model.data_const_in_time
72 operators_required = 0;
76 rhs = L_E * U{t}.dofs + b;
77 if isequal(L_I, speye(size(L_I)))
81 % disp(
'check symmetry and choose solver accordingly!');
83 % nonsymmetric solvers:
84 % [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
85 % [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
87 % symmetric solver, non pd:
89 % see bug_symmlq.mat
for a very strange bug: cannot solve identity system!
90 % reported to matlab central, but no solution up to now.
91 % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
93 % [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
94 % symmetric solver, pd:
95 %[U(:,t+1), flag] = pcg(L_I,rhs,[],1000);
96 % bicgstab works also quite well:
97 %[U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
98 U{t+1}.dofs = L_I \ rhs;
101 % disp([
'error in system solution, solver return flag = ', ...
108 if model.verbose == 9