properly test fixed h, refined H

parent d66de049
......@@ -2,13 +2,17 @@ import yaml
yaml.warnings({'YAMLLoadWarning': False})
import os
import numpy as np
np.warnings.filterwarnings('error')
np.warnings.filterwarnings('ignore', message="numpy.ufunc size changed")
np.warnings.filterwarnings('ignore', message="Using or importing the ABCs from 'collections'")
np.warnings.filterwarnings('ignore', message="The 'warn' method is deprecated, use 'warning' instead")
from simdb.run import new_dataset, append_values
from pymor.core.logger import set_log_levels, getLogger
from pymor.operators.block import BlockOperator
from thermalblock_problem import domain, init_problem
from dune.gdt.usercode import DomainDecomposition, ContinuousLagrangePartitionOfUnity
from dune.gdt.usercode import DomainDecomposition, ContinuousLagrangePartitionOfUnity, interpolate_pou
from distributed_elliptic_discretizer import discretize_elliptic_block_swipdg
from distributed_elliptic_reductor import PouLocalizedCoerciveReductor as Reductor
......@@ -18,6 +22,7 @@ set_log_levels({
FILENAME: 'INFO',
'distributed_elliptic_discretizer': 'INFO',
'pymor.discretizations': 'WARN',
'pymor.algorithms.gram_schmidt': 'WARN',
})
......@@ -177,6 +182,105 @@ def test_fixed_h_refined_H():
assert_results(run_fixed_h_refined_H(**kwargs), expected_results_fixed_h_refined_H(**kwargs))
def run_fixed_h_refined_H_with_fine_PoU_in_basis(
global_grid_refines=3,
num_H_refs_less=0,
initial_macro_grid=[1, 1],
make_pou=lambda dd: ContinuousLagrangePartitionOfUnity(dd),
pou_id='Q1PoU',
DATASET_ID=FILENAME + '.run_fixed_h_refined_H_with_fine_PoU_in_basis'):
config = dict(base_config)
config['local_space_type'] = 'dg_p1'
macro_grid = list(initial_macro_grid)
config['initial_macro_grid'] = initial_macro_grid
config['global_grid_refines'] = global_grid_refines
new_dataset(DATASET_ID, **config)
results = {'num_subdomains': [], 'num_DoFs': [], 'error': [], 'estimate': []}
# compute the finest PoU
finest_dd = DomainDecomposition([macro_grid[0]*(2**(global_grid_refines - num_H_refs_less)),
macro_grid[1]*(2**(global_grid_refines - num_H_refs_less))],
num_H_refs_less,
domain[0], domain[1])
finest_PoU = make_pou(finest_dd)
def compute_on_single_refinement(cfg):
logger = getLogger(f'{FILENAME}.compute_on_single_refinement')
def init_dd():
return DomainDecomposition(cfg['num_macro_elements'],
cfg['num_refinements_per_subdomain'],
domain[0],
domain[1])
dd = init_dd()
p = make_problem()
num_subdomains = dd.num_subdomains
fom = discretize_elliptic_block_swipdg(make_problem, init_dd, cfg['local_space_type'])
fom.globalize() # We need this later on, only here to trigger the output
num_DoFs = fom.solution_space.dim
append_values(num_DoFs=num_DoFs)
logger.info('reducing with coarse and fine PoU ...')
reductor = Reductor(
fom, init_dd, cfg, make_pou=make_pou, kappa_lower=init_problem()['kappa_lower'], pou_id=pou_id)
for ss in np.arange(dd.num_subdomains):
interpolated_fine_pou = interpolate_pou(
finest_PoU, dd, ss, cfg['local_space_type'], clean_up_since_grids_match=True)
reductor.extend_basis(ss, fom.local_spaces[ss].from_data(interpolated_fine_pou))
append_values(reductor=str(reductor))
append_values(basis_size=[len(b) for b in reductor.bases])
rom = reductor.reduce()
logger.info('building error product ...')
product = np.full((dd.num_subdomains, dd.num_subdomains), None)
for (ii, jj), local_prod in fom.products['broken_h1'].items():
product[ii, jj] = local_prod
product = BlockOperator(product)
logger.info('estimating ...')
mu = rom.parse_parameter(cfg['mu'])
u = rom.solve(mu=mu)
estimate = rom.estimate(u, mu=mu)[0]
logger.info('... done (is {})'.format(estimate))
append_values(estimate=estimate)
logger.info('computing error ...')
diff = fom.solve(mu) - reductor.reconstruct(u)
error = np.sqrt(product.apply2(diff, diff))[0][0]
logger.info('... done (is {})'.format(error))
append_values(error=error)
return num_subdomains, num_DoFs, error, estimate
for _ in np.arange(global_grid_refines - num_H_refs_less + 1):
cfg = dict(config)
cfg['num_macro_elements'] = macro_grid
cfg['num_refinements_per_subdomain'] = global_grid_refines
nsd, nD, er, es = compute_on_single_refinement(cfg)
results['num_subdomains'].append(nsd)
results['num_DoFs'].append(nD)
results['error'].append(er)
results['estimate'].append(es)
print('')
global_grid_refines -= 1
macro_grid[0] *= 2
macro_grid[1] *= 2
return results
def expected_results_fixed_h_refined_H_with_fine_PoU_in_basis(global_grid_refines, num_H_refs_less, initial_macro_grid):
if global_grid_refines == 3 and num_H_refs_less == 1 and initial_macro_grid == [2, 2]:
return {'num_subdomains': [1, 4, 16],
'num_DoFs': [256, 256, 256],
'error': [0.04924311575284581, 0.0492431157528458, 0.049243115752845845],
'estimate': [0.19600366888010526, 0.20314200704603694, 0.2167643126698387]}
else:
assert False, (f'missing expected results for global_grid_refines={global_grid_refines}, '
+ f'num_H_refs_less={num_H_refs_less}, initial_macro_grid={initial_macro_grid}')
def test_fixed_h_refined_H_with_fine_PoU_in_basis():
kwargs = {'global_grid_refines': 3, 'num_H_refs_less': 1, 'initial_macro_grid': [1, 1]}
assert_results(run_fixed_h_refined_H_with_fine_PoU_in_basis(**kwargs),
expected_results_fixed_h_refined_H_with_fine_PoU_in_basis(**kwargs))
if __name__ == '__main__':
r = run_fixed_H_refined_h()
for kk, vv in r.items():
......@@ -186,4 +290,9 @@ if __name__ == '__main__':
r = run_fixed_h_refined_H()
for kk, vv in r.items():
print(kk, vv)
print('')
print('')
r = run_fixed_h_refined_H_with_fine_PoU_in_basis()
for kk, vv in r.items():
print(kk, vv)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment