Skip to content
Snippets Groups Projects
Commit cd6a9249 authored by Stephan Rave's avatar Stephan Rave
Browse files

add rb-intro sheet 3

parent 5f1b44df
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
# pyMOR RB-Intro
# Exercise Sheet 3 - Solutions
%% Cell type:code id: tags:
``` python
from pymor.basic import *
set_log_levels({'pymor': 'WARN',
'pymor.functions.bitmap': 'CRITICAL'})
```
%% Cell type:markdown id: tags:
We will continue to work with the models defined on exercise sheet 1. To simplify working with these models, you can find their definition in `problems.py` in the same directory as this notebook:
%% Cell type:code id: tags:
``` python
%cat problems.py
```
%% Output
from pymor.basic import *
def problem_1_1_a():
return StationaryProblem(
domain=RectDomain([[0., 0.], [1., 1.]]),
diffusion=ConstantFunction(1, 2),
rhs=ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),
name='Sheet1_Problem_1a'
)
def problem_1_1_b():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),
name='Sheet1_Problem_1b'
)
def problem_1_1_c():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
),
name='Sheet1_Problem_1c'
)
def problem_1_1_d():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=BitmapFunction('RB.png', range=[0.001, 1]),
name='Sheet1_Problem_1d'
)
def problem_1_2_a():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
),
parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),
name='Sheet1_Problem_2a'
)
def problem_1_2_b():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
2, (),
values={'K': 10},
parameter_type={'diffu': ()}
),
parameter_space=CubicParameterSpace(
{'neum': (), 'diffu': ()},
ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}
),
name='Sheet1_Problem_2b'
)
def problem_1_2_c():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=LincombFunction(
[ConstantFunction(1., 2),
BitmapFunction('R.png', range=[1, 0]),
BitmapFunction('B.png', range=[1, 0])],
[1.,
ExpressionParameterFunctional('-(1 - R)', {'R': ()}),
ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
),
parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),
name='Sheet1_Problem_2c'
)
%% Cell type:markdown id: tags:
The file can be imported as usual:
%% Cell type:code id: tags:
``` python
import problems
```
%% Cell type:markdown id: tags:
We have four non-parametric problems:
%% Cell type:code id: tags:
``` python
for f in [problems.problem_1_1_a, problems.problem_1_1_b,
problems.problem_1_1_c, problems.problem_1_1_d]:
p = f()
m, _ = discretize_stationary_cg(p, diameter=1/100)
m.visualize(m.solve(), legend=m.name)
```
%% Cell type:markdown id: tags:
The parametric problems have already been assigned a `ParameterSpace` which is inherited by the generated models:
%% Cell type:code id: tags:
``` python
for f in [problems.problem_1_2_a, problems.problem_1_2_b,
problems.problem_1_2_c]:
p = f()
m, _ = discretize_stationary_cg(p, diameter=1/100)
mu = m.parameter_space.sample_randomly(1)[0]
m.visualize(m.solve(mu), legend=str(mu))
```
%% Cell type:markdown id: tags:
## Problem 1
%% Cell type:markdown id: tags:
## Problem 1 a)
%% Cell type:markdown id: tags:
Write a method
```python
def orthogonal_projection(U, basis, product=None,
basis_is_orthonormal=False,
return_projection_matrix=False):
...
```
that encapsulates the orthogonal projection algorithm you implemented in Sheet 2, Problem 2.
%% Cell type:markdown id: tags:
## Problem 1 b)
%% Cell type:markdown id: tags:
Repeat the error calculations from Sheet 2, Problem 2 with an orthonormal basis obtained from `U`
using `pymor.algorithms.gram_schmidt`. Compare your results.
%% Cell type:markdown id: tags:
## Problem 2
%% Cell type:markdown id: tags:
## Problem 2 a)
%% Cell type:markdown id: tags:
Extend the `greedy` function from Sheet 2, Problem 3 to a function
```python
def greedy(U, basis_size, product=None,
rtol=0,
orthonormalize=False):
...
```
The algorithm should stop, when the relative approximation error goes below the prescribed value `rtol`. When `orthonormalize` is set to `True`, the algorithm shall return an orthonormal basis w.r.t. to `product` by orthonormalizing each new snapshot vector w.r.t. `product`.
**Hint:** When using `gram_schmidt` for orthonormalization, the `offset` parameter can be used to avoid re-orthonormalizing the old basis.
%% Cell type:markdown id: tags:
## Problem 2 b)
%% Cell type:markdown id: tags:
Repeat the experiments from Sheet 2, Problem 3, using `strong_greedy` with `orthonormalize=True`.
%% Cell type:markdown id: tags:
## Bonus Problem
%% Cell type:markdown id: tags:
Write a function `word_problem` which, given an arbitrary string, generates a `StationaryProblem` similar to Sheet 1, Problem 2 (c). Each letter should have a corresponding parameter component.
**Hint:** Use the modules `PIL.ImageDraw` und `PIL.ImageFont`.
from pymor.basic import *
def problem_1_1_a():
return StationaryProblem(
domain=RectDomain([[0., 0.], [1., 1.]]),
diffusion=ConstantFunction(1, 2),
rhs=ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),
name='Sheet1_Problem_1a'
)
def problem_1_1_b():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),
name='Sheet1_Problem_1b'
)
def problem_1_1_c():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
),
name='Sheet1_Problem_1c'
)
def problem_1_1_d():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=BitmapFunction('RB.png', range=[0.001, 1]),
name='Sheet1_Problem_1d'
)
def problem_1_2_a():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
),
parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),
name='Sheet1_Problem_2a'
)
def problem_1_2_b():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),
diffusion=ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
2, (),
values={'K': 10},
parameter_type={'diffu': ()}
),
parameter_space=CubicParameterSpace(
{'neum': (), 'diffu': ()},
ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}
),
name='Sheet1_Problem_2b'
)
def problem_1_2_c():
return StationaryProblem(
domain=RectDomain(bottom='neumann'),
neumann_data=ConstantFunction(-1., 2),
diffusion=LincombFunction(
[ConstantFunction(1., 2),
BitmapFunction('R.png', range=[1, 0]),
BitmapFunction('B.png', range=[1, 0])],
[1.,
ExpressionParameterFunctional('-(1 - R)', {'R': ()}),
ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
),
parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),
name='Sheet1_Problem_2c'
)
This diff is collapsed.
../problems.py
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment